Model residual is calculatable. Now Jacobian

CMakeLists.txt

@@ -91,22 +91,6 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/cmake_tools)
 # Python #####################
 # ##########################
 
-include_directories(/usr/include/PythonQt)
-
-
-if(TaSMET_PY_VERSION)
-  # Find major version from version string
-  set(PYTHON_LIBRARY "/usr/lib/libpython${TaSMET_PY_VERSION}.so")
-  set(PYTHON_INCLUDE_DIR "/usr/include/python${TaSMET_PY_VERSION}")  
-  set(PYTHON_INCLUDE_DIRS "/usr/include/python${TaSMET_PY_VERSION}")  
-
-endif(TaSMET_PY_VERSION)
-message("Python include dirs: ${PYTHON_INCLUDE_DIRS}")
-
-find_package(PythonLibs REQUIRED)
-string(REGEX MATCH "^." TaSMET_PY_MAJOR_VERSION ${PYTHONLIBS_VERSION_STRING})
-MESSAGE("Python major version: ${TaSMET_PY_MAJOR_VERSION}")  
-
 # Find the site_packages directory of python
 execute_process(COMMAND python${TaSMET_PY_MAJOR_VERSION} -c "from distutils.sysconfig import get_python_lib; print(get_python_lib())" OUTPUT_VARIABLE PYTHON_SITE_PACKAGES OUTPUT_STRIP_TRAILING_WHITESPACE)
 
@@ -120,7 +104,6 @@ add_definitions(-DARMA_USE_SUPERLU -DARMA_USE_CXX11)
 # link_directories(common)
 
 include_directories(
-  ${PYTHON_INCLUDE_DIRS}
   src
   src/solver
   src/material
@@ -135,5 +118,5 @@ add_subdirectory(testing)
 
 add_executable(tasmet src/main.cpp src/gui/tasmet_resources.qrc)
 add_executable(tasmet_solvemodel src/tasmet_solvemodel.cpp)
-target_link_libraries(tasmet tasmet_gui tasmet_src messages PythonQt Qt5::Widgets chaiscript_stdlib-5.7.0 openblas)
+target_link_libraries(tasmet tasmet_gui tasmet_src messages Qt5::Widgets chaiscript_stdlib-5.7.0 openblas)
-target_link_libraries(tasmet_solvemodel tasmet_src messages PythonQt openblas chaiscript_stdlib-5.7.0)
+target_link_libraries(tasmet_solvemodel tasmet_src messages openblas chaiscript_stdlib-5.7.0)

src/duct/duct.cpp

@@ -9,7 +9,7 @@
 #include "tasystem.h"
 #include "tasmet_assert.h"
 #include "tasmet_evalscript.h"
-
+#include "perfectgas.h"
 
 Duct::Duct(const us id,const pb::Duct& ductpb):
     Segment(id,ductpb.name()),
@@ -62,11 +62,127 @@ Duct::~Duct() {
 void Duct::residual(const TaSystem& sys,arma::subview_col<d> && residual) const {
 
     TRACE(15,"Duct::residual()");
+    
     const arma::subview_col<d> sol = sys.getSolution(_id);
 
-    VARTRACE(15,sol(0));
+    vd rho,u,T,p,Ts;            // Solution at this gp
+    vd rhop,up,Tp,pp,Tsp;       // Solution at next gp
 
+    // Continuity eq residual, momentum, energy, state, solid energy
+    vd cont,mom,en,st,sen;
     
+    // When we have to solve a solid heat balance
+    bool has_solideq = _ductpb.stempmodel() == pb::HeatBalance;
+    
+    us number_eqs = 4;
+    number_eqs += (has_solideq) ? 1 : 0;
+
+    VARTRACE(15,number_eqs);
+
+    us Ns = sys.Ns();
+    us eq_offset = 0;           // Equation offset for current gp
+    us res_offset = 0;          // Residual offset for current gp
+    us res_offsetp = 0;         // Residual offset for next gp
+
+    us gp_jump = number_eqs * Ns; // The jump per gp
+    
+    rhop = getvart(sys,constants::rho,0);
+    up = getvart(sys,constants::u,0);
+    Tp = getvart(sys,constants::T,0);
+    pp   = getvart(sys,constants::p,0);
+
+    const Gas& gas = sys.gas();
+
+    for(us gp=0;gp<ngp()-1;gp++) {
+
+        eq_offset = gp*Ns*number_eqs;
+        res_offset = eq_offset;
+        res_offsetp = res_offset + gp_jump;
+        
+        d dx = x(gp+1)-x(gp);
+
+        // Update the current gp solution
+        rho = rhop; u=up; T=Tp; p = pp; Ts=Tsp;
+
+        // Update the next gp solution
+        rhop = getvart(sys,constants::rho,gp+1);
+        up = getvart(sys,constants::u,gp+1);
+        Tp = getvart(sys,constants::T,gp+1);
+        pp   = getvart(sys,constants::p,gp+1);
+
+        cont  = ((rhop%up)-(rho%u))/dx;
+
+        mom = (rhop%up%up - rho%u%u + pp - p)/dx;
+        
+        switch (_ductpb.htmodel()) {
+        case pb::Isentropic: {
+
+           d T0 = gas.T0();
+           d p0 = gas.p0();
+           d rho0 = gas.rho0();
+           d gamma0 = gas.gamma(T0,p0);
+
+           en = p/p0 - pow(rho/rho0,gamma0);
+
+        }
+            break;
+        default:
+            tasmet_assert(false,"Not implemented htmodel");
+        }
+
+        st = gas.rho(T,p) - rho;
+
+        residual.subvec(eq_offset+0*Ns,eq_offset+1*Ns-1) = cont;
+        residual.subvec(eq_offset+1*Ns,eq_offset+2*Ns-1) = mom;
+        residual.subvec(eq_offset+2*Ns,eq_offset+3*Ns-1) = en;
+        residual.subvec(eq_offset+3*Ns,eq_offset+4*Ns-1) = st;
+
+    }
+
+    eq_offset += number_eqs*Ns;
+
+    // Equation of state for the last node
+    st = gas.rho(Tp,pp) - rhop;
+    residual.subvec(eq_offset,eq_offset+Ns-1) = st;
+
+    // Two more equations for the last grid point in case
+    // the heat transfer model is not a transport equation.
+    if(_ductpb.htmodel() == pb::Isentropic) {
+
+        eq_offset += Ns;
+
+        d T0 = gas.T0();
+        d p0 = gas.p0();
+        d rho0 = gas.rho0();
+        d gamma0 = gas.gamma(T0,p0);
+
+        en = p/p0 - pow(rho/rho0,gamma0);
+
+        residual.subvec(eq_offset,eq_offset+Ns-1) = en;
+    }
+
+}
+vd Duct::getvart(const TaSystem& sys,int varnr,int gp) const {
+    TRACE(15,"Duct::getvart()");
+    const arma::subview_col<d> sol = sys.getSolution(_id);
+
+    us Ns = sys.Ns();
+
+    // Wraparound
+    if(gp<0) gp+=ngp();
+
+    us vars_per_gp = 4;
+    vars_per_gp+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
+    
+    return sol.subvec((gp*vars_per_gp+varnr)*Ns,
+                      (gp*vars_per_gp+varnr+1)*Ns-1);
+}
+vd Duct::getvarx(const TaSystem& sys,int varnr,int t) const {
+    vd res(ngp());
+    for(us i=0;i<ngp();i++){
+        res(i) = getvart(sys,varnr,i)(t);
+    }
+    return res;
 }
 vd Duct::initialSolution(const TaSystem& sys) const {
     
@@ -87,7 +203,7 @@ vd Duct::initialSolution(const TaSystem& sys) const {
         // Initial density
 
         initsol.subvec((i*vars_per_gp+0)*Ns,(i*vars_per_gp+1)*Ns-1) =
-            gas.rho0();
+            gas.rho0()+.01;
         
         // Initial velocity
         initsol.subvec((i*vars_per_gp+1)*Ns,(i*vars_per_gp+2)*Ns-1) =
@@ -115,20 +231,35 @@ vd Duct::initialSolution(const TaSystem& sys) const {
 us Duct::getNEqs(const TaSystem& sys) const {
     TRACE(15,"Duct::getNEqs()");    
     us Ns = sys.Ns();
+
+    // The number of equations per gridpoint. We have: continuity,
+    // momentum, energy, and state
     us number_eqs = 4;
 
     // When we have to solve a solid heat balance
-    number_eqs+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
+    number_eqs+= (_ductpb.stempmodel() == pb::HeatBalance ?  : 0);
 
-    return Ns*number_eqs*(ngp()-1);
+    us neqs = Ns*number_eqs*(ngp()-1);
 
+    // For the last gridpoint, we also have an equation of state
+    neqs += Ns;
+
+    // We also have an extra equation for isentropic. For the energy
+    // transport equation, this would result in a boundary condition
+    if(_ductpb.htmodel() == pb::Isentropic) {
+        neqs += Ns;
+    }
+
+    VARTRACE(15,neqs);
+    return neqs;
 }
 us Duct::getNDofs(const TaSystem& sys) const {
     TRACE(15,"Duct::getNDofs()");        
     us Ns = sys.Ns();  
 
-    // rho,u,T,p,Ts
+    // rho,u,T,p
     us nvars_per_gp = 4;
+    // Ts maybe
     nvars_per_gp += (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
 
     return Ns*nvars_per_gp*ngp();

src/duct/duct.h

@@ -11,6 +11,7 @@
 #include "segment.h"
 #include "duct.pb.h"
 #include "geom.h"
+#include "tasmet_constants.h" // For the variable nrs
 
 class Equation;
 class Drag;
@@ -43,6 +44,28 @@ public:
     
     const pb::Duct& getDuctPb() const { return _ductpb;}
 
+    // Obtain values as a function of position, for a given time
+    // instance.
+    vd rhox(const TaSystem& sys,int t) const { return getvarx(sys,constants::rho,t);}
+    vd ux(const TaSystem& sys,int t) const { return getvarx(sys,constants::u,t); }
+    vd Tx(const TaSystem& sys,int t) const { return getvarx(sys,constants::T,t); }
+    vd px(const TaSystem& sys,int t) const { return getvarx(sys,constants::p,t); }
+    vd Tsx(const TaSystem& sys,int t) const { return getvarx(sys,constants::Ts,t); }
+
+    vd rhot(const TaSystem& sys,int gp) const { return getvart(sys,constants::rho,gp); }
+    vd ut(const TaSystem& sys,int gp) const { return getvart(sys,constants::u,gp); }
+    vd Tt(const TaSystem& sys,int gp) const { return getvart(sys,constants::T,gp); }
+    vd pt(const TaSystem& sys,int gp) const { return getvart(sys,constants::p,gp); }
+    vd Tst(const TaSystem& sys,int gp) const { return getvart(sys,constants::Ts,gp); }
+
+    /// Obtain variable as a function of time for a given grid point
+    vd getvart(const TaSystem& sys,int varnr,int gp) const;
+
+    /// Obtain variable as a function of position, for a given time
+    /// instance
+    vd getvarx(const TaSystem& sys,int varnr,int t) const;
+
+    d getvartx(const TaSystem& sys,int t,int gp) const;
     // Solving
     virtual void residual(const TaSystem&,arma::subview_col<d>&& residual) const;
 

src/ductbc/adiabaticwall.cpp

@@ -44,16 +44,44 @@ void AdiabaticWall::residual(const TaSystem& sys,
                              ) const {
 
     TRACE(15,"AdiabaticWall::residual()");
+    const pb::Duct& dpb = getDuct(sys).getDuctPb();
+    us Ns = sys.Ns();
+
+    const Duct& duct = getDuct(sys);
+    if(_side == pb::left) {
+        residual.subvec(0,Ns-1) = duct.ut(sys,0);
+        if(dpb.htmodel() != pb::Isentropic) {
+            // TODO: Put the boundary condition of zero heat flux here
+            // residual.subvec(Ns,2*Ns-1) = 
+            tasmet_assert(false,"");
+        }
+    }
+    else {
+        residual.subvec(0,Ns-1) = duct.ut(sys,-1);
+        if(dpb.htmodel() != pb::Isentropic) {
+            // TODO: Put the boundary condition of zero heat flux here
+            // residual.subvec(Ns,2*Ns-1) = duct.Tt(sys,-1) - _T;
+            tasmet_assert(false,"");
+        }
+    }
 
 }
 
 us AdiabaticWall::getNEqs(const TaSystem& sys) const {
     TRACE(15,"AdiabaticWall::getNEqs()");    
     // u = 0
-    // dT/dx = 0
+    // dT/dx = 0 --> if htmodel is not Isentropic
-    // dTs/dx = 0 => 3 eqs
+    // dTs/dx = 0 => if stempmodel is not Prescribed
-    bool has_solideq = getDuct(sys).getDuctPb().stempmodel() != pb::Prescribed;
+
-    return sys.Ns()*(has_solideq ? 3: 2);
+    const pb::Duct& dpb = getDuct(sys).getDuctPb();
+
+    bool has_solideq = dpb.stempmodel() != pb::Prescribed;
+
+    us neqs = sys.Ns()*(has_solideq ? 2: 1);
+    if(dpb.htmodel() != pb::Isentropic) neqs+= sys.Ns();
+
+    VARTRACE(15,neqs);
+    return neqs;
 }
 void AdiabaticWall::show(const TaSystem&,us verbosity_level) const {
     TRACE(15,"AdiabaticWall::show()");    

src/ductbc/pressurebc.cpp

@@ -11,6 +11,9 @@
 #include "tasmet_tracer.h"
 #include "tasystem.h"
 #include "duct.h"
+#include "tasmet_evalscript.h"
+#include "perfectgas.h"
+#include "adiabatictemp.h"
 
 PressureBc::PressureBc(const us id,
                        const TaSystem& sys,
@@ -18,13 +21,45 @@ PressureBc::PressureBc(const us id,
     DuctBc(id,dbc)
 {
     TRACE(15,"PressureBc(id,sys,dbc)");
+    vd time = sys.timeInstances();
+    VARTRACE(15,time);
+    
+    EvaluateFun pfun(dbc.pressure(),
+                     "Error in evaluating prescribed pressure",
+                     "t");
+                     
+    pfun.addGlobalDef("f",sys.getfreq());
+    pfun.addGlobalDef("omg",sys.getomg());
+    pfun.addGlobalDef("p0",sys.gas().p0());
+    pfun.addGlobalDef("T0",sys.gas().T0());
+
+    _p = pfun(time);
+
+    if(dbc.isentropic()) {
+        auto& gas = dynamic_cast<const PerfectGas&>(sys.gas());
+        _T = adiabaticTemp(gas,_p);
+    }
+
+    else {
+        EvaluateFun Tfun(dbc.temperature(),
+                         "Error in evaluating prescribed temperature",
+                         "t");
+                     
+        Tfun.addGlobalDef("f",sys.getfreq());
+        Tfun.addGlobalDef("omg",sys.getomg());
+        Tfun.addGlobalDef("p0",sys.gas().p0());
+        Tfun.addGlobalDef("T0",sys.gas().T0());
+        _T = Tfun(time);
+
+    }
 
 }
 PressureBc::PressureBc(const PressureBc& o):
-    DuctBc(o) {
+    DuctBc(o),_p(o._p),_T(o._T) {
 
     TRACE(15,"PressureBc(o)");
 
+
 }
 PressureBc::~PressureBc() {
     
@@ -36,23 +71,49 @@ vd PressureBc::initialSolution(const TaSystem& sys) const {
     return vd();
 }
 
-void PressureBc::residual(const TaSystem&,
+void PressureBc::residual(const TaSystem& sys,
                           arma::subview_col<d>&& residual
                           ) const {
 
     TRACE(15,"PressureBc::residual()");
 
+    const pb::Duct& dpb = getDuct(sys).getDuctPb();
+    us Ns = sys.Ns();
 
+    const Duct& duct = getDuct(sys);
+    if(_side == pb::left) {
+        residual.subvec(0,Ns-1) = duct.pt(sys,0) - _p;
+
+        if(dpb.htmodel() != pb::Isentropic) {
+            residual.subvec(Ns,2*Ns-1) = duct.Tt(sys,0) - _T;
+        }
+    }
+    else {
+        residual.subvec(0,Ns-1) = duct.pt(sys,-1) - _p;
+
+        if(dpb.htmodel() != pb::Isentropic) {
+            residual.subvec(Ns,2*Ns-1) = duct.Tt(sys,-1) - _T;
+        }
+    }
 }
 
 us PressureBc::getNEqs(const TaSystem& sys) const {
     TRACE(15,"PressureBc::getNEqs()");    
-    // p = x
+
-    // T = x
+    // We provide equations for
-    // This one only if the duct solves for solid
+    // p = prescribed
-    // Ts = x => 3 equations
+    // T = prescribed, if htmodel of duct is not Isentropic
-    bool has_solideq = getDuct(sys).getDuctPb().stempmodel() != pb::Prescribed;
+    // Ts = prescribed, if stempmodel of duct is not Prescribed
-    return sys.Ns()*(has_solideq ? 3: 2);
+
+    const pb::Duct& dpb = getDuct(sys).getDuctPb();
+
+    bool has_solideq = dpb.stempmodel() != pb::Prescribed;
+
+    us neqs = sys.Ns()*(has_solideq ? 2: 1);
+    if(dpb.htmodel() != pb::Isentropic) neqs+= sys.Ns();
+
+    VARTRACE(15,neqs);
+    return neqs;
 }
 void PressureBc::show(const TaSystem&,us verbosity_level) const {
     TRACE(15,"PressureBc::show()");    

src/ductbc/pressurebc.h

@@ -14,9 +14,17 @@
 class TaSystem;
 class Variable;
 
-class PressureBc: public DuctBc {
+/**
-    Variable *_p,*_T,*_Ts;
+ * A PressureBc is a Duct boundary condition, where the pressure is
+ * prescribed on one side of a Duct. Besides the pressure it is
+ * necessary to also prescribe the temperature, and if a model is used
+ * to thermally interact with the solid, also a solid temperature is
+ * prescribed.
+ */
 
+class PressureBc: public DuctBc {
+    vd _p,_T,_Ts;               /**< Prescribed values for pressure,
+                                   temperature and solid temperature */
     pb::DuctSide _side;         /**< Duct side at which this b.c. works */
 protected:
     PressureBc(const PressureBc&);

src/material/adiabatictemp.cpp

@@ -88,15 +88,13 @@ namespace {
         }
     }
 }
-Variable adiabaticTemp(const PerfectGas& gas,
+vd adiabaticTemp(const PerfectGas& gas,
-                       const Variable& pressure) {
+                 const vd& p) {
                        
     
     TRACE(10,"adiabaticTemp()");
-    Variable temp(pressure);
-    const vd p = pressure.tdata();
 
-    vd T = temp.tdata();
+    vd T(p.size());
 
     AdiabaticTemp adt(gas,p(0)); // Equation
 
@@ -115,8 +113,7 @@ Variable adiabaticTemp(const PerfectGas& gas,
         VARTRACE(15,T(i));
 
     }
-    temp.settdata(T);
+    return T;
-    return temp;
 }
 //////////////////////////////////////////////////////////////////////
 

src/material/adiabatictemp.h

@@ -3,17 +3,24 @@
 // 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
-
+#include "tasmet_types.h"
 class PerfectGas;
-class Variable;
+/** 
-
+ * Compute the adiabatic/isentropic temperature as a function of the pressure
-Variable adiabaticTemp(const PerfectGas& gas,const Variable& pressure);
+ * for a thermally perfect gas.
+ *
+ * @param gas: The gas to compute the temperature for. Be aware of the
+ * value for T0 and p0 in the Gas.
+ * @param pressure the pressure to compute the temperature for
+ *
+ * @return the temperature
+ */
+vd adiabaticTemp(const PerfectGas& gas,const vd& pressure);
 
 #endif // ADIABATICTEMP_H
 //////////////////////////////////////////////////////////////////////

src/material/gas.h

@@ -78,6 +78,16 @@ public:
         element_wise(rho);
     }
 
+    virtual d drhodT(d T,d p) const=0;
+    vd drhodT(const vd& T,const vd& p) const {
+        element_wise(drhodT);
+    }
+
+    virtual d drhodp(d T,d p) const=0;
+    vd drhodp(const vd& T,const vd& p) const {
+        element_wise(drhodp);
+    }
+
     // Adiabatic speed of sound
     virtual d cm(d T,d p) const=0;
     vd cm(const vd& T,const vd& p) const {

src/material/perfectgas.h

@@ -63,6 +63,13 @@ public:
         checkzero(T);
         return p/Rs()/T;
     }
+    d drhodT(d T,d p) const {
+        return -p/Rs()/pow(T,2);
+    }
+    d drhodp(d T,d p) const {
+        return Rs()/T;
+    }
+
     d cv(d T,d p) const {
         return cp(T,p)-Rs();
     }

src/solver/newton_raphson.cpp

@@ -9,6 +9,8 @@
 #include "newton_raphson.h"
 #include "tasmet_tracer.h"
 
+#define DEBUG_TASMET_SYSTEM
+
 void NewtonRaphson::start_implementation(GradientNonlinearSystem& system,
                                          progress_callback* callback) {
 
@@ -23,6 +25,17 @@ void NewtonRaphson::start_implementation(GradientNonlinearSystem& system,
 
     vd residual=system.residual();
 
+    #ifdef DEBUG_TASMET_SYSTEM
+    cout << "Initial solution: " << endl;
+    cout << guess << endl;
+    #endif  // DEBUG_TASMET_SYSTEM
+    
+    #ifdef DEBUG_TASMET_SYSTEM
+    cout << "Initial residual: " << endl;
+    cout << residual << endl;
+    #endif  // DEBUG_TASMET_SYSTEM
+
+
     SolverProgress progress;
     SolverAction action;
     
@@ -40,6 +53,11 @@ void NewtonRaphson::start_implementation(GradientNonlinearSystem& system,
         system.updateSolution(guess);
 
         residual = system.residual();
+        #ifdef DEBUG_TASMET_SYSTEM
+        cout << "Residual: ";
+        cout << residual << endl;
+        #endif  // DEBUG_TASMET_SYSTEM
+
     
         progress.rel_err = norm(dx);
         progress.fun_err = norm(residual);

src/sys/globalconf.cpp

@@ -48,6 +48,11 @@ GlobalConf::GlobalConf(us Nf,d freq):
     setfreq(freq);
 
 }
+vd GlobalConf::timeInstances() const {
+    us Ns = this->Ns();
+    d lastt = (Ns-1)/(Ns*getfreq());
+    return arma::linspace<vd>(0,lastt,Ns);
+}
 void GlobalConf::show() const {
     cout << "------- Global configuration ------ \n";        
     cout << "------- Number of harmonics to solve for: "<< _Nf <<"\n";

src/sys/globalconf.h

@@ -44,6 +44,8 @@ public:
     d getfreq() const {return _omg/2/number_pi;}
     // d meshPeclet(const Gas& gas,d dx,d u) const {return u*dx*gas.rho0()*gas().cp(T0())/gas().kappa(T0());}
 
+    vd timeInstances() const;
+
     void setfreq(d freq){setomg(2*number_pi*freq);}
     void setomg(d omg);
 

src/sys/tasystem.cpp

@@ -43,7 +43,7 @@ TaSystem::TaSystem(const pb::System& sys):
             throw e;
         }
     }
-    // Create all ducts
+    // Create all ductbcs
     for(const auto& d : sys.ductbcs()) {
         // d.first: id
         // d.second: duct description
@@ -61,6 +61,33 @@ TaSystem::TaSystem(const pb::System& sys):
         }
     }
     
+    // Create the initial solution from the segments and store it
+    // here. Please be careful not to call any virtual functions!
+    vus ndofs = getNDofs();
+    vus dofend = arma::cumsum(ndofs);
+    us total_dofs = arma::sum(ndofs);
+    vd solution = vd(total_dofs);
+
+    us i=0;
+    const Segment* seg;
+    for(auto& seg_: _segs) {
+        seg = seg_.second;
+
+        if(ndofs(i)>0) {
+            if(i==0) {
+                solution.subvec(0,ndofs(0)-1) =
+                    seg->initialSolution(*this);
+            }
+            else {
+                solution.subvec(dofend(i-1),dofend(i)-1) =
+                    seg->initialSolution(*this);
+            }
+            i++;
+        }
+    }
+    // TODO: work directly on the final solution array
+    _solution = solution;
+
     // Copy solution vector, if valid
     // const auto& sol = sys.solution();
     // us size = sol.size(), i=0;
@@ -75,7 +102,8 @@ TaSystem::TaSystem(const pb::System& sys):
 }
 TaSystem::TaSystem(const TaSystem& o):
     GlobalConf(o),                 // Share a ptr to the Global conf
-    _gas(o._gas->copy())
+    _gas(o._gas->copy()),
+    _solution(o._solution)
 {
     TRACE(25,"TaSystem::TaSystem(TaSystem&) copy");
 
@@ -150,6 +178,10 @@ vd TaSystem::residual() const {
 
     vd residual(total_neqs);
 
+    #ifdef TASMET_DEBUG
+    residual = arma::datum::nan*ones(total_neqs);
+    #endif
+
     us i=0;
     const Segment* seg;
     for(auto seg_: _segs) {
@@ -172,6 +204,8 @@ vd TaSystem::residual() const {
         i++;
     }
 
+    TRACE(15,"Obtained residual from all Segments");
+
     #ifdef TASMET_DEBUG
     assert(arbitrateMassEq< (int) total_neqs);    
     #endif // TASMET_DEBUG
@@ -185,43 +219,16 @@ vd TaSystem::residual() const {
 }
 vd TaSystem::getSolution() const {
 
-    if(_solution.size() == 0) {
-        // Create the initial solution from the segments
-        // and store it here
-
-        vus ndofs = getNDofs();
-        vus dofend = arma::cumsum(ndofs);
-        us total_dofs = arma::sum(ndofs);
-        vd solution = vd(total_dofs);
-
-        us i=0;
-        const Segment* seg;
-        for(auto& seg_: _segs) {
-            seg = seg_.second;
-
-            if(ndofs(i)>0) {
-                if(i==0) {
-                    solution.subvec(0,ndofs(0)-1) =
-                        seg->initialSolution(*this);
-                }
-                else {
-                    solution.subvec(dofend(i-1),dofend(i)-1) =
-                        seg->initialSolution(*this);
-                }
-                i++;
-            }
-        }
-        return solution;
-    } // if the solution did not yet exist
-
     return _solution;
 }
 const arma::subview_col<d> TaSystem::getSolution(const us seg_id) const {
 
     vus ndofs = getNDofs();
     vus dofsend = arma::cumsum(ndofs);
-
+    VARTRACE(15,dofsend);
+    VARTRACE(15,seg_id);
     if(seg_id == 0) {
+        VARTRACE(15,_solution.size());
         return _solution.subvec(0,dofsend(0)-1);
     }
     else {

src/tasmet_constants.h

@@ -96,14 +96,11 @@ namespace constants {
 
     // These variable numbers are important, as they determine the
     // position of these variables in the array in cell.h
-    // const int rho=1;
+    const int rho=0;
-    // const int m=2;
+    const int u=1;
-    // const int T=3;
+    const int T=2;
-    // const int p=4;
+    const int p=3;
-    // const int Ts=5;
+    const int Ts=4;
-    // Number of variables
-    const int nvars_reserve=7;
-    const int neqs_reserve=7;
   
     const char* const model_fileext = ".tasmet";
 

src/tasmet_evalscript.cpp

@@ -47,7 +47,8 @@ inline void wrap_eval<void>(ChaiScript* chai,const string& script) {
 }
 
 EvaluateFun::EvaluateFun(const string& fun_return,
-                         const string& err_msg):
+                         const string& err_msg,
+                         const string& vars):
     _err_msg(err_msg),
     _fun_return(fun_return)
 {
@@ -55,7 +56,7 @@ EvaluateFun::EvaluateFun(const string& fun_return,
     _chai = getChaiScriptInstance();
     if(!_chai) throw TaSMETBadAlloc();
 
-    string script = "def myfun(x) {\n";
+    string script = string("def myfun(") + vars + ") {\n";
     script += "return " + fun_return + "; }\n";
     
     wrap_eval<void>(_chai.get(),script);

src/tasmet_evalscript.h

@@ -23,7 +23,7 @@ class EvaluateFun {
     string _err_msg;
     string _fun_return;
 public:
-    EvaluateFun(const string& fun_return,const string& err_msg = "Script error");
+    EvaluateFun(const string& fun_return,const string& err_msg = "Script error",const string& vars = "x");
     
     // Add a global definition to the namespace
     void addGlobalDef(const string& name,

testing/solver/brent_test.cpp

@@ -66,8 +66,8 @@ int main(){
     Helium helium(293.15,101325);
 
     gc_ptr gc(new GlobalConf(10,100));
-    Variable pressure(gc,10*101325);
+    vd pressure(1); pressure(0) = 10*101325;
-    Variable temperature = adiabaticTemp(helium,pressure);
+    vd temperature = adiabaticTemp(helium,pressure);
 }