Working code

This commit is contained in:
J.A. de Jong @ vulgaris 2015-02-22 11:36:23 +01:00
parent 55170494e3
commit bc4a8941b7
6 changed files with 104 additions and 175 deletions

106
main.cpp
View File

@ -1,106 +0,0 @@
#include "tube.h"
#include "tracer.h"
#include <sys/types.h>
#include <sys/stat.h>
#include <boost/program_options.hpp>
#include <ctime>
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<y?x:y;}
int main(int argc,char *argv[]){
d f=85.785;
int loglevel=20;
po::options_description desc("Allowed options");
desc.add_options()
("help", "Produce help message")
("frequency",po::value<d>(&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;
}

View File

@ -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()

View File

@ -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()

62
run.py Executable file
View File

@ -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<imax):
update_progress(i/imax)
sol=tube.getSol()
uarr[i]=sol.u()
parr[i]=sol.p()
# Update solution
t=sol.getTime()
i+=1
tube.DoIntegration(dt,intsteps)
print("")
print("Time integration done. Now saving data...")
save('udata.npy',uarr)
save('pdata.npy',parr)
print("Saving data done. Exiting.")

View File

@ -107,7 +107,7 @@ namespace td{
const vd& rho() const;
const vd& m() const;
const vd& rhoE() const;
d getTime() const;
void setTime(d t);
void setrho(d rho);
};

View File

@ -1,7 +1,7 @@
{
"metadata": {
"name": "",
"signature": "sha256:4b695d53d8ce9715401825830b83ae0060692d114cf9b5b8e49402b7bbd96c0a"
"signature": "sha256:6cbbbd8e5e8262619c7495c93758f9e0e660ffd345f9e5bb3915bf21352002f6"
},
"nbformat": 3,
"nbformat_minor": 0,
@ -12,8 +12,9 @@
"cell_type": "code",
"collapsed": false,
"input": [
"from timedomaineuler import *\n",
"#import timedomaineuler"
"from os import system\n",
"system('cmake . && make')\n",
"%run -i run.py"
],
"language": "python",
"metadata": {},
@ -23,18 +24,12 @@
"cell_type": "code",
"collapsed": false,
"input": [
"#Globalconf accessable with cvar.gc\n",
"f=85.785\n",
"T=1/f\n",
"loglevel=20\n",
"L=1.\n",
"gp=300\n",
"gc=cvar.gc #Reference!\n",
"dx=L/(gp-1); # One left and right gp, so\n",
"CFL=0.5;\n",
"gc.setfreq(f)\n",
"tube=TubeLF(L,gp)\n",
"dt=min(CFL*dx/gc.c0(),T/50)"
"uarr=load('udata.npy')\n",
"\n",
"uleft=zeros(uarr.shape[0],dtype=float)\n",
"for i in range(uarr.shape[0]):\n",
" uleft[i]=uarr[i][0]\n",
"plot(uleft)"
],
"language": "python",
"metadata": {},
@ -44,27 +39,13 @@
"cell_type": "code",
"collapsed": false,
"input": [
"nr_p_period=20\n",
"intsteps=int(floor(1./(gc.getfreq()*dt)/nr_p_period))"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"tube.DoIntegration(dt,intsteps)\n",
"sol=tube.getSol()\n",
"u=sol.u()\n",
"p=sol.p()\n",
"#rho=sol.rho()\n",
"figure(figsize=(9,3))\n",
"subplot(121)\n",
"plot(p)\n",
"subplot(122)\n",
"plot(u)"
"parr=load('pdata.npy')\n",
"\n",
"pend=zeros(parr.shape[0],dtype=float)\n",
"for i in range(parr.shape[0]):\n",
" pend[i]=parr[i][-1]\n",
"plot(pend)\n",
"show()"
],
"language": "python",
"metadata": {},