diff --git a/CMakeLists.txt b/CMakeLists.txt index 434c862..5a04da5 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b8f55d1..66842d8 100644 --- a/src/CMakeLists.txt +++ b/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 ) diff --git a/src/solver/bracket_root.h b/src/solver/bracket_root.h index 82a3526..4ff86a6 100644 --- a/src/solver/bracket_root.h +++ b/src/solver/bracket_root.h @@ -24,7 +24,13 @@ inline std::pair bracket_root(NoGradientNonlinearSystem& 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::epsilon(); diff --git a/src/solver/brent.cpp b/src/solver/brent.cpp index cd88330..064384f 100644 --- a/src/solver/brent.cpp +++ b/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& 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"); diff --git a/src/sys/segment.h b/src/sys/segment.h new file mode 100644 index 0000000..22a72f1 --- /dev/null +++ b/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 +////////////////////////////////////////////////////////////////////// diff --git a/src/sys/tasystem.cpp b/src/sys/tasystem.cpp new file mode 100644 index 0000000..421d1a4 --- /dev/null +++ b/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;kresetHarmonics(); + } +} +dmat TaSystem::showJac(){ + + TRACE(15,"TaSystem::showJac()"); + + return dmat(jacobian()); +} + +TaSystem::~TaSystem() { + TRACE(25,"~TaSystem()"); + + for(auto& seg: _segs){ + delete seg.second; + } + delete _gas; +} + + +////////////////////////////////////////////////////////////////////// diff --git a/src/sys/tasystem.h b/src/sys/tasystem.h new file mode 100644 index 0000000..6ce8bc1 --- /dev/null +++ b/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 + +// Inherit all global configuration members +class TaSystem : public GradientNonlinearSystem { +protected: + + d _mass = -1; + + std::map _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_ */ + diff --git a/src/tasmet_exception.h b/src/tasmet_exception.h index 44a455c..cec5052 100644 --- a/src/tasmet_exception.h +++ b/src/tasmet_exception.h @@ -11,10 +11,12 @@ #include #include +#include // 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 { diff --git a/src/tasmet_types.h b/src/tasmet_types.h index aa276dd..fdb21bf 100644 --- a/src/tasmet_types.h +++ b/src/tasmet_types.h @@ -74,6 +74,7 @@ using vci = arma::Col::fixed; typedef arma::Col vd; /* Column vector of doubles */ typedef arma::Col vc; /* Column vector of complex numbers */ +typedef arma::Col vus; /* Column vector of unsigned integers */ typedef arma::Mat dmat; /* (Dense) Matrix of doubles */ typedef arma::Mat cmat; /* (Dense) matrix of complex numbers */ diff --git a/testing/solver/brent_test.cpp b/testing/solver/brent_test.cpp index 57f5fce..06ab724 100644 --- a/testing/solver/brent_test.cpp +++ b/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); }