Working! New way of obtaining Jacobian succesfully implemented. Newton-Raphson solver working

src/CMakeLists.txt

@@ -65,9 +65,8 @@ add_library(tasmet_src
   # 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
 )

src/duct/duct.cpp

@@ -60,12 +60,11 @@ Duct::~Duct() {
     //     delete eq;
     // }
 }
-void Duct::residualJac(arma::subview_col<d> && residual) const {
+void Duct::residualJac(SegPositionMapper& residual,
+                       SegJacRows& jacrows) const {
 
-    TRACE(15,"Duct::residual()");
+    TRACE(15,"Duct::residualJac()");
     
-    const arma::subview_col<d> sol = sys.getSolution(_id);
-
     vd rho,u,T,p,Ts;            // Solution at this gp
     vd rhop,up,Tp,pp,Tsp;       // Solution at next gp
 
@@ -78,28 +77,26 @@ void Duct::residualJac(arma::subview_col<d> && residual) const {
     us number_eqs = 4;
     number_eqs += (has_solideq) ? 1 : 0;
 
+    // rho,u,T,p
+    // Ts maybe
+    us nvars_per_gp = (_ductpb.stempmodel()
+                       == pb::HeatBalance ? 5 : 4);
+
     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(constants::rho,0);
     up = getvart(constants::u,0);
     Tp = getvart(constants::T,0);
-    pp   = getvart(constants::p,0);
+    pp = getvart(constants::p,0);
 
     const Gas& gas = sys.gas();
 
+    us i=0;
+
     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
@@ -112,8 +109,35 @@ void Duct::residualJac(arma::subview_col<d> && residual) const {
         pp   = getvart(constants::p,gp+1);
 
         cont  = ((rhop%up)-(rho%u))/dx;
+        *residual.at(i) = cont;
+
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::u}] =
+            -diagmat(rho)/dx;
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] =
+            -diagmat(u)/dx;
+        jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::rho}] =
+            diagmat(up)/dx;
+        jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::u}] =
+            diagmat(rhop)/dx;
+        TRACE(15,"SFSG");
+        i++;
 
         mom = (rhop%up%up - rho%u%u + pp - p)/dx;
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::u}] =
+            -diagmat(rho%u)/dx;
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] =
+            -diagmat(u%u)/dx;
+        jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::rho}] =
+            diagmat(u%u)/dx;
+        jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::u}] =
+            diagmat(up%rhop)/dx;
+        jacrows.at(i)[{_id,(gp)*nvars_per_gp+constants::p}] =
+            -eye(Ns,Ns)/dx;
+        jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::p}] =
+            eye(Ns,Ns)/dx;
+        *residual.at(i) = mom;
+        TRACE(15,"SFSG");
+        i++;
         
         switch (_ductpb.htmodel()) {
         case pb::Isentropic: {
@@ -125,33 +149,53 @@ void Duct::residualJac(arma::subview_col<d> && residual) const {
 
            en = p/p0 - pow(rho/rho0,gamma0);
 
+           *residual.at(i) = en;
+
+           jacrows.at(i)[{_id,gp*nvars_per_gp+constants::p}] =
+               eye(Ns,Ns)/p0;
+           jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] =
+               -(gamma0/rho0)*diagmat(pow(rho/rho0,gamma0-1));
+           TRACE(15,"SFSG");
         }
             break;
         default:
             tasmet_assert(false,"Not implemented htmodel");
         }
 
+        i++;
+        TRACE(15,"SFSG");
+        // Equation of state
         st = gas.rho(T,p) - rho;
-
+        TRACE(15,"SFSG");
-        residual.subvec(eq_offset+0*Ns,eq_offset+1*Ns-1) = cont;
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] =
-        residual.subvec(eq_offset+1*Ns,eq_offset+2*Ns-1) = mom;
+            -eye(Ns,Ns);
-        residual.subvec(eq_offset+2*Ns,eq_offset+3*Ns-1) = en;
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::p}] =
-        residual.subvec(eq_offset+3*Ns,eq_offset+4*Ns-1) = st;
+            diagmat(gas.drhodp(T,p));
-
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::T}] =
+            diagmat(gas.drhodT(T,p));
+           
+        *residual.at(i) = st; 
+        TRACE(15,"SFSG");
+        i++;
     }
 
-    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;
+    *residual.at(i) = st; 
+
+    jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::rho}] =
+        -eye(Ns,Ns);
+    jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::p}] =
+        diagmat(gas.drhodp(T,p));
+    jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::T}] =
+        diagmat(gas.drhodT(T,p));
+
+    i++;
 
     // 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();
@@ -159,24 +203,43 @@ void Duct::residualJac(arma::subview_col<d> && residual) const {
 
         en = p/p0 - pow(rho/rho0,gamma0);
 
-        residual.subvec(eq_offset,eq_offset+Ns-1) = en;
+        jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::p}] =
+            eye(Ns,Ns)/p0;
+        jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::rho}] =
+            -(gamma0/rho0)*diagmat(pow(rho/rho0,gamma0-1));
+
+        *residual.at(i) = en; i++;
     }
 
 }
+PosId Duct::getDof(int varnr,int gp) const {
+    // Wraparound
+    if(gp<0) gp+=ngp();
+
+    // rho,u,T,p
+    // Ts maybe
+    us nvars_per_gp = (_ductpb.stempmodel()
+                       == pb::HeatBalance ? 5 : 4);
+
+
+    return {_id,nvars_per_gp*gp+varnr};
+}
+
 vd Duct::getvart(int varnr,int gp) const {
-    TRACE(15,"Duct::getvart()");
+    TRACE(14,"Duct::getvart()");
-    const arma::subview_col<d> sol = sys.getSolution(_id);
+    const SegPositionMapper& sol = sys.getSolution(_id);
 
     us Ns = sys.Ns();
 
     // Wraparound
     if(gp<0) gp+=ngp();
 
-    us vars_per_gp = 4;
+    // rho,u,T,p
-    vars_per_gp+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
+    // Ts maybe
+    us nvars_per_gp = (_ductpb.stempmodel()
+                       == pb::HeatBalance ? 5 : 4);
     
-    return sol.subvec((gp*vars_per_gp+varnr)*Ns,
+    return *sol.at((gp*nvars_per_gp+varnr));
-                      (gp*vars_per_gp+varnr+1)*Ns-1);
 }
 vd Duct::getvarx(int varnr,int t) const {
     vd res(ngp());
@@ -185,48 +248,44 @@ vd Duct::getvarx(int varnr,int t) const {
     }
     return res;
 }
-vd Duct::initialSolution() const {
+void Duct::initialSolution(SegPositionMapper& sol) const {
     
     TRACE(15,"Duct::initialSolution()");
 
-    vd initsol(getNDofs());
+    VARTRACE(15,sol.size());
-    VARTRACE(15,initsol.size());
 
-    us vars_per_gp = 4;
+    us nvars_per_gp = 4;
-    vars_per_gp+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
+    nvars_per_gp+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
-
-    us Ns = sys.Ns();
 
     const Gas& gas = sys.gas();
-    
+
+    us segdof = 0;
+
     for(us i=0;i<ngp();i++){
         VARTRACE(15,i);
         // Initial density
-
+        *sol.at(segdof) += gas.rho0()+0.1;
-        initsol.subvec((i*vars_per_gp+0)*Ns,(i*vars_per_gp+1)*Ns-1) =
+        segdof++;
-            gas.rho0()+.01;
         
         // Initial velocity
-        initsol.subvec((i*vars_per_gp+1)*Ns,(i*vars_per_gp+2)*Ns-1) =
+        *sol.at(segdof) += 0.1;
-            0;
+        segdof++;
 
         // Initial Temperature
-        initsol.subvec((i*vars_per_gp+2)*Ns,(i*vars_per_gp+3)*Ns-1) =
+        *sol.at(segdof) += _Tsprescribed(i);
-            _Tsprescribed(i);
+        segdof++;
 
         // Initial pressure
-        initsol.subvec((i*vars_per_gp+3)*Ns,(i*vars_per_gp+4)*Ns-1) =
+        *sol.at(segdof) += gas.p0();
-            gas.p0();
+        segdof++;
-        
+
         // Initial solid temperature, if not prescribed
         if(_ductpb.stempmodel() != pb::Prescribed) {
-            initsol.subvec((i*vars_per_gp+4)*Ns,(i*vars_per_gp+5)*Ns-1) =
+            *sol.at(segdof) += _Tsprescribed(i);
-                _Tsprescribed(i);
+            segdof++;
         }
     }
 
-    return initsol;
-
 }
 
 us Duct::getNEqs() const {
@@ -257,9 +316,9 @@ us Duct::getNDofs() const {
     TRACE(15,"Duct::getNDofs()");        
  
     // rho,u,T,p
-    us nvars_per_gp = 4;
     // Ts maybe
-    nvars_per_gp += (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
+    us nvars_per_gp = (_ductpb.stempmodel()
+                       == pb::HeatBalance ? 5 : 4);
 
     return nvars_per_gp*ngp();
 }

src/duct/duct.h

@@ -12,6 +12,7 @@
 #include "duct.pb.h"
 #include "geom.h"
 #include "tasmet_constants.h" // For the variable nrs
+#include "jacobian.h"
 
 class Equation;
 class Drag;
@@ -34,8 +35,9 @@ public:
 
     vd Tsprescribed() const { return _Tsprescribed;}
 
-    // Initialize the solution to something sensible
+    void initialSolution(SegPositionMapper&) const;
-    vd initializeSolution(const TaSystem& sys);
+
+
 
     virtual Duct* copy(const TaSystem&) const;
     const Geom& geom() const;
@@ -70,10 +72,10 @@ public:
     vd getvarx(int varnr,int t) const;
 
     d getvartx(int t,int gp) const;
-    // Solving
-    virtual void residualJac(arma::subview_col<d>&& residual) const;
 
-    vd initialSolution() const;
+    // Solving
+    virtual void residualJac(SegPositionMapper& residual,
+                             SegJacRows& jac) const;
 
     // virtual void getSolution(const TaSystem&,const us insertion_start) const;
 

src/ductbc/adiabaticwall.cpp

@@ -35,20 +35,22 @@ AdiabaticWall::~AdiabaticWall() {
 AdiabaticWall* AdiabaticWall::copy(const TaSystem& sys) const {
     return new AdiabaticWall(sys,*this);
 }
-vd AdiabaticWall::initialSolution() const {
-    return vd();
-}
 
-void AdiabaticWall::residualJac(arma::subview_col<d>&& residual
-                             ) const {
 
-    TRACE(15,"AdiabaticWall::residual()");
+void AdiabaticWall::residualJac(SegPositionMapper& residual,
+                                SegJacRows& jac) const {
+
+    TRACE(15,"AdiabaticWall::residualJac()");
     const Duct& duct = getDuct();
     const pb::Duct& dpb = duct.getDuctPb();
     us Ns = sys.Ns();
 
+    us i=0;
+
     if(_side == pb::left) {
-        residual.subvec(0,Ns-1) = duct.ut(0);
+        *residual.at(i) = duct.ut(0);
+        jac.at(i)[duct.getDof(constants::u,0)].eye(Ns,Ns);
+
         if(dpb.htmodel() != pb::Isentropic) {
             // TODO: Put the boundary condition of zero heat flux here
             // residual.subvec(Ns,2*Ns-1) = 
@@ -56,7 +58,8 @@ void AdiabaticWall::residualJac(arma::subview_col<d>&& residual
         }
     }
     else {
-        residual.subvec(0,Ns-1) = duct.ut(-1);
+        *residual.at(i) = duct.ut(-1);
+        jac.at(i)[duct.getDof(constants::u,-1)].eye(Ns,Ns);
         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;

src/ductbc/adiabaticwall.h

@@ -34,10 +34,8 @@ public:
     ~AdiabaticWall();
     AdiabaticWall* copy(const TaSystem& sys) const;
 
-    vd initialSolution() const;
+    virtual void residualJac(SegPositionMapper& residual,
-
+                             SegJacRows& jac) const;
-    virtual void residualJac(arma::subview_col<d>&& residual
-                          ) const;
 
     // Return the total number of equations in this segment
     virtual us getNEqs() const;
@@ -50,11 +48,6 @@ public:
     // Reset amplitude data in higher harmonics
     // virtual void resetHarmonics() = 0;
 
-    // Fill Jacobian with values from the equations in this
-    // segment/connector.
-    virtual void jac(const TaSystem&,Jacobian&,us dof_start,us eq_start) const;
-                  
-
 };
 
 #endif // ADIABATICWALL_H

src/ductbc/ductbc.cpp

@@ -31,10 +31,10 @@ DuctBc* DuctBc::newDuctBc(const us id,
         return new PressureBc(sys,id,dbc);
         break;
     }
-    // case pb::AdiabaticWall: {
+    case pb::AdiabaticWall: {
-    //     return new AdiabaticWall(id,sys,dbc);
+        return new AdiabaticWall(sys,id,dbc);
-    //     break;
+        break;
-    // }
+    }
     default:
         tasmet_assert(false,"Not implemented DuctBc");
         break;

src/ductbc/ductbc.h

@@ -33,6 +33,8 @@ public:
 
     const Duct& getDuct() const;
 
+    virtual void initialSolution(SegPositionMapper&) const {}
+
 };
 
 

src/ductbc/pressurebc.cpp

@@ -69,31 +69,28 @@ PressureBc::~PressureBc() {
 PressureBc* PressureBc::copy(const TaSystem& sys) const {
     return new PressureBc(sys,*this);
 }
-vd PressureBc::initialSolution() const {
-    return vd();
-}
 
-void PressureBc::residualJac(arma::subview_col<d>&& residual
+void PressureBc::residualJac(SegPositionMapper& residual,
-                          ) const {
+                             SegJacRows& jacrows) const {
 
-    TRACE(15,"PressureBc::residual()");
+    TRACE(15,"PressureBc::residualJac()");
 
     const pb::Duct& dpb = getDuct().getDuctPb();
     us Ns = sys.Ns();
 
     const Duct& duct = getDuct();
     if(_side == pb::left) {
-        residual.subvec(0,Ns-1) = duct.pt(0) - _p;
+        *residual.at(0) = duct.pt(0) - _p;
-
+        jacrows.at(0)[duct.getDof(constants::p,0)] = eye(Ns,Ns);
         if(dpb.htmodel() != pb::Isentropic) {
-            residual.subvec(Ns,2*Ns-1) = duct.Tt(0) - _T;
+            // residual.subvec(Ns,2*Ns-1) = duct.Tt(0) - _T;
         }
     }
     else {
-        residual.subvec(0,Ns-1) = duct.pt(-1) - _p;
+        *residual.at(0) = duct.pt(-1) - _p;
-
+        jacrows.at(0)[duct.getDof(constants::p,-1)] = eye(Ns,Ns);
         if(dpb.htmodel() != pb::Isentropic) {
-            residual.subvec(Ns,2*Ns-1) = duct.Tt(-1) - _T;
+            // residual.subvec(Ns,2*Ns-1) = duct.Tt(-1) - _T;
         }
     }
 }

src/ductbc/pressurebc.h

@@ -38,10 +38,8 @@ public:
     
     PressureBc* copy(const TaSystem&) const;
 
-    vd initialSolution() const;
+    virtual void residualJac(SegPositionMapper& residual,
-
+                             SegJacRows& jac) const;
-    virtual void residualJac(arma::subview_col<d>&& residual
-                          ) const;
 
     // Return the total number of equations in this segment
     virtual us getNEqs() const;

src/material/gas.h

@@ -21,11 +21,11 @@ using pb::air;
 using pb::helium;
 using pb::nitrogen;
 
-#define element_wise(varname) \
+#define element_wise(function) \
-    vd varname_(T.size());     \
+    vd function_(T.size());     \
     for(us i=0;i<T.size();i++) \
-        varname_ = varname(T(i),p(i));\
+        function_(i) = function(T(i),p(i));     \
-    return varname_
+    return function_
 
 
 

src/solver/newton_raphson.cpp

@@ -34,27 +34,23 @@ void NewtonRaphson::start_implementation(GradientNonlinearSystem& system,
     ResidualJac resjac;
     system.residualJac(resjac);
 
-    const sdmat& jac = resjac.jacobian;
-    const vd& residual = resjac.residual;
     vd dx;
 
-    #ifdef DEBUG_TASMET_SYSTEM
+    assert(resjac.jacobian.n_cols==resjac.residual.size());
-    cout << "Initial residual: " << endl;
+    assert(resjac.jacobian.n_rows==resjac.residual.size());
-    cout << residual << endl;
-    #endif  // DEBUG_TASMET_SYSTEM
-
-    assert(jac.n_cols==residual.size());
-    assert(jac.n_rows==residual.size());
 
     while (true) {
 
         #ifdef DEBUG_TASMET_SYSTEM
         cout << "Residual: ";
-        cout << residual << endl;
+        cout << resjac.residual << endl;
+        cout << dmat(resjac.jacobian) << endl;        
         #endif  // DEBUG_TASMET_SYSTEM
     
-        dx = -1*_dampfac*arma::spsolve(jac,residual,"superlu");
+        TRACE(15,"Solving system of eqs");
-
+        dx = -1*_dampfac*arma::spsolve(resjac.jacobian,resjac.residual,"superlu");
+        TRACE(15,"Solving system of eqs done");
+        
         progress.rel_err = norm(dx);
 
         guess += dx;
@@ -64,12 +60,14 @@ void NewtonRaphson::start_implementation(GradientNonlinearSystem& system,
         // Obtain a new residual and Jacobian matrix
         system.residualJac(resjac);
 
-        progress.fun_err = norm(residual);
+        progress.fun_err = norm(resjac.residual);
         progress.iteration++;
 
         action = (*callback)(progress);
 
         if(action==Stop){
+            TRACE(15,"Solution at stop:");
+            TRACE(15,guess);
             break;
         }
             

src/sys/jaccol.cpp

@@ -1,13 +0,0 @@
-// jaccol.cpp
-//
-// last-edit-by: J.A. de Jong 
-// 
-// Description:
-//
-//////////////////////////////////////////////////////////////////////
-
-#include "jaccol.h"
-
-
-
-//////////////////////////////////////////////////////////////////////

src/sys/jaccol.h

@@ -1,50 +0,0 @@
-// jaccol.h
-//
-// Author: J.A. de Jong 
-//
-// Description:
-// A column block of the Jacobian matrix
-//////////////////////////////////////////////////////////////////////
-#pragma once
-#ifndef JACCOL_H
-#define JACCOL_H
-#include "tasmet_types.h"
-
-/**
- * A PosId is a specifier for the position in the global Jacobian
- * matrix. Both considering the column as well as the row
- * 
- */
-struct PosId {
-    int segno;
-    int segoffset;
-};
-
-class JacCol{                // Column block of Jacobian
-
-public:
-
-    // For more ease of use, we make the members public.
-
-    PosId id;                  /**< The Id of this DOF */
-    dmat data;                 /**< The matrix block */
-
-    JacCol(PosId id): id(id){}
-    JacCol(PosId id,const dmat& data):id(id),data(data){}
-
-    JacCol& operator=(const JacCol& o) { id = o.id; data=o.data; return *this; }
-    JacCol(const JacCol& o): id(o.id),data(o.data){}
-
-    void prePostMultiply(const dmat& pre,const dmat& post){ data=pre*data*post;}
-
-    // Some operations on the data
-    JacCol operator-() { return JacCol(id,-data);}
-
-    JacCol& operator+=(const dmat& add) {data+=add; return *this;}
-    JacCol& operator*=(const double d) {data*=d; return *this;}
-
-};
-  
-
-#endif // JACCOL_H
-//////////////////////////////////////////////////////////////////////

src/sys/jacobian.cpp

@@ -2,14 +2,6 @@
 #include "triplets.h"
 #include "tasmet_tracer.h"
 
-void Jacobian::operator+=(const Jacobian& other){
-    TRACE(2,"Jacobian::append()");
-    jacrows.insert(jacrows.end(),other.jacrows.begin(),other.jacrows.end());
-}
-void Jacobian::operator+=(const JacRow& other){
-    TRACE(2,"Jacobian::append()");
-    jacrows.push_back(other);
-}
 // Jacobian::operator TripletList() const{
 
 //     TRACE(18,"Jacobian::operator Tripletlist()");

src/sys/jacobian.h

@@ -8,25 +8,59 @@
 #pragma once
 #ifndef JACOBIAN_H
 #define JACOBIAN_H
+#include "tasmet_types.h"
 
-#include "jacrow.h"
+#include <unordered_map>
+#include <map>
 
 class TripletList;
+/**
+ * A PosId is a specifier for the position in the global Jacobian
+ * matrix. Both considering the column as well as the row
+ * 
+ */
+struct PosId {
+    us segno;                  /**< This is the segment number */
+    us segoffset;              /**< This is the offset in the segment */
 
-class Jacobian{
+    bool operator==(const PosId& other) const {
-    us _ndofs;
+        return (segno==other.segno) &&
-public:
+            (segoffset == other.segoffset);
-    Jacobian(us ndofs):_ndofs(ndofs){}
+    }
-    vector<JacRow> jacrows;
 
-    void operator+=(const Jacobian&);
+    bool operator<(const PosId& other) const {
-    void operator+=(const JacRow&);
+        if(other.segno < segno) return true;
+        else if(other.segno > segno) return false;
+        else return (other.segoffset < segoffset);
 
-    operator TripletList() const;
+    }
+};
+/**
+ * To use the PosId in an unordered_map, it requires a custom hash
+ * function. This struct has a custom function to generate a hash from
+ * the given two integers.
+ * 
+ */
+
+struct hash_PosId
+{
+    std::size_t operator() ( const PosId& m ) const
+    {
+        const us a = m.segno+1 ;
+        const us b = m.segoffset+1 ;
+        return std::hash<us>()( a*b + (a+b) ) ;
+    }
 };
 
-  
 
+typedef std::unordered_map<PosId,dmat,hash_PosId> JacRow;
+
+typedef std::vector<JacRow> SegJacRows;
+
+typedef std::map<us,SegJacRows> Jacobian;
+
+typedef std::vector<vd*> SegPositionMapper;
+typedef std::map<us, SegPositionMapper> GlobalPositionMapper;
 
 #endif // JACOBIAN_H
 //////////////////////////////////////////////////////////////////////

src/sys/jacrow.cpp

@@ -1,40 +0,0 @@
-#include "jacobian.h"
-#include "tasmet_variable.h"
-#include "globalconf.h"
-
-
-JacRow& JacRow::operator+=(const JacCol& j){
-    TRACE(10,"JacRow::operator+=(const JacCol& j)");
-    jaccols.push_back(j);
-    return *this;
-}
-
-JacRow& JacRow::operator*=(const d val){
-    TRACE(2,"Jacobian::operator*=()");
-    for(JacCol& col: jaccols)
-        col*=val;
-    return *this;
-}
-JacRow& JacRow::operator+=(const JacRow& other){
-    // Reserve some space
-    this->jaccols.reserve(this->jaccols.capacity()+other.jaccols.size()-this->jaccols.size());
-
-    for(const JacCol& col :other.jaccols)
-        (*this)+=col;
-
-    return *this;
-}
-JacRow JacRow::operator-() const{
-    TRACE(15,"JacRow::operator-()");
-
-    JacRow result(*this);
-    for (JacCol& jaccol : result.jaccols) 
-        jaccol*=-1;
-
-    return result;
-}
-void JacRow::prePostMultiply(const dmat& pre,const dmat& post) {
-    for(JacCol& jaccol: jaccols)
-        jaccol.prePostMultiply(pre,post);
-}
-

src/sys/jacrow.h

@@ -1,55 +0,0 @@
-// jacrow.h
-//
-// Author: J.A. de Jong 
-//
-// Description:
-//
-//////////////////////////////////////////////////////////////////////
-#pragma once
-#ifndef JACROW_H
-#define JACROW_H
-#include <vector>
-#include "jaccol.h"
-
-class JacRow{                 // Row in Jacobian matrix
-
-public:
-
-    PosId id;                   /**< Specifies the row for this
-                                   JacRow */
-    std::vector<JacCol> jaccols; /**< The columns in this row */
-
-    /** 
-     * Initialize a Jacobian row with a single column block
-     *
-     * @param id The equation ID
-     * @param col The added column
-     */
-    JacRow(PosId id,const JacCol& col):JacRow(id){jaccols.push_back(col);}
-
-    JacRow(PosId id): id(id){}
-
-    JacRow operator-() const;
-
-    JacRow& operator+=(const JacCol&);
-    JacRow& operator+=(const JacRow& jacrow);
-    JacRow& operator*=(const d val); // Multiply all terms with constant value
-
-    /** 
-     * Pre- and postmultiply all Column blocks with the given
-     * matrices. Hence the new columns will be col_new = pre*col_old*post
-     *
-     * @param pre : The pre-multiply matrix
-     * @param post: The post-multiply matrix
-     */
-    void prePostMultiply(const dmat& pre,const dmat& post);
-
-
-};
-
-  
-
-
-
-#endif // JACROW_H
-//////////////////////////////////////////////////////////////////////

src/sys/segment.h

@@ -15,9 +15,10 @@
 
 #include "tasmet_types.h"
 #include "tasmet_exception.h"
+#include "jacobian.h"
 // #include "phaseconstraint.h"
 
-class Jacobian;
+
 class GlobalConf;
 class TaSystem;
 
@@ -51,9 +52,10 @@ public:
      */    
     virtual Segment* copy(const TaSystem&) const = 0;
     
-    virtual vd initialSolution() const = 0;
+    virtual void initialSolution(SegPositionMapper&) const = 0;
 
-    virtual void residualJac(arma::subview_col<d>&& residual // Here we store the residual
+    virtual void residualJac(SegPositionMapper& residual,
+                             SegJacRows& jacrows// Here we store the residual
                              ) const=0;
 
     // Get and set name

src/sys/tasystem.cpp

@@ -60,47 +60,8 @@ TaSystem::TaSystem(const pb::System& sys):
             throw e;
         }
     }
-    
-    // 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;
+    initSolRes();
-
-
-
-    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();
-            }
-            else {
-                solution.subvec(dofend(i-1),dofend(i)-1) =
-                    seg->initialSolution();
-            }
-            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;
-    // if(size>0) {
-    //     _solution = vd(size);
-    //     for(auto& val: sol) {
-    //         _solution(i) = val;
-    //         i++;
-    //     }
-    // }
 
 }
 TaSystem::TaSystem(const TaSystem& o):
@@ -117,7 +78,89 @@ TaSystem::TaSystem(const TaSystem& o):
         seg.second = seg.second->copy(*this);
         if(!seg.second) throw TaSMETBadAlloc();
     }
+
+    initSolRes();
+
 }
+ void TaSystem::initSolRes() {
+     // Create the initial solution from the segments and store it
+     // here. Please be careful not to call any virtual functions!
+     vus ndofs = getNDofs();
+     vus neqs = getNEqs();
+     
+     us Ns = this->Ns();
+
+     us total_dofs = arma::sum(ndofs);
+     us total_eqs = arma::sum(neqs);
+
+     if(total_dofs != total_eqs) {
+         throw TaSMETError("Number of equations does not match number of DOFS. "
+                           "This is probably due to missing, or too much boundary conditions");
+     }
+
+     if(total_dofs*Ns > constants::maxndofs) {
+        stringstream error;
+        error << "Too many DOFS required. Problem too large. "
+            "Number of equations computed: ";
+        error << total_dofs*Ns;
+        throw TaSMETError(error);
+     }
+
+
+     _solution = zeros(total_dofs*Ns);
+     _residual = zeros(total_dofs*Ns)+arma::datum::nan;
+    
+     d* solptr = _solution.memptr();
+     d* resptr = _residual.memptr();
+     
+     us globalposdof=0;
+     us globalposeq=0;
+
+     // Loop over the segments, and for each dof in the segment, we
+     // allocate a vd, with its memory pointed to from the global
+     // solution pointer.
+
+     for(auto& seg_: _segs) {
+
+         us segid = seg_.first;
+         VARTRACE(15,segid);
+
+         for(us i=0;i< ndofs(segid); i++) {
+            
+             // false: own memory (means the vd does not own the memory
+             // true: strict true means bound to the memory for its
+             // lifetime. Does not allow memory adjustments.
+             _solution_dofs[segid].push_back(new vd(&solptr[globalposdof*Ns],Ns,false,true));
+             globalposdof++;
+             VARTRACE(15,i);
+         }
+
+         // Set the initial solution
+         if(ndofs(segid)>0) 
+             seg_.second->initialSolution(_solution_dofs.at(segid));
+
+
+
+         for(us i=0;i< neqs(segid); i++) {
+            
+             // false: own memory (means the vd does not own the memory
+             // true: strict true means bound to the memory for its
+             // lifetime. Does not allow memory adjustments.
+             _residual_eqs[segid].push_back(new vd(&resptr[globalposeq*Ns],Ns,false,true));
+
+             _jac[segid].push_back(JacRow());
+
+             globalposeq++;
+         }
+
+     }
+
+     tasmet_assert(globalposdof == total_dofs,"Bug in assigning solution vectors");
+     tasmet_assert(globalposeq == total_eqs,"Bug in assigning residual vectors");
+
+ }
+
+
 int TaSystem::getArbitrateMassEq(const vus& neqs) const {
     // Tells the TaSystem which Dof should be overwritten with the
     // mass arbitration equation. The special value of -1 means, that
@@ -156,56 +199,85 @@ int TaSystem::getArbitrateMassEq(const vus& neqs) const {
 void TaSystem::residualJac(ResidualJac& resjac) const {
     TRACE(15,"TaSystem::residual()");
 
-    vus neqs = getNEqs();
+    int arbitrateMassEq = -1;//getArbitrateMassEq(neqs);
-    us total_neqs = arma::sum(neqs);
-
-    if(total_neqs>constants::maxndofs)      {
-        stringstream error;
-        error << "Too many DOFS required. Problem too large. Number of equations computed: ";
-        error << total_neqs;
-        throw TaSMETError(error);
-    }
-
-    // This vector of indices stores the last equation number + 1 for
-    // each equation set in a Segment
-    vus eqsend = arma::cumsum(neqs);
-
-    int arbitrateMassEq = getArbitrateMassEq(neqs);
 
     // This is the mass in the sytem. Only counted when there is an
     // equation present which arbitrates the mass in the system
     d mass=0;
-    
-    VARTRACE(25,total_neqs);
-    VARTRACE(25,eqsend);
-
-    vd residual(total_neqs);
 
     #ifdef TASMET_DEBUG
-    residual = arma::datum::nan*ones(total_neqs);
+    _residual *= arma::datum::nan;
     #endif
 
-    us i=0;
+    us segid;
     const Segment* seg;
     for(auto seg_: _segs) {
-
+        segid = seg_.first;
         seg = seg_.second;
 
-        if(i==0) {
             // Put the residual of the segment in the over-all residual
-            seg->residualJac(residual.subvec(0,eqsend(0)-1));
+        seg->residualJac(_residual_eqs[segid],_jac[segid]);
-        }
-        else {
-            seg->residualJac(residual.subvec(eqsend(i-1),eqsend(i)-1));
-        }
 
         // Count the mass, add it
         if(arbitrateMassEq!=-1) {
             mass += seg->getMass();
         }
         
-        i++;
+     }
-    }
+
+    us Ns = this->Ns();
+    vus neqs = getNEqs();
+    vus neqs_sum = arma::cumsum(neqs);
+    us total_neqs = arma::sum(neqs);
+
+    vus ndofs = getNDofs();
+    vus ndofs_sum = arma::cumsum(ndofs);
+    us total_ndofs = arma::sum(ndofs);
+
+    // Now we create the triplets
+    TripletList triplets(total_neqs*Ns);
+
+    for(auto& segjacrows : _jac) {
+
+        // The current segment number
+        us segno = segjacrows.first;
+
+        // The row offset for this segment
+        us rowoffset = (segno==0) ? 0 : neqs_sum(segno-1)*Ns;
+
+        for(auto& jacrow: segjacrows.second) {
+            
+            for(auto& jaccol : jacrow) {
+                const PosId& id = jaccol.first;
+
+                us coloffset = (id.segno==0) ? 0 : ndofs_sum(segno-1)*Ns;
+                coloffset += id.segoffset;
+                coloffset*=Ns;
+
+                dmat& colmat = jaccol.second;
+                
+                for(us i=0;i<Ns;i++){
+
+                    us row = rowoffset + i;
+
+                    for(us j=0;j<Ns;j++) {
+
+                        us col = coloffset + j;
+                        triplets.addTriplet(Triplet(row,col,colmat(i,j)));
+
+                    }
+
+                    
+                }
+                
+                
+            } // For loop over the columns
+
+            rowoffset += Ns;
+
+        } // For loop over the rows for this segment
+
+    } // For loop over the segments
 
     TRACE(15,"Obtained residual from all Segments");
 
@@ -214,48 +286,33 @@ void TaSystem::residualJac(ResidualJac& resjac) const {
     #endif // TASMET_DEBUG
 
     // Exchange equation if we need to arbitrate mass
-    if(arbitrateMassEq!=-1) {
+    if(arbitrateMassEq>=0) {
-        residual(arbitrateMassEq)=mass - _mass;
+        _residual(arbitrateMassEq)=mass - _mass;
     }
 
-    resjac.residual = residual;
+    resjac.jacobian = sdmat(triplets);
+    resjac.residual = _residual;
+
 }
 vd TaSystem::getSolution() const {
-
     return _solution;
 }
-const arma::subview_col<d> TaSystem::getSolution(const us seg_id) const {
+const SegPositionMapper& TaSystem::getSolution(const us seg_id) const {
-
+    return _solution_dofs.at(seg_id);
-    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 {
-        return _solution.subvec(dofsend(seg_id-1),dofsend(seg_id)-1);
-    }
 }
 vus TaSystem::getNDofs() const  {
     TRACE(0,"TaSystem::getNDofs()");
     vus Ndofs(_segs.size());
-    us i=0;
     for (auto seg : _segs) {
-        Ndofs(i)=seg.second->getNDofs();
+        Ndofs(seg.first)=seg.second->getNDofs();
-        i++;
     }
     return Ndofs;
 }
 vus TaSystem::getNEqs() const  {
     TRACE(15,"TaSystem::getNEqs()");
     vus Neqs(_segs.size());
-    us i=0;
     for (auto seg :_segs) {
-        Neqs(i)=seg.second->getNEqs();
+        Neqs(seg.first)=seg.second->getNEqs();
-        VARTRACE(15,Neqs(i));
-        i++;
     }
     return Neqs;
 }
@@ -341,7 +398,19 @@ TaSystem::~TaSystem() {
     cleanup();
 }
 void TaSystem::cleanup() {
+
     purge(_segs);
+    for(auto& segdofs: _solution_dofs) {
+        for(vd* dof : segdofs.second){
+            delete dof;
+        }
+    }
+    for(auto& segdofs: _residual_eqs) {
+        for(vd* dof : segdofs.second){
+            delete dof;
+        }
+    }
+
 }
 
 //////////////////////////////////////////////////////////////////////

src/sys/tasystem.h

@@ -11,9 +11,9 @@
 #include "system.h"
 #include "globalconf.h"
 #include "triplets.h"
+#include "jacobian.h"
 #include "gas.h"
 #include "protobuf/system.pb.h"
-#include <map>
 #include <memory>
 
 class Segment;
@@ -28,21 +28,48 @@ protected:
 
     std::unique_ptr<Gas> _gas;
 
+    mutable
     vd _solution;               /**< This column vector stores the
                                    current solution. */
 
-    std::vector<vd> _dofs;      /**< This vector of column vectors
+    mutable
-                                   stores columns of doubles that do
+    GlobalPositionMapper _solution_dofs; /**< This vector of column vectors
-                                   not own the memory they point
+                                                      stores columns of doubles that do
-                                   to. Instead, they use the memory of
+                                                      not own the memory they point
-                                   _solution.*/
+                                                      to. Instead, they use the memory of
+                                                      _solution.*/
 
-    vus _startdof;              // Store the start DOFS for each segment
+    
+    mutable
+    vd _residual;
+
+    mutable
+    GlobalPositionMapper _residual_eqs;
+
+    mutable
+    Jacobian _jac;              /**< Storing the Jacobian brings the
+                                   hugh benefit of only allocating the
+                                   memory once. Once this is done, all
+                                   values can be updated during each
+                                   iteration. */
     
     TaSystem& operator=(const TaSystem& other)=delete;
+
+private:
+    /** 
+     * Initializes the _solution and _solution_dofs vector. This
+     * should be done during construction of the object. The method is
+     * called from all constructors. Should be done after construction
+     * of the segments.
+     *
+     */
+    void initSolRes();
 public:
     TaSystem(const TaSystem& o);
     TaSystem(const pb::System&);
+
+
+
     static TaSystem testSystem() {
         return TaSystem(pb::System::default_instance());
     }
@@ -67,7 +94,7 @@ public:
     vd getSolution() const;
 
     // Obtain the solution vector for the Segment with given id
-    const arma::subview_col<d> getSolution(const us seg_id) const;
+    const SegPositionMapper& getSolution(const us seg_id) const;
     
     virtual void updateSolution(const vd& sol) {_solution = sol; }	// Update the solution