From 0de55ea7fd56ae7696db9eaa2292f45464a2445d Mon Sep 17 00:00:00 2001 From: "J.A. de Jong @ vulgaris" Date: Fri, 2 Dec 2016 22:47:07 +0100 Subject: [PATCH] Tested 1D-solver, added adiabatictemp, which solves the adiabatic temperature change due to a pressure change --- CMakeLists.txt | 1 + doc/usg.lyx | 4 +- src/CMakeLists.txt | 48 +++++++------ src/funcs/CMakeLists.txt | 7 -- src/material/CMakeLists.txt | 7 -- src/material/adiabatictemp.cpp | 122 +++++++++++++++++++++++++++++++++ src/material/adiabatictemp.h | 19 +++++ src/material/air.h | 2 - src/solver/CMakeLists.txt | 7 -- src/solver/brent.cpp | 37 +++++----- src/solver/solver.cpp | 35 ++++++---- src/solver/solver.h | 2 +- src/sys/CMakeLists.txt | 10 --- src/sys/jaccol.cpp | 2 +- src/sys/tasmet_variable.cpp | 1 + src/sys/tasmet_variable.h | 5 +- testing/solver/CMakeLists.txt | 4 +- testing/solver/brent_test.cpp | 6 ++ 18 files changed, 225 insertions(+), 94 deletions(-) delete mode 100644 src/funcs/CMakeLists.txt delete mode 100644 src/material/CMakeLists.txt create mode 100644 src/material/adiabatictemp.cpp create mode 100644 src/material/adiabatictemp.h delete mode 100644 src/solver/CMakeLists.txt delete mode 100644 src/sys/CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt index a6f0183..434c862 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -105,6 +105,7 @@ include_directories( src src/solver src/material + src/sys ) # Add the code subdirectory add_subdirectory(src) diff --git a/doc/usg.lyx b/doc/usg.lyx index 6928247..cebbed8 100644 --- a/doc/usg.lyx +++ b/doc/usg.lyx @@ -1881,7 +1881,6 @@ Create a PhaseConstraint object: \end_layout \begin_layout Verbatim - pc=PhaseConstraint(Varnr, freqnr, left) \end_layout @@ -1890,7 +1889,6 @@ Apply this contraint to a segment which can accept them, for example a Tube: \end_layout \begin_layout Verbatim - t1.setPhaseConstraint(pc) \end_layout @@ -3817,7 +3815,7 @@ In that case \begin_layout Standard \begin_inset Formula \begin{equation} -\frac{1}{R_{s}}\left(c_{p,0}\ln\left(\frac{T_{p}}{T_{0}}\right)+\sum_{i=1}^{N_{c_{p}}}\frac{c_{p,i}\left(T_{p}-T_{0}\right)^{i}}{i}\right)-\ln\left(\frac{p_{p}-p_{0}}{p_{0}}\right)=0.\label{eq:T-p-adiabatic} +\frac{1}{R_{s}}\left(c_{p,0}\ln\left(\frac{T_{p}}{T_{0}}\right)+\sum_{i=1}^{N_{c_{p}}}\frac{c_{p,i}\left(T_{p}^{i}-T_{0}^{i}\right)}{i}\right)-\ln\left(\frac{p_{p}-p_{0}}{p_{0}}\right)=0.\label{eq:T-p-adiabatic} \end{equation} \end_inset diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9cb5666..b8f55d1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -24,30 +24,38 @@ include_directories( # volume ) -add_subdirectory(funcs) -add_subdirectory(material) -# add_subdirectory(duct) -# add_subdirectory(mech) -# add_subdirectory(seg) -add_subdirectory(solver) -add_subdirectory(sys) -# add_subdirectory(var) - -add_library(tasmet_src_src +add_library(tasmet_src tasmet_tracer.cpp tasmet_exception.cpp tasmet_assert.cpp - ) -add_library(tasmet_src) -target_link_libraries(tasmet_src - funcs - material - solver - sys - # This one should be last as other parts link to these utilities - tasmet_src_src - ) + funcs/bessel.cpp + funcs/cbessj.cpp + funcs/rottfuncs.cpp + funcs/skewsine.cpp + funcs/sph_j.cpp + + solver/solver.cpp + solver/system.cpp + solver/newton_raphson.cpp + solver/brent.cpp + + material/gas.cpp + material/air.cpp + material/helium.cpp + material/nitrogen.cpp + material/solid.cpp + material/adiabatictemp.cpp + + # sys/enginesystem.cpp + sys/globalconf.cpp + sys/tasmet_variable.cpp + sys/jaccol.cpp + sys/jacobian.cpp + sys/jacrow.cpp + # sys/tasystem.cpp + sys/triplets.cpp +) # set_source_files_properties(swig/nonlin.i # PROPERTIES CPLUSPLUS ON) diff --git a/src/funcs/CMakeLists.txt b/src/funcs/CMakeLists.txt deleted file mode 100644 index 1c75ddb..0000000 --- a/src/funcs/CMakeLists.txt +++ /dev/null @@ -1,7 +0,0 @@ -add_library(funcs - bessel.cpp - cbessj.cpp - rottfuncs.cpp - skewsine.cpp - sph_j.cpp - ) diff --git a/src/material/CMakeLists.txt b/src/material/CMakeLists.txt deleted file mode 100644 index 96e6462..0000000 --- a/src/material/CMakeLists.txt +++ /dev/null @@ -1,7 +0,0 @@ -add_library(material - gas.cpp - air.cpp - helium.cpp - nitrogen.cpp - solid.cpp - ) diff --git a/src/material/adiabatictemp.cpp b/src/material/adiabatictemp.cpp new file mode 100644 index 0000000..ac8aeca --- /dev/null +++ b/src/material/adiabatictemp.cpp @@ -0,0 +1,122 @@ +// adiabatictemp.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#define TRACERPLUS 10 +#include "adiabatictemp.h" +#include "perfectgas.h" +#include "tasmet_variable.h" +#include "tasmet_tracer.h" +#include "tasmet_exception.h" +#include "brent.h" + +struct AdiabaticTemp : public NoGradientNonlinearSystem +{ + d p,p0,Rs,T0; + d T; + const PerfectGas& gas; // Careful here, reference! + const vd& cpc; // Here too! + + d gamma0; // Used for initial guess + + AdiabaticTemp(const PerfectGas& gas,d p): + gas(gas), + cpc(gas.cpc()) + { + + p0 = gas.p0(); + T0 = gas.T0(); + Rs = gas.Rs(); + + gamma0 = gas.gamma(T0,p0); + + #ifdef TASMET_DEBUG + if(cpc.size() < 1) + throw TaSMETError("Invalid array of heat capacity coefficients"); + #endif + + setP(p); + + } + + d getSolution() const { return T; } + void updateSolution(const d& T) { this->T=T;} + + AdiabaticTemp* copy() const { return new AdiabaticTemp(gas,p);} + + // Set a pressure to solve the temperature for + void setP(d p) { + + this->p = p; + + // Good initial guess + T = T0*pow(p/p0,(gamma0-1)/gamma0); + + } + + + + // Final function to solve for + d residual() const { + + // Do the integration for the temperature + d lhs = (cpc(0)/Rs)*log(T/T0); + for(us i=1;i cb = solver_callback; + + for(us i=0;i& system, d fa = system.residual(a); d fb = system.residual(b); - + if((fb) == 0) { + TRACE(15,"Found root during bracketing"); + return; + } d c = a; d fc = fa; @@ -85,7 +88,7 @@ void Brent::start_implementation(NoGradientNonlinearSystem& system, if((fa!=fc) && (fb!=fc)){ // Inverse quadratic interpolation - TRACE(15,"IQI"); + TRACE(5,"IQI"); t = a*fb*fc/((fa-fb)*(fa-fc)); u = b*fa*fc/((fb-fa)*(fb-fc)); v = c*fa*fb/((fc-fa)*(fc-fb)); @@ -94,7 +97,7 @@ void Brent::start_implementation(NoGradientNonlinearSystem& system, } else { // Secant method - TRACE(15,"Secant"); + TRACE(5,"Secant"); s = b-fb*(b-a)/(fb-fa); } @@ -104,19 +107,19 @@ void Brent::start_implementation(NoGradientNonlinearSystem& system, d absbmc = abs(b-c); d abscmd = abs(c-d_); - VARTRACE(15,s); + VARTRACE(5,s); if((bisec_flag |= (!is_between(s,(3*a+b)/4,b)))) goto bflag; - TRACE(15,"Survived 1"); + TRACE(5,"Survived 1"); if(bisec_flag |= (mflag && (abssmb >= absbmc/2))) goto bflag; - TRACE(15,"Survived 2"); + TRACE(5,"Survived 2"); if(bisec_flag |= ((!mflag) && (abssmb >= abscmd/2))) goto bflag; - TRACE(15,"Survived 3"); + TRACE(5,"Survived 3"); if(bisec_flag |= (mflag && (absbmc < abs(_reltol)))) goto bflag;; - TRACE(15,"Survived 4"); + TRACE(5,"Survived 4"); bflag: if(bisec_flag || ((!mflag) && (abscmd < abs(_reltol)))) { - TRACE(15,"Bisection"); + TRACE(5,"Bisection"); s = (a+b)/2; mflag = true; } @@ -148,14 +151,14 @@ void Brent::start_implementation(NoGradientNonlinearSystem& system, progress.rel_err = abs(b-a); progress.iteration++; - VARTRACE(15,s); - VARTRACE(15,a); - VARTRACE(15,b); - VARTRACE(15,c); - VARTRACE(15,fa); - VARTRACE(15,fb); - VARTRACE(15,fc); - VARTRACE(15,fs); + VARTRACE(5,s); + VARTRACE(5,a); + VARTRACE(5,b); + VARTRACE(5,c); + VARTRACE(5,fa); + VARTRACE(5,fb); + VARTRACE(5,fc); + VARTRACE(5,fs); SolverAction action = (*callback)(progress); diff --git a/src/solver/solver.cpp b/src/solver/solver.cpp index 267c0ba..d4e61e3 100644 --- a/src/solver/solver.cpp +++ b/src/solver/solver.cpp @@ -5,11 +5,12 @@ // Description: // ////////////////////////////////////////////////////////////////////// - +#include "tasmet_tracer.h" #include "solver.h" #include "tasmet_exception.h" #include + template static void SolverThread(Solver* solver, system_T* system, @@ -44,19 +45,19 @@ void Solver::start(progress_callback* callback,bool wait){ assert(_solver_thread == nullptr); - this->_solver_thread = new std::thread(SolverThread, - this, - _sys, - callback); - if(!_solver_thread) - throw TasMETBadAlloc(); + if(!wait) { + this->_solver_thread = new std::thread(SolverThread, + this, + _sys, + callback); + if(!_solver_thread) + throw TasMETBadAlloc(); + } + else { + + TRACE(15,"Waiting for solver..."); + start_implementation(*_sys,callback); - if(wait){ - while (_running){ - std::this_thread::sleep_for(std::chrono::milliseconds(1)); - } - // Cleanup resources - stop(); } } @@ -64,6 +65,7 @@ template void Solver::stop() { _running = false; + if(_solver_thread){ _solver_thread->join(); @@ -73,12 +75,15 @@ void Solver::stop() { } } template -result_T Solver::getSolution() const{ +result_T Solver::getSolution() { + if(_running){ throw TaSMETError("Solver is running"); } + // Cleanup thread resources stop(); + return _sys->getSolution(); } template @@ -94,7 +99,7 @@ void SolverThread(Solver* solver,system_T* system,progress_ca } -// Explicit instantiation for a +// Explicit instantiation for some types of systems and results template class Solver,vd>; template class Solver; template class Solver,d>; diff --git a/src/solver/solver.h b/src/solver/solver.h index 25085af..c668387 100644 --- a/src/solver/solver.h +++ b/src/solver/solver.h @@ -49,7 +49,7 @@ public: void stop(); // Stops the solver // Returns the solution of the problem - result_T getSolution() const; + result_T getSolution(); virtual ~Solver(); template diff --git a/src/sys/CMakeLists.txt b/src/sys/CMakeLists.txt deleted file mode 100644 index cc24fc7..0000000 --- a/src/sys/CMakeLists.txt +++ /dev/null @@ -1,10 +0,0 @@ -add_library(sys - # enginesystem.cpp - globalconf.cpp - tasmet_variable.cpp - jaccol.cpp - jacobian.cpp - jacrow.cpp - # tasystem.cpp - triplets.cpp - ) diff --git a/src/sys/jaccol.cpp b/src/sys/jaccol.cpp index 7c4915e..6202b88 100644 --- a/src/sys/jaccol.cpp +++ b/src/sys/jaccol.cpp @@ -12,7 +12,7 @@ JacCol::JacCol(const Variable& thevar): coldof_(thevar.getDofNr()), - data_(thevar.getGc().Ns(),thevar.getGc().Ns(),fillwith::zeros) + data_(thevar.getGc()->Ns(),thevar.getGc()->Ns(),fillwith::zeros) { } JacCol::JacCol(us coldof,const GlobalConf& gc): coldof_(coldof), diff --git a/src/sys/tasmet_variable.cpp b/src/sys/tasmet_variable.cpp index 56f3660..49a7d8e 100644 --- a/src/sys/tasmet_variable.cpp +++ b/src/sys/tasmet_variable.cpp @@ -44,6 +44,7 @@ Variable Variable::operator*(const Variable& Variable2) const { // Multiply two return Variable(this->_gc,_tdata%Variable2._tdata,false); } //***************************************** The Variable class +Variable::Variable(const gc_ptr& gc): Variable(gc,0){} Variable::Variable(const gc_ptr& gc,double initval): _gc(gc) { TRACE(0,"Variable::Variable(us ndofs, double initval)"); diff --git a/src/sys/tasmet_variable.h b/src/sys/tasmet_variable.h index 7a23923..7000c62 100644 --- a/src/sys/tasmet_variable.h +++ b/src/sys/tasmet_variable.h @@ -45,7 +45,8 @@ public: ~Variable(){} void setGc(const gc_ptr&); - const GlobalConf& getGc() const {return *_gc;} + const gc_ptr& getGc() const {return _gc;} + void updateNf(); // Get methods @@ -56,7 +57,7 @@ public: // Extract data const vd& tdata() const {return _tdata; } //Get time data - const vd& adata() const {return _adata; } //Get time data //vector + const vd& adata() const {return _adata; } //Get amplitude data //vector // Obtain a time response vector vd timeResponse(us nperiod=2,us ninst=100) const; diff --git a/testing/solver/CMakeLists.txt b/testing/solver/CMakeLists.txt index e81a9f9..cc38c09 100644 --- a/testing/solver/CMakeLists.txt +++ b/testing/solver/CMakeLists.txt @@ -1,4 +1,4 @@ -add_executable(newton_raphson_test newton_raphson_test.cpp) +# add_executable(newton_raphson_test newton_raphson_test.cpp) add_executable(brent_test brent_test.cpp) target_link_libraries(brent_test tasmet_src armadillo openblas pthread) -target_link_libraries(newton_raphson_test tasmet_src armadillo openblas pthread) +# target_link_libraries(newton_raphson_test tasmet_src tasmet_src armadillo openblas pthread) diff --git a/testing/solver/brent_test.cpp b/testing/solver/brent_test.cpp index 2deaab7..57f5fce 100644 --- a/testing/solver/brent_test.cpp +++ b/testing/solver/brent_test.cpp @@ -12,6 +12,9 @@ #include "helium.h" #include "nitrogen.h" #include "air.h" +#include "globalconf.h" +#include "tasmet_variable.h" +#include "adiabatictemp.h" SolverAction solver_callback(SolverProgress p){ @@ -62,6 +65,9 @@ int main(){ Air air(293.15,101325); + gc_ptr gc(new GlobalConf(10,100)); + Variable pressure(gc,10*101325); + Variable temperature = adiabaticTemp(air,pressure); }