diff --git a/main.cpp b/main.cpp deleted file mode 100644 index 6a8491a..0000000 --- a/main.cpp +++ /dev/null @@ -1,106 +0,0 @@ -#include "tube.h" -#include "tracer.h" - -#include -#include -#include -#include -TRACETHIS; - -using geom::Geom; -using geom::Grid; -namespace po = boost::program_options; -using namespace td; -SPOILNAMESPACE; - -// Directory for storage -const char* dir="results"; - - -inline double min(d x,d y){return x(&f)->default_value(100)) - ; - d T=1/f; - d L=1; - us gp=500; - - d dx=L/(gp-1); // One left and right gp, so - d CFL=0.95; - // This is questionable - d dt=min(CFL*dx/gc.c0(),T/50); - cout << "dt:" << dt <<"\n"; - cout << "dx:" << dx <<"\n"; - - td::gc=tasystem::Globalconf::airSTP(0,f,1); - cout << "omg:" << td::gc.getomg() << endl; - inittrace(loglevel); - initothertrace(loglevel,nltracer); - - - // Create results dir, put all stuff in - mkdir(dir,S_IRWXU | S_IRWXG); - // Empty this directory - - td::Tube t(L,gp); - - - int nsaves_per_period=50; - - int approx_instances_per_period=round(T/dt); - if(nsaves_per_period>approx_instances_per_period){ // Not good.. - WARN("Error: to many saves for number of time instances available. Try adjusting nsaves_per_period Exiting.."); - exit(192); - } - - - d time=0; - - us loopnr=1; - us savenr=0; - us saves_per_period=20; - - d tbetweensaves=T/saves_per_period; - - d tlastsave=0; - - std::stringstream pfilename,ufilename; - pfilename << dir; ufilename << dir; - pfilename << "/p.dat"; - ufilename << "/u.dat"; - std::ofstream p(pfilename.str()); - std::ofstream u(ufilename.str()); - p << "#Pressure data vs gridpoint. gp=" << gp << "\n"; - u << "#Velocity data vs gridpoint. gp=" << gp << "\n"; - - t.getSol().save(p,data::p); - t.getSol().save(u,data::u); - const clock_t begin_time = std::clock(); - while(time<180*T){ - t.DoIntegration(dt); - time=t.getTime(); - - if(time-tlastsave>tbetweensaves){ - tlastsave=time; - savenr++; - t.getSol().save(p,data::p); - t.getSol().save(u,data::u); - cout << "Saving for t= " << time << " s.\n"; - } // ti%10 - loopnr++; - } - d comptime=std::clock()-begin_time; - p << "#Total time required to obtain solution: " << comptime/CLOCKS_PER_SEC <<" \n"; - return 0; -} - - - diff --git a/pmovie.py b/pmovie.py index 82e5474..a1eceb1 100755 --- a/pmovie.py +++ b/pmovie.py @@ -3,48 +3,53 @@ from numpy import * from matplotlib.pyplot import * from matplotlib import animation +from run import * import sys inst=0 -print("Loading data...") -pdat=loadtxt('results/p.dat') + +parr=load('pdata.npy') +uarr=load('udata.npy') # save(pdat,'p.npy') +imax=parr.shape[0] -gp=300 # Need to be taken from file -xmax=pdat[-1,0] -x=linspace(0,1,300) - +amplitude=1 fig = figure(figsize=(12,8)) -ax = axes(xlim=(0,1),ylim=(-180,180)) +ax = axes(xlim=(0,x[-1]),ylim=(-amplitude,amplitude)) grid('on') line, = plot([],[],lw=2) -ponly=pdat[:,1] -imax=int(floor(pdat.shape[0]/gp)) - -p=array_split(ponly,imax) - def init(): line.set_data([], []) return line, # print(imax) -def update_progress(progress): - sys.stdout.write('\r[{0}{1}] {2}%'.format('#'*int(floor(progress)),' '*int(ceil((100-progress))),progress)) -def animate(i): - update_progress(int(round(100*i/imax))) - line.set_data(x,p[i]) +def animatep(i): + # update_progress(int(round(100*i/imax))) + line.set_data(x,parr[i]) + return line, + +# animp = animation.FuncAnimation(fig, animatep, init_func=init, +# frames=imax, interval=20, blit=True,repeat=False) +def animateu(i): + # update_progress(int(round(100*i/imax))) + line.set_data(x,uarr[i]) return line, print("Creating animation...") -anim = animation.FuncAnimation(fig, animate, init_func=init, - frames=imax, interval=20, blit=True) +animu = animation.FuncAnimation(fig, animateu, init_func=init, + frames=imax, interval=20, blit=True,repeat=False) + # print("Saving animation...") # anim.save('p.mp4', fps=20, extra_args=['-vcodec', 'libx264']) print("Showing animation...") show() + + + + diff --git a/ptime.py b/ptime.py deleted file mode 100755 index 02e6920..0000000 --- a/ptime.py +++ /dev/null @@ -1,13 +0,0 @@ -#!/usr/bin/python - -from numpy import * -from matplotlib.pyplot import * - -pfile=loadtxt('results/p.dat') -gp=300 -pend=zeros(floor(pfile.shape[0]/gp)) -for i in range(pend.size): - pend[i]=pfile[(i+1)*gp-1,1] - -plot(pend) -show() diff --git a/run.py b/run.py new file mode 100755 index 0000000..97e179c --- /dev/null +++ b/run.py @@ -0,0 +1,62 @@ +#!/usr/bin/python + +from timedomaineuler import * +from numpy import * +import sys +#Globalconf accessable with cvar.gc +f=85.785 #Frequency of oscillation + + +L=1. #Length of tube +gp=300 #Number of gridpoints + +nr_p_period=50 #Number of save per oscillation period +CFL=0.5; # CFL number +periods=100 #Number of periods to compute + +T=1/f +gc=cvar.gc #Reference! +dx=L/(gp-1); # One left and right gp, so + +gc.setfreq(f) +tube=TubeLF(L,gp) +dt=min(CFL*dx/gc.c0(),T/nr_p_period) + +intsteps=int(floor(1./(gc.getfreq()*dt)/nr_p_period)) +x=linspace(0,L,gp) + +# Create tube instance +tube=TubeLF(L,gp) + +# To create a nice progress bar + +def update_progress(progress): + nbars=60 + text1="\r[{0}{1}] {2}%".format('#'*int(floor(nbars*progress)),' '*int(nbars-floor(nbars*progress)),"%0.0f" %(100*progress)) + # text=",' '*(80-int(80*floor(progress))),progress) + sys.stdout.write(text1) + sys.stdout.flush() +if(__name__=="__main__"): + + # Create data storage + uarr=zeros((nr_p_period*periods+1,gp),dtype=float) + parr=zeros((nr_p_period*periods+1,gp),dtype=float) + t=0.0 + i=0 + # Now we create the results + imax=nr_p_period*periods + print("Start time integration...") + while(i