From 1a699b6bc0a707509ea5a54a4cc0750ec6f4a626 Mon Sep 17 00:00:00 2001 From: Anne_ctwtm Date: Tue, 27 Jan 2015 10:00:35 +0100 Subject: [PATCH] Going right --- CMakeLists.txt | 24 ++--- model.lyx | 277 +++++++++++++++++++++++++++++++++++++++++++++++++ src/solution.h | 54 ++++++++-- src/tube.cpp | 87 ++++++++++++++-- src/tube.h | 24 +++-- timedomain.cpp | 84 +++++++++++++-- 6 files changed, 499 insertions(+), 51 deletions(-) create mode 100644 model.lyx diff --git a/CMakeLists.txt b/CMakeLists.txt index a4a5497..f64fa2f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/model.lyx b/model.lyx new file mode 100644 index 0000000..bdba2df --- /dev/null +++ b/model.lyx @@ -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 diff --git a/src/solution.h b/src/solution.h index 5f5684c..9f216e9 100644 --- a/src/solution.h +++ b/src/solution.h @@ -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 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_ */ + diff --git a/src/tube.cpp b/src/tube.cpp index 58fa2a9..5183536 100644 --- a/src/tube.cpp +++ b/src/tube.cpp @@ -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< 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_ */ diff --git a/timedomain.cpp b/timedomain.cpp index 5782164..93f752b 100644 --- a/timedomain.cpp +++ b/timedomain.cpp @@ -1,17 +1,79 @@ -#include "geom.h" -#include "grid.h" #include "tube.h" +#include "tracer.h" + +#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=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; } + + +