ASCEE/timedomainresonace
5
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Going right

Browse Source
This commit is contained in:
Anne_ctwtm 2015-01-27 10:00:35 +01:00
parent a534e46833
commit 1a699b6bc0
6 changed files with 499 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
CMakeLists.txt
View File

@ -13,37 +13,31 @@ add_definitions(-DTRACERNAME=timedomaineuler -DARMA_USE_BLAS -DARMA_USE_LAPACK)
# read from a pipe; but the GNU assembler has no troubl.
# ${tubedrag}
set (CMAKE_GCC " -Wno-unused-function -ffunction-sections -fdata-sections -Wno-unused-local-typedefs -Wno-empty-body")
set (CMAKE_GENERAL "${CMAKE_GENERAL} -std=c++11 -pipe -fPIC -Wall \
set (CMAKE_GENERAL "${CMAKE_GENERAL} -std=c++11 -pipe -O2 -fPIC -Wall \
-Wextra -Wno-unused-parameter \
-Wno-unused-variable -Wno-unused-but-set-variable \
-Wno-return-local-addr -Wno-cpp -Wno-address")
set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CMAKE_GENERAL} ${CMAKE_GCC}")
link_directories(${link_directories} /usr/local/lib)
# link_directories(${link_directories} /usr/local/lib)
include_directories(
../tmtubes/src
../tmtubes/src/sys
../tmtubes/src/seg
../tmtubes/src/nonlinear/sys
../tmtubes/src/nonlinear/seg
../tmtubes/src/nonlinear/geom
../tmtubes/src/nonlinear/geom/grid
../tmtubes/src/common
../tmtubes/src/common/gas
../tmtubes/src/common/bessel
../tmtubes/src/common/fsolve
../tmtubes/src/common/solid
../tmtubes/src/common/rottfuncs
src
)
AUX_SOURCE_DIRECTORY(src src)
link_directories(tmtubes/src/nonlinear tmtubes/src/common)
# set(src ${src} src/seg/geom.cpp src/sys/globalconf.cpp)
link_directories(/home/anne/bin/lib)
add_executable(timedomaineuler timedomain.cpp ${src})
target_link_libraries(timedomaineuler nonlin math_common blas lapack boost_iostreams boost_serialization)
target_link_libraries(timedomaineuler nonlin armadillo boost_program_options)

277
model.lyx Normal file
View File

@ -0,0 +1,277 @@
#LyX 2.1 created this file. For more info see http://www.lyx.org/
\lyxformat 474
\begin_document
\begin_header
\textclass article
\use_default_options true
\maintain_unincluded_children false
\language english
\language_package default
\inputencoding auto
\fontencoding global
\font_roman default
\font_sans default
\font_typewriter default
\font_math auto
\font_default_family default
\use_non_tex_fonts false
\font_sc false
\font_osf false
\font_sf_scale 100
\font_tt_scale 100
\graphics default
\default_output_format default
\output_sync 0
\bibtex_command default
\index_command default
\paperfontsize default
\use_hyperref false
\papersize default
\use_geometry false
\use_package amsmath 1
\use_package amssymb 1
\use_package cancel 1
\use_package esint 1
\use_package mathdots 1
\use_package mathtools 1
\use_package mhchem 1
\use_package stackrel 1
\use_package stmaryrd 1
\use_package undertilde 1
\cite_engine basic
\cite_engine_type default
\biblio_style plain
\use_bibtopic false
\use_indices false
\paperorientation portrait
\suppress_date false
\justification true
\use_refstyle 1
\index Index
\shortcut idx
\color #008000
\end_index
\secnumdepth 3
\tocdepth 3
\paragraph_separation indent
\paragraph_indentation default
\quotes_language english
\papercolumns 1
\papersides 1
\paperpagestyle default
\tracking_changes false
\output_changes false
\html_math_output 0
\html_css_as_file 0
\html_be_strict false
\end_header
\begin_body
\begin_layout Title
1D Euler time domain
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{eqnarray*}
\frac{\partial\rho}{\partial t}+\frac{\partial m}{\partial x} & = & 0\\
\frac{\partial m}{\partial t}+\frac{\partial\left(\rho u^{2}+p_{0}+p\right)}{\partial x} & = & 0\\
\frac{\partial E}{\partial t}+\frac{\partial\left[\left(E+p\right)u\right]}{\partial x} & = & 0
\end{eqnarray*}
\end_inset
\end_layout
\begin_layout Standard
with
\begin_inset Note Note
status open
\begin_layout Plain Layout
\begin_inset Formula $\rho e=\rho c_{v}T=\frac{c_{v}}{R}\rho R_{s}T=\frac{c_{v}}{R}\left(p_{0}+p\right)=\frac{1}{\gamma-1}\left(p_{0}+p\right)$
\end_inset
\end_layout
\end_inset
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
E=\rho e+\rho\tfrac{1}{2}u^{2}=\frac{p_{0}+p}{\left(\gamma-1\right)}+\rho\tfrac{1}{2}u^{2}
\]
\end_inset
\end_layout
\begin_layout Standard
If we replace
\begin_inset Formula $E$
\end_inset
with
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
E=\hat{E}+\frac{p_{0}}{\gamma-1}
\end{equation}
\end_inset
then
\begin_inset Formula
\[
\frac{\partial\hat{E}}{\partial t}+\frac{\partial\left(\hat{E}+p_{0}+p\right)u}{\partial x}=0
\]
\end_inset
so
\begin_inset Formula $E$
\end_inset
is internal energy per unit volume!
\end_layout
\begin_layout Section
Scheme
\end_layout
\begin_layout Standard
If we say:
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\frac{\partial U}{\partial t}+\frac{\partial}{\partial x}\left(F(U)\right)=0
\]
\end_inset
which is
\begin_inset Formula
\begin{equation}
\frac{U_{i}^{n+1}-\left(\frac{1}{2}\left(U_{i+1}+U_{i-1}^{^{2}}\right)\right)}{\Delta t}+\frac{F(U_{i+1}^{n})-F(U_{i-1}^{n})}{2\Delta x}=0
\end{equation}
\end_inset
in its discrete form
\end_layout
\begin_layout Standard
then
\begin_inset Formula
\[
U=\left\{ \begin{array}{c}
\rho\\
\rho u\\
\rho E
\end{array}\right\}
\]
\end_inset
and Lax-Friedrichs says for a middle node
\begin_inset Formula
\[
U_{i}^{n+1}=\frac{1}{2}\left(U_{i+1}^{n}+U_{i-1}^{n}\right)-\lambda\left(F\left(U_{i+1}^{n}\right)-F(U_{i-1}^{n})\right)
\]
\end_inset
with
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\lambda=\frac{\Delta t}{2\Delta x}
\]
\end_inset
Total domain length:
\begin_inset Formula $L$
\end_inset
.
If number of gridpoints is 3, one left, one right, than
\begin_inset Formula $dx=L/(gp-1$
\end_inset
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
L
\]
\end_inset
\end_layout
\begin_layout Section
Right wall bc
\end_layout
\begin_layout Standard
At a wall:
\begin_inset Formula
\[
m=0
\]
\end_inset
and
\begin_inset Formula
\[
\frac{\partial}{\partial x}\left(F(U)\right)\approx\frac{F\left(U_{i}^{n}\right)-F(U_{i-1}^{n})}{\Delta x}
\]
\end_inset
but
\end_layout
\begin_layout Section
Left pressure bc
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\frac{\partial}{\partial x}\left(F(U)\right)\approx\frac{F\left(U_{i+1}^{n}\right)-F(U_{i}^{n})}{\Delta x}
\]
\end_inset
for continuity, energy for momentum:
\begin_inset Formula
\[
\frac{\partial}{\partial x}\left(F(U)\right)\approx\frac{\left[\rho u^{2}+p\right]_{i+1}-(\rho u^{2})|_{i}-P_{pres}}{\Delta x}
\]
\end_inset
\end_layout
\end_body
\end_document

54
src/solution.h
View File

@ -3,23 +3,59 @@
#define _SOLUTION_H_
#include "vtypes.h"
#include "geom.h"
#include "globalconf.h"
namespace td{
SPOILNAMESPACE
class Solution{
extern tasystem::Globalconf gc;
vd rho;
vd u;
vd e;
vd p;
class SolutionAtGp
{
d rho_,m_=0,rhoE_=0;
public:
Solution(){}
~Solution(){}
// Dependent variables
const d& rho() const {return rho_;}
const d& m() const {return m_;} // which is rho*u
const d& rhoE() const {return rhoE_;} // which is rho*E
// Derived from dependent
d u() const {return m()/rho();}
d ekin() const {return 0.5*rho()*pow(u(),2);}
d estat() const {return rhoE()-ekin();}
d p() const {return estat()*(gc.gas.gamma(gc.T0)-1);}
d Cflux() const {return m();}
d Mflux() const {return pow(m(),2)/rho()+p();}
d Eflux() const {return (rhoE()+gc.p0/(gc.gas.gamma(gc.T0)-1)+p()+gc.p0)*u();}
SolutionAtGp();
~SolutionAtGp(){}
void set(d rho,d m=0,d rhoE=0){
rho_=rho;
m_=m;
rhoE_=rhoE;
}
};
class SolutionInstance{
vector<SolutionAtGp> gps;
d time=0;
int nr=0;
public:
SolutionInstance(d TimeInstance,int gp);
~SolutionInstance(){}
d getTime() const {return time;}
void setTime(d t) {time=t;}
void set(us i,SolutionAtGp& sol){ gps.at(i)=sol;}
SolutionAtGp& get(us i){ return gps.at(i); }
void setrho(d rho);
int safe(const char* dir);
};
} // namespace td
#endif /* _SOLUTION_H_ */

87
src/tube.cpp
View File

@ -2,13 +2,88 @@
namespace td{
// Definition of a Globalconf in this namespace. Is initialized in main()
tasystem::Globalconf gc;
Tube::Tube(const segment::Geom& g):
geom(g),
gc(tasystem::Globalconf::airSTP(0,100,1))
{
d Tube::pleft(d t){
d pleft=std::sin(gc.getomg()*t);
// cout << "t:"<< t<<endl;
// cout << "pleft:"<< pleft<<endl;
return pleft;
}
Tube::Tube(d L,int gp): L(L),gp(gp),sol(SolutionInstance(0,gp)) {
TRACE(15,"Tube::Tube()");
dx=L/(gp-1);
VARTRACE(15,gp);
sol.setrho(gc.rho0());
}
};
void Tube::DoIntegration(d dt){
TRACE(14,"Tube::DoIntegration()");
// Define new time
d newt=t+dt;
// Define new solutioninstance
SolutionInstance newsol=sol;
SolutionAtGp newsolgp;
d rho,m,rhoE;
newsol.setTime(newt);
int i=0;
// Left boundary
{
SolutionAtGp& ip1=sol.get(i+1);
SolutionAtGp& ip0=sol.get(i);
d la=dt/(dx);
rho=ip0.rho()-la*(ip1.Cflux()-ip0.Cflux());
d momfluxl=ip0.rho()*pow(ip0.u(),2)+pleft(t);
m=ip0.m()-la*(ip1.Mflux()-momfluxl);
rhoE=ip0.rhoE()-la*(ip1.Eflux()-ip0.Eflux());
// rhoE=ip0.rhoE()-la*(ip1.Eflux()-.rhoE()*ip0.u()-ip0.u()*(gc.p0+pleft(t)));
newsolgp.set(rho,m,rhoE);
// cout << "pleft: " << pleft(t) <<"\n";
// cout << "pleft: " << newsolgp.p() <<"\n";
newsol.set(i,newsolgp);
} // Leftmost node
// TRACE(15,"leftmost done");
// Inner nodes
for(i=1;i<gp-2;i++){
SolutionAtGp& im1=sol.get(i-1);
SolutionAtGp& ip1=sol.get(i+1);
d lambda=dt/(2*dx);
rho=0.5*(im1.rho()+ip1.rho())-lambda*(ip1.Cflux()-im1.Cflux());
m=0.5*(im1.m()+ip1.m())-lambda*(ip1.Mflux()-im1.Mflux());
rhoE=0.5*(im1.rhoE()+ip1.rhoE())-lambda*(ip1.Eflux()-im1.Eflux());
newsolgp.set(rho,m,rhoE);
newsol.set(i,newsolgp);
} // for over all gridpoints in mid
{
i++;
SolutionAtGp& im1=sol.get(i-1);
SolutionAtGp& ip0=sol.get(i);
d la=dt/(dx);
rho=ip0.rho()-la*(0-im1.Cflux());
m=0;
rhoE=ip0.rhoE()-la*(0-im1.Eflux());
newsolgp.set(rho,m,rhoE);
newsol.set(i,newsolgp);
}
// Finally, update time and solution
sol=newsol;
t=newt;
}
} // namespace td

24
src/tube.h
View File

@ -6,21 +6,25 @@
#include "globalconf.h"
#include "solution.h"
namespace td{
extern tasystem::Globalconf gc;
class Tube{
segment::Geom geom;
tasystem::Globalconf gc;
vector<Solution> sols;
d dx,L; // Grid spacing, total length
int gp; // Number of gridpoints
SolutionInstance sol; // Solutions at time instances
d t=0; // Current time
d pleft(d t); // Compute pressure bc
public:
Tube(const segment::Geom& g);
Solution& getSol(int i=-1);
};
Tube(d L,int gp);
SolutionInstance& getSol() { return sol;}
void setSol(const SolutionInstance& sol) {this->sol=sol;}
void DoIntegration(d dt);
d getTime(){return t;}
};
} // namespace td
#endif /* _TUBE_H_ */

84
timedomain.cpp
View File

@ -1,17 +1,79 @@
#include "geom.h"
#include "grid.h"
#include "tube.h"
#include "tracer.h"
#include <sys/types.h>
#include <sys/stat.h>
#include <boost/program_options.hpp>
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;}
using segment::Geom;
using segment::Grid;
int main(int argc,char *argv[]){
SPOILNAMESPACE
us gp=100;
d L=1;
Grid grid(gp,L);
d radius=1.0;
Geom g=Geom::CylinderBlApprox(grid,radius);
td::Tube t(g);
d f=100;
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=0.01;
us gp=10;
d dx=L/gp;
d CFL=0.5;
// 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(15);
initothertrace(10,nltracer);
// Create results dir, put all stuff in
// Create a directory for results
mkdir(dir,S_IRWXU | S_IRWXG);
// Empty this directory
system("rm results/*");
td::Tube t(L,gp);
int nsaves_per_period=50;
int nr=0;
d time=0;
while(time<2*T){
t.DoIntegration(dt);
time=t.getTime();
// cout << "Iteration done. Saving data...\n";
if(true){
// Create filename to store data to
} // ti%10
}
return 0;
}