Tested 1D-solver, added adiabatictemp, which solves the adiabatic temperature change due to a pressure change

CMakeLists.txt

@@ -105,6 +105,7 @@ include_directories(
   src
   src/solver
   src/material
+  src/sys
   )
 # Add the code subdirectory
 add_subdirectory(src)

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

src/CMakeLists.txt

@@ -24,30 +24,38 @@ include_directories(
   # volume
   )
 
-add_subdirectory(funcs)
+add_library(tasmet_src
-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
   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)

src/funcs/CMakeLists.txt

@@ -1,7 +0,0 @@
-add_library(funcs
-  bessel.cpp
-  cbessj.cpp
-  rottfuncs.cpp
-  skewsine.cpp
-  sph_j.cpp
-  )

src/material/CMakeLists.txt

@@ -1,7 +0,0 @@
-add_library(material
-  gas.cpp
-  air.cpp
-  helium.cpp
-  nitrogen.cpp
-  solid.cpp
-  )

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> 
+{
+    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<cpc.size();i++){
+            lhs+= cpc(i)*(pow(T,d(i))-pow(T0,d(i)))/(i*Rs);
+        }
+
+        d rhs = log(p/p0);
+
+        return lhs-rhs;
+    }
+
+};
+namespace {
+    SolverAction solver_callback(SolverProgress p){
+
+        cout << "Iteration: " << p.iteration << ". Function error: " << p.fun_err <<". Relative error: " << p.rel_err << endl;
+        if(p.fun_err < 1e-7){
+            TRACE(9, "Returning stop");
+            return Stop;
+        }
+        else {
+            TRACE(9, "Returning continue");
+            return Continue;
+        }
+    }
+}
+Variable adiabaticTemp(const PerfectGas& gas,
+                       const Variable& pressure) {
+                       
+    
+    TRACE(10,"adiabaticTemp()");
+    Variable temp(pressure);
+    const vd p = pressure.tdata();
+
+    vd T = temp.tdata();
+
+    AdiabaticTemp adt(gas,p(0)); // Equation
+
+    std::function<SolverAction(SolverProgress)> cb = solver_callback;
+
+    for(us i=0;i<p.size();i++){
+
+        adt.setP(p(i));
+        Brent brent(adt);           // Solver
+
+        VARTRACE(15,p(i));
+        VARTRACE(15,adt.getSolution());
+        brent.start(&cb);
+
+        T(i) = brent.getSolution();
+        VARTRACE(15,T(i));
+
+    }
+    temp.settdata(T);
+    return temp;
+}
+//////////////////////////////////////////////////////////////////////
+

src/material/adiabatictemp.h

@@ -0,0 +1,19 @@
+// adiabatictemp.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Compute the adiabatic/isentropic temperature as a function of the pressure
+// for a thermally perfect gas
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef ADIABATICTEMP_H
+#define ADIABATICTEMP_H
+
+class PerfectGas;
+class Variable;
+
+Variable adiabaticTemp(const PerfectGas& gas,const Variable& pressure);
+
+#endif // ADIABATICTEMP_H
+//////////////////////////////////////////////////////////////////////

src/material/air.h

@@ -28,8 +28,6 @@ protected:
 public:
     Air(d T0,d p0):PerfectGas(air,T0,p0){}
     const vd& cpc() const {return _cpc;}
-    d cp(d T,d p) const;
-    d h(d T,d p) const;
     d mu(d T,d p) const;
     d kappa(d T,d p) const;
     ~Air(){}

src/solver/CMakeLists.txt

@@ -1,7 +0,0 @@
-include_directories(${PROJECT_SOURCE_DIR}/src/solver)
-add_library(solver
-  solver.cpp
-  system.cpp
-  newton_raphson.cpp
-  brent.cpp
-  )

src/solver/brent.cpp

@@ -57,7 +57,10 @@ void Brent::start_implementation(NoGradientNonlinearSystem<d>& 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<d>& 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<d>& 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<d>& 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<d>& system,
         progress.rel_err = abs(b-a);
         progress.iteration++;
 
-        VARTRACE(15,s);
+        VARTRACE(5,s);
-        VARTRACE(15,a);
+        VARTRACE(5,a);
-        VARTRACE(15,b); 
+        VARTRACE(5,b); 
-        VARTRACE(15,c);
+        VARTRACE(5,c);
-        VARTRACE(15,fa);
+        VARTRACE(5,fa);
-        VARTRACE(15,fb); 
+        VARTRACE(5,fb); 
-        VARTRACE(15,fc);
+        VARTRACE(5,fc);
-        VARTRACE(15,fs);
+        VARTRACE(5,fs);
 
         SolverAction action = (*callback)(progress);
 

src/solver/solver.cpp

@@ -5,11 +5,12 @@
 // Description:
 //
 //////////////////////////////////////////////////////////////////////
-
+#include "tasmet_tracer.h"
 #include "solver.h"
 #include "tasmet_exception.h"
 #include <chrono>
 
+
 template<typename system_T,typename result_T>
 static void SolverThread(Solver<system_T,result_T>* solver,
                          system_T* system,
@@ -44,19 +45,19 @@ void Solver<system_T,result_T>::start(progress_callback* callback,bool wait){
 
     assert(_solver_thread == nullptr);
 
-    this->_solver_thread = new std::thread(SolverThread<system_T,result_T>,
+    if(!wait) {
-                                           this,
+        this->_solver_thread = new std::thread(SolverThread<system_T,result_T>,
-                                           _sys,
+                                               this,
-                                           callback);
+                                               _sys,
-    if(!_solver_thread)
+                                               callback);
-        throw TasMETBadAlloc();
+        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<typename system_T,typename result_T>
 void Solver<system_T,result_T>::stop() {
 
     _running = false;
+
     if(_solver_thread){
 
         _solver_thread->join();
@@ -73,12 +75,15 @@ void Solver<system_T,result_T>::stop() {
     }
 }
 template<typename system_T,typename result_T>
-result_T Solver<system_T,result_T>::getSolution() const{
+result_T Solver<system_T,result_T>::getSolution() {
+
     if(_running){
         throw TaSMETError("Solver is running");
     }
+
     // Cleanup thread resources
     stop();
+
     return _sys->getSolution();
 }
 template<typename system_T,typename result_T>
@@ -94,7 +99,7 @@ void SolverThread(Solver<system_T,result_T>* solver,system_T* system,progress_ca
 }
 
 
-// Explicit instantiation for a
+// Explicit instantiation for some types of systems and results
 template class Solver<NoGradientNonlinearSystem<vd>,vd>;
 template class Solver<GradientNonlinearSystem,vd>;
 template class Solver<NoGradientNonlinearSystem<d>,d>;

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<typename Y,typename rT>

src/sys/CMakeLists.txt

@@ -1,10 +0,0 @@
-add_library(sys
-  # enginesystem.cpp
-  globalconf.cpp
-  tasmet_variable.cpp
-  jaccol.cpp
-  jacobian.cpp
-  jacrow.cpp
-  # tasystem.cpp
-  triplets.cpp
-  )

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),

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)");

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;

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)

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);
 }