Added improvements on solver, added TaSystem in new form

CMakeLists.txt

@@ -19,7 +19,7 @@ add_definitions(-DTASMET_DEBUG=1)
 #====================================================
 
 # Always required make flags in case of both compilers
-set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -pipe -fPIC -Wall -Wextra -Wno-unused-parameter -Wno-unused-variable  ")
+set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -pipe -fPIC -Wfatal-errors -Wall -Wextra -Wno-unused-parameter -Wno-unused-variable  ")
 
 
 # Stop letting Numpy complain about its API

src/CMakeLists.txt

@@ -53,7 +53,7 @@ add_library(tasmet_src
   sys/jaccol.cpp
   sys/jacobian.cpp
   sys/jacrow.cpp
-  # sys/tasystem.cpp
+  sys/tasystem.cpp
   sys/triplets.cpp
 )
  

src/solver/bracket_root.h

@@ -24,7 +24,13 @@ inline std::pair<d,d> bracket_root(NoGradientNonlinearSystem<d>& sys,const d gue
     d xa = guess;
     d fa = sys.residual(xa);
 
+    if(fa == 0) {
+        // Guess was exactly right
+        return std::make_pair(xa,xa);
+    }
+    
     d fb = fa;
+
     // We add here a small amount, such that, in case the first guess
     // is 0, we 
     d xb = xa+std::numeric_limits<d>::epsilon();

src/solver/brent.cpp

@@ -6,7 +6,7 @@
 // Brent's root finding algorithm as implemented from Wikipedia:
 // https://en.wikipedia.org/wiki/Brent's_method
 //////////////////////////////////////////////////////////////////////
-
+#define TRACERPLUS (10)
 #include "brent.h"
 #include "tasmet_tracer.h"
 #include "tasmet_exception.h"
@@ -56,6 +56,11 @@ void Brent::start_implementation(NoGradientNonlinearSystem<d>& system,
     d b = brackets.second;
 
     d fa = system.residual(a);    
+    if((fa) == 0) {
+        TRACE(15,"Found root during bracketing");
+        return;
+    }
+    
     d fb = system.residual(b);
     if((fb) == 0) {
         TRACE(15,"Found root during bracketing");

src/sys/segment.h

@@ -0,0 +1,79 @@
+// segment.h
+//
+// Author: J.A. de Jong 
+//
+// Description: The Segments class is the base class for
+// both the normal segments and the connectors. The main difference
+// between a segment and a connector is, that a connector only provide
+// equations. A normal segment contains both. All
+// common stuff, such as a name, initialization
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef SEGMENT_H
+#define SEGMENT_H
+
+#include "tasmet_types.h"
+#include "tasmet_exception.h"
+// #include "phaseconstraint.h"
+
+class Jacobian;
+class GlobalConf;
+class TaSystem;
+
+class Segment{
+
+    us _id; // Unique number for each segment in a
+    // TaSystem
+
+    // User identifier
+    std::string _name;
+
+protected:
+    Segment(us id,const std::string& name): _id(id),_name(name) {}
+    Segment(const Segment& o): Segment(o._id,o._name){}
+    Segment& operator=(const Segment&)=delete;
+public:
+    virtual ~Segment(){}
+
+    virtual Segment* copy() const = 0;
+    
+    // Get and set name
+    const std::string& getName() const{return _name;} // This one is just the name
+    us getid() const{return _id;} // This one is just the name
+    void setName(const std::string& name){ _name = name; } // This one is just the name
+    void setid(const us id){ _id=id;} // Set ID
+    
+    // Tell a TaSystem whether this Segment arbitrates Mass or
+    // not. The special return value of -1 tells it does not. If it
+    // does, the derived class should return which equation should be
+    // overwritten with the mass arbitration equation.
+    virtual int arbitrateMassEq() const {return -1;}
+    virtual void residual(const TaSystem&,vd&,const us insertion_start) const=0;
+
+    virtual void updateSolution(const TaSystem&,const vd&) const = 0;
+    virtual vd getSolution(const TaSystem&,vd& sol,const us insertion_start) const = 0;
+
+    // Return the total number of equations in this segment
+    virtual us getNEqs(const TaSystem&) const { return 0;}
+    // Return the total number of DOFS in this segment
+    virtual us getNDofs(const TaSystem&) const { return 0;}
+
+    // Return the current mass in this segment
+    virtual d getMass(const TaSystem&) const = 0;
+    virtual void dmtotdx(const TaSystem&,vd& dmtotdx,us dof_start) const {}
+    
+    virtual void show(const TaSystem&,us verbosity_level) const=0;
+
+    // Reset amplitude data in higher harmonics
+    void resetHarmonics();
+
+    // Fill Jacobian with values from the equations in this
+    // segment/connector.
+    virtual void jac(const TaSystem&,Jacobian&,us dof_start,us eq_start) const=0;
+
+};
+
+
+#endif // SEGMENT_H
+//////////////////////////////////////////////////////////////////////

src/sys/tasystem.cpp

@@ -0,0 +1,297 @@
+// tasystem.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#include "segment.h"
+#include "tasystem.h"
+#include "triplets.h"
+#include "jacobian.h"
+#include "tasmet_assert.h"
+#include "tasmet_exception.h"
+#include "tasmet_constants.h"
+
+TaSystem::TaSystem(const GlobalConf& gc,Gas::GasType gastype):
+    _gc(new GlobalConf(gc)),
+    _gas(Gas::newGas(gastype))
+{
+    TRACE(14,"TaSystem::TaSystem(gc,gastype)");
+}
+TaSystem::TaSystem(const TaSystem& o):
+    _gc(o._gc),                 // Share a ptr to the Global conf
+    _gas(Gas::newGas(Gas::GasType(*o._gas)))
+{
+    TRACE(25,"TaSystem::TaSystem(TaSystem&) copy");
+
+    // Take over the pointers to the segments
+    _segs = o._segs;
+
+    for(auto& seg: _segs){
+        seg.second = seg.second->copy();
+    }
+
+}
+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
+    // mass is not arbitrated. This can be the case if for example a
+    // pressure boundary condition is present.
+    int arbitrateMassEq=-1;
+    bool masseq_found = false;
+
+    us total_neqs = 0, i=0;
+
+    // Iterate over all segments to find a single one segment Dof that is
+    // willing to arbitrate the mass in the system. This can, of
+    // course, be only be one segment.
+    int local_mass_eq = -1;
+    
+    const Segment* seg;
+    for(auto& seg_: _segs) {
+        seg = seg_.second;
+
+        if(masseq_found && (local_mass_eq = seg->arbitrateMassEq())) {
+            throw TaSMETError("Error: only one Segment is allowed to"
+                              " arbitrate the system mass");
+        }
+        else if(local_mass_eq != -1){
+
+            arbitrateMassEq = total_neqs+local_mass_eq;
+
+            masseq_found = true;
+        }
+
+        total_neqs += neqs(i);
+        i++;
+    }
+    return arbitrateMassEq;
+}
+TaSystem& TaSystem::operator+=(const Segment& s){
+
+    TRACE(24,"TaSystem::operator+=(Seg)");
+
+    if(_segs.find(s.getid())!=_segs.end()){
+        std::stringstream error;
+        error << "Segment with id " << s.getid() <<
+            "already present in the system";
+        throw TaSMETError(error);
+    }
+    _segs[s.getid()]=s.copy();
+    return *this;
+}
+void TaSystem::updateSolution(const vd& sol) {
+    TRACE(15,"TaSystem::updateSolution()");
+
+    us firstdof = 0;
+    us ndofs;
+    Segment* seg;
+    for(auto& seg_: _segs) {
+        seg = seg_.second;
+        ndofs = seg->getNDofs(*this);
+
+        firstdof += ndofs;
+
+        seg->updateSolution(*this,sol.subvec(firstdof,firstdof+ndofs-1));
+    }
+}
+vd TaSystem::getSolution() const {
+    us firstdof = 0,i=0;
+    vus ndofs = getNDofs();
+    us total_ndofs = arma::sum(ndofs);
+    vd sol(total_ndofs);
+
+    #ifdef TASMET_DEBUG
+    sol.zeros();
+    #endif
+
+    Segment* seg;
+    for(auto& seg_: _segs) {
+        seg = seg_.second;
+        seg->getSolution(*this,sol,firstdof);
+        firstdof += ndofs(i);
+        i++;
+    }
+    return sol;
+}
+vd TaSystem::residual() const {
+    TRACE(15,"TaSystem::residual()");
+
+    vus neqs = getNEqs();
+    us total_neqs = arma::sum(neqs);
+    vus eqstart(neqs.size());
+
+    us i=0;
+
+    for(us& eqs : neqs) {
+
+        if(i>0) {
+            eqstart(i) = eqstart(i-1) + eqs;
+        }
+        else {
+            eqstart(i) = 0;
+        }
+        i++;                 
+    }
+    assert(i==_segs.size()-1);
+
+    if(total_neqs>constants::maxndofs)      {
+        throw TaSMETError("Too many DOFS required."
+                          " Problem too large.");
+    }
+
+    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,eqstart);
+
+    vd residual(total_neqs);
+
+    i=0;
+    const Segment* seg;
+    for(auto seg_: _segs) {
+
+        seg = seg_.second;
+
+        // Put the residual of the segment in the over-all residual
+        seg->residual(*this,residual,eqstart(i));
+
+        // Count the mass, add it
+        if(arbitrateMassEq!=-1) {
+            mass += seg->getMass(*this);
+        }
+        
+        i++;
+    }
+
+    #ifdef TASMET_DEBUG
+    assert(arbitrateMassEq< (int) total_neqs);    
+    #endif // TASMET_DEBUG
+
+    // Exchange equation if we need to arbitrate mass
+    if(arbitrateMassEq!=-1) {
+        residual(arbitrateMassEq)=mass - _mass;
+    }
+
+    return residual;
+}
+
+vus TaSystem::getNDofs() const  {
+    TRACE(0,"TaSystem::getNDofs()");
+    vus Ndofs(_segs.size());
+    us i=0;
+    for (auto seg : _segs) {
+        Ndofs(i)=seg.second->getNDofs(*this);
+        i++;
+    }
+    return Ndofs;
+}
+vus TaSystem::getNEqs() const  {
+    TRACE(0,"TaSystem::getNDofs()");
+    vus Neqs(_segs.size());
+    us i=0;
+    for (auto seg :_segs) {
+        Neqs(i)=seg.second->getNEqs(*this);
+        i++;
+    }
+    return Neqs;
+}
+void TaSystem::show(us detailnr){
+    TRACE(25,"TaSystem::show()");
+
+    cout << "########################## Showing TaSystem...\n";
+    cout << "Showing Global configuration...\n";
+    _gc->show();
+    
+    if(detailnr>0){
+        for(auto seg:_segs){
+            cout << "Showing segment with ID " << seg.first << "\n";
+            seg.second->show(*this,detailnr);
+        }
+    } // detailnr>0
+}
+
+TripletList TaSystem::jacTriplets() const {
+
+    TRACE(14,"TaSystem::jacobian()");
+    vus ndofs=getNDofs();
+    vus neqs=getNEqs();
+    us total_ndofs = arma::sum(ndofs);
+    us total_neqs = arma::sum(neqs);
+    us i=0,dofnr=0,eqnr=0;
+    Jacobian j(total_ndofs);
+
+    int arbitrateMassEq = getArbitrateMassEq(neqs);
+    vd dmtotdx;
+    if(arbitrateMassEq > -1)
+        dmtotdx = vd(total_ndofs);
+
+    const Segment* seg;
+    for(auto& seg_ :_segs) {
+        seg = seg_.second;
+        seg->jac(*this,j,dofnr,eqnr);
+        if(arbitrateMassEq > -1) {
+            seg->dmtotdx(*this,dmtotdx,dofnr);
+        }
+        eqnr += neqs(i);
+        dofnr += ndofs(i);
+        i++;
+    }
+
+    if(dofnr!=eqnr){
+        throw TaSMETError("System of equations is over/underdetermined");
+    }
+    
+    // Convert to tripletlist
+    TripletList jac=j;
+
+    assert(arbitrateMassEq< (int) total_neqs);
+
+    // Exchange equation if we need to arbitrate mass
+    if(arbitrateMassEq!=-1) {
+        // Replace this equation with global mass conservation
+        jac.zeroOutRow(arbitrateMassEq);
+
+        us dmtotdxsize=dmtotdx.size();
+        for(us k=0;k<dmtotdxsize;k++)
+            if(dmtotdx(k)!=0){
+                // TRACE(20,"k: " << k);
+                // TRACE(20,"dmtotdx:"<< dmtotdx(k));
+                jac.addTriplet(Triplet(arbitrateMassEq,k,dmtotdx(k)));
+            }
+    }
+
+    return jac;
+}
+sdmat TaSystem::jacobian() const {
+    TRACE(14,"TaSystem::Jac()");
+    return jacTriplets();            // Implicitly converts to sdmat
+}
+void TaSystem::resetHarmonics(){
+    for(auto seg: _segs) {
+        seg.second->resetHarmonics();
+    }
+}
+dmat TaSystem::showJac(){
+
+    TRACE(15,"TaSystem::showJac()");
+
+    return dmat(jacobian());
+}
+
+TaSystem::~TaSystem() {
+    TRACE(25,"~TaSystem()");
+
+    for(auto& seg: _segs){
+        delete seg.second;
+    }
+    delete _gas;
+}
+
+
+//////////////////////////////////////////////////////////////////////

src/sys/tasystem.h

@@ -0,0 +1,89 @@
+// tasystem.h, created March 19th 2014
+// Author: J.A. de Jong
+
+// A system contains one or more (un)connected (thermoacoustic)
+// segments. Some segments require boundary conditions, these boundary
+// conditions are provided as segments as well.
+
+#pragma once
+#ifndef _TASYSTEM_H_
+#define _TASYSTEM_H_
+#include "system.h"
+#include "segment.h"
+#include "globalconf.h"
+#include "triplets.h"
+#include "gas.h"
+#include <map>
+
+// Inherit all global configuration members
+class TaSystem : public GradientNonlinearSystem {
+protected:
+
+    d _mass = -1;
+    
+    std::map<us,Segment*> _segs;		
+
+    gc_ptr _gc;
+
+    Gas* _gas = nullptr;
+    
+    TaSystem& operator=(const TaSystem& other)=delete;
+    TaSystem(const TaSystem& o);
+public:
+    TaSystem(const GlobalConf& g,Gas::GasType gastype);
+
+    // Set globalconf configuration. Applies updateNf as well.
+    void setGc(const GlobalConf& gc);
+    const Gas& getGas() const {return *_gas;}
+    void setGas(Gas::GasType gastype);
+    
+    // Set and get the mass in the system. If the mass is not set
+    // before initializing, the mass is computed from the segment's
+    // intial configuration.
+    void setMass(d mass) { _mass = mass; }
+    d getMass() const;
+    
+    virtual ~TaSystem();
+    virtual TaSystem* copy() const {return new TaSystem(*this);}
+
+    us nSegments() const {return _segs.size();}
+
+    TaSystem& operator+=(const Segment& s);	// Add a segment to the
+    // system. It creates a copy
+
+    dmat showJac();
+    virtual void show(us detailnr=0);
+
+    vd residual() const; 
+
+    vd getSolution() const;
+    
+    virtual void updateSolution(const vd& resvec);	// Update the solution
+
+    // Compute Jacobian matrix. The dampfac value is used in an
+    // EngineSystem
+    sdmat jacobian() const;		// Return Jacobian matrix
+
+    // Change Nf in the system, while keeping the results.
+    void updateNf(us);
+
+    // Reset amplitude data in higher harmonics
+    void resetHarmonics();
+
+    // void delseg(us n); // Not yet implemented.  Delete a segment
+    // from the system (we have to determine how elaborated the API
+    // has to be.)
+    vus getNDofs() const;	// Compute DOFS in all segments
+    vus getNEqs() const;        // Compute Equations in all segments
+
+    const Segment& getSegment(const us id) const {return *_segs.at(id);}
+
+protected:
+    virtual int getArbitrateMassEq(const vus& neqs) const;
+    virtual TripletList jacTriplets() const;
+
+}; // class System
+  
+
+#endif /* _TASYSTEM_H_ */
+

src/tasmet_exception.h

@@ -11,10 +11,12 @@
 
 #include <string>
 #include <stdexcept>
+#include <sstream>              // stringstream
 
 class TaSMETError : public std::runtime_error {
 public:
     TaSMETError(const std::string& msg = "") : std::runtime_error(msg) {}
+    TaSMETError(const std::stringstream& stream) : std::runtime_error(stream.str()){}
 };
 
 class TasMETBadAlloc: public std::bad_alloc {

src/tasmet_types.h

@@ -74,6 +74,7 @@ using vci = arma::Col<c>::fixed<i>;
 
 typedef arma::Col<d> vd;		  /* Column vector of doubles */		
 typedef arma::Col<c> vc;	  /* Column vector of complex numbers */	
+typedef arma::Col<us> vus;	  /* Column vector of unsigned integers */	
 typedef arma::Mat<d> dmat; /* (Dense) Matrix of doubles */		
 typedef arma::Mat<c> cmat;	  /* (Dense) matrix of complex numbers */	
 

testing/solver/brent_test.cpp

@@ -63,11 +63,11 @@ int main(){
 
     Nitrogen nit(293.15,101325);
     Air air(293.15,101325);
-
+    Helium helium(293.15,101325);
 
     gc_ptr gc(new GlobalConf(10,100));
     Variable pressure(gc,10*101325);
-    Variable temperature = adiabaticTemp(air,pressure);
+    Variable temperature = adiabaticTemp(helium,pressure);
 }