diff --git a/src/duct/duct.cpp b/src/duct/duct.cpp index 2081fea..8c08eb9 100644 --- a/src/duct/duct.cpp +++ b/src/duct/duct.cpp @@ -11,8 +11,10 @@ #include "tasmet_evalscript.h" #include "perfectgas.h" -Duct::Duct(const us id,const pb::Duct& ductpb): - Segment(id,ductpb.name()), +Duct::Duct(const TaSystem& sys, + const us id, + const pb::Duct& ductpb): + Segment(sys,id,ductpb.name()), Geom(ductpb), _ductpb(ductpb) { @@ -38,8 +40,8 @@ Duct::Duct(const us id,const pb::Duct& ductpb): } } -Duct::Duct(const Duct& other): - Segment(other), +Duct::Duct(const TaSystem& sys,const Duct& other): + Segment(sys,other), Geom(other), _ductpb(other._ductpb), _Tsprescribed(other._Tsprescribed) @@ -47,11 +49,10 @@ Duct::Duct(const Duct& other): // Do something with the equations here TRACE(15,"Duct::~Duct"); - } -Duct* Duct::copy() const { - return new Duct(*this); +Duct* Duct::copy(const TaSystem& sys) const { + return new Duct(sys,*this); } Duct::~Duct() { TRACE(15,"Duct::~Duct"); @@ -59,7 +60,7 @@ Duct::~Duct() { // delete eq; // } } -void Duct::residual(const TaSystem& sys,arma::subview_col && residual) const { +void Duct::residualJac(arma::subview_col && residual) const { TRACE(15,"Duct::residual()"); @@ -86,10 +87,10 @@ void Duct::residual(const TaSystem& sys,arma::subview_col && residual) const 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); + rhop = getvart(constants::rho,0); + up = getvart(constants::u,0); + Tp = getvart(constants::T,0); + pp = getvart(constants::p,0); const Gas& gas = sys.gas(); @@ -105,10 +106,10 @@ void Duct::residual(const TaSystem& sys,arma::subview_col && residual) const 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); + rhop = getvart(constants::rho,gp+1); + up = getvart(constants::u,gp+1); + Tp = getvart(constants::T,gp+1); + pp = getvart(constants::p,gp+1); cont = ((rhop%up)-(rho%u))/dx; @@ -162,7 +163,7 @@ void Duct::residual(const TaSystem& sys,arma::subview_col && residual) const } } -vd Duct::getvart(const TaSystem& sys,int varnr,int gp) const { +vd Duct::getvart(int varnr,int gp) const { TRACE(15,"Duct::getvart()"); const arma::subview_col sol = sys.getSolution(_id); @@ -177,18 +178,18 @@ vd Duct::getvart(const TaSystem& sys,int varnr,int gp) const { 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 Duct::getvarx(int varnr,int t) const { vd res(ngp()); for(us i=0;i _eqs; + pb::Duct _ductpb; vd _Tsprescribed; - - protected: - Duct(const Duct&); + Duct(const TaSystem&,const Duct&); public: - Duct(const us id,const pb::Duct&); + Duct(const TaSystem&,const us id,const pb::Duct&); ~Duct(); vd Tsprescribed() const { return _Tsprescribed;} @@ -39,58 +37,60 @@ public: // Initialize the solution to something sensible vd initializeSolution(const TaSystem& sys); - virtual Duct* copy() const; + virtual Duct* copy(const TaSystem&) const; const Geom& geom() const; 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); } + /** + * Returns the dof nr for the given variable at given gridpoint. + * + * @return The index of the degree of freedom + */ + Dof getDof(int varnr,int gp) const; - 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); } + /// Postprocessing functions + vd rhox(int t) const { return getvarx(constants::rho,t);} + vd ux(int t) const { return getvarx(constants::u,t); } + vd Tx(int t) const { return getvarx(constants::T,t); } + vd px(int t) const { return getvarx(constants::p,t); } + vd Tsx(int t) const { return getvarx(constants::Ts,t); } + + vd rhot(int gp) const { return getvart(constants::rho,gp); } + vd ut(int gp) const { return getvart(constants::u,gp); } + vd Tt(int gp) const { return getvart(constants::T,gp); } + vd pt(int gp) const { return getvart(constants::p,gp); } + vd Tst(int gp) const { return getvart(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; + vd getvart(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; + vd getvarx(int varnr,int t) const; - d getvartx(const TaSystem& sys,int t,int gp) const; + d getvartx(int t,int gp) const; // Solving - virtual void residual(const TaSystem&,arma::subview_col&& residual) const; + virtual void residualJac(arma::subview_col&& residual) const; - vd initialSolution(const TaSystem&) const; + vd initialSolution() const; // virtual void getSolution(const TaSystem&,const us insertion_start) const; // Return the total number of equations in this segment - virtual us getNEqs(const TaSystem&) const; + virtual us getNEqs() const; // Return the total number of DOFS in this segment - virtual us getNDofs(const TaSystem&) const; + virtual us getNDofs() const; // Return the current mass in this segment - virtual d getMass(const TaSystem&) const; - virtual void dmtotdx(const TaSystem&,vd& dmtotdx,us dof_start) const {} + virtual d getMass() const; + virtual void dmtotdx(vd& dmtotdx,us dof_start) const {} - virtual void show(const TaSystem&,us verbosity_level) const; + virtual void show(us verbosity_level) const; // 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; - }; #endif // DUCT_H diff --git a/src/ductbc/adiabaticwall.cpp b/src/ductbc/adiabaticwall.cpp index b23437b..10c5b88 100644 --- a/src/ductbc/adiabaticwall.cpp +++ b/src/ductbc/adiabaticwall.cpp @@ -12,19 +12,19 @@ #include "tasystem.h" #include "duct.h" -AdiabaticWall::AdiabaticWall(const us id, - const TaSystem& sys, - const pb::DuctBc& dbc): - DuctBc(id,dbc) +AdiabaticWall::AdiabaticWall(const TaSystem& sys, + const us id, + const pb::DuctBc& dbc): + DuctBc(sys,id,dbc) { TRACE(15,"AdiabaticWall(id,sys,dbc)"); tasmet_assert(dbc.type() == pb::AdiabaticWall,"Wrong type given to constructor"); - } -AdiabaticWall::AdiabaticWall(const AdiabaticWall& o): - DuctBc(o) +AdiabaticWall::AdiabaticWall(const TaSystem& sys, + const AdiabaticWall& o): + DuctBc(sys,o) { TRACE(15,"AdiabaticWall(o)"); @@ -32,15 +32,14 @@ AdiabaticWall::AdiabaticWall(const AdiabaticWall& o): AdiabaticWall::~AdiabaticWall() { } -AdiabaticWall* AdiabaticWall::copy() const { - return new AdiabaticWall(*this); +AdiabaticWall* AdiabaticWall::copy(const TaSystem& sys) const { + return new AdiabaticWall(sys,*this); } -vd AdiabaticWall::initialSolution(const TaSystem& sys) const { +vd AdiabaticWall::initialSolution() const { return vd(); } -void AdiabaticWall::residual(const TaSystem& sys, - arma::subview_col&& residual +void AdiabaticWall::residualJac(arma::subview_col&& residual ) const { TRACE(15,"AdiabaticWall::residual()"); @@ -67,7 +66,7 @@ void AdiabaticWall::residual(const TaSystem& sys, } -us AdiabaticWall::getNEqs(const TaSystem& sys) const { +us AdiabaticWall::getNEqs() const { TRACE(15,"AdiabaticWall::getNEqs()"); // u = 0 // dT/dx = 0 --> if htmodel is not Isentropic @@ -83,13 +82,8 @@ us AdiabaticWall::getNEqs(const TaSystem& sys) const { VARTRACE(15,neqs); return neqs; } -void AdiabaticWall::show(const TaSystem&,us verbosity_level) const { +void AdiabaticWall::show(us verbosity_level) const { TRACE(15,"AdiabaticWall::show()"); } -void AdiabaticWall::jac(const TaSystem&,Jacobian&, - us dof_start,us eq_start) const { - TRACE(15,"AdiabaticWall::jac()"); - -} ////////////////////////////////////////////////////////////////////// diff --git a/src/ductbc/adiabaticwall.h b/src/ductbc/adiabaticwall.h index 63bd04f..648ae97 100644 --- a/src/ductbc/adiabaticwall.h +++ b/src/ductbc/adiabaticwall.h @@ -24,27 +24,28 @@ class AdiabaticWall: public DuctBc { us _duct_id; /**< ID of Duct for this b.c. */ protected: - AdiabaticWall(const AdiabaticWall&); + AdiabaticWall(const TaSystem& sys, + const AdiabaticWall&); public: - AdiabaticWall(const us id, - const TaSystem& sys, - const pb::DuctBc&); + AdiabaticWall(const TaSystem& sys, + const us id, + const pb::DuctBc&); + ~AdiabaticWall(); - AdiabaticWall* copy() const; + AdiabaticWall* copy(const TaSystem& sys) const; - vd initialSolution(const TaSystem&) const; + vd initialSolution() const; - virtual void residual(const TaSystem&, - arma::subview_col&& residual + virtual void residualJac(arma::subview_col&& residual ) const; // Return the total number of equations in this segment - virtual us getNEqs(const TaSystem&) const; + virtual us getNEqs() const; // Return the current mass in this segment - virtual d getMass(const TaSystem&) const { return 0;}; + virtual d getMass() const { return 0;}; - virtual void show(const TaSystem&,us verbosity_level) const; + virtual void show(us verbosity_level) const; // Reset amplitude data in higher harmonics // virtual void resetHarmonics() = 0; diff --git a/src/ductbc/ductbc.cpp b/src/ductbc/ductbc.cpp index 5ebd4fa..deaa74f 100644 --- a/src/ductbc/ductbc.cpp +++ b/src/ductbc/ductbc.cpp @@ -14,7 +14,7 @@ #include "tasystem.h" #include "duct.h" -const Duct& DuctBc::getDuct(const TaSystem& sys) const { +const Duct& DuctBc::getDuct() const { return dynamic_cast(sys.getSegment(_dbc.duct_id())); @@ -28,13 +28,13 @@ DuctBc* DuctBc::newDuctBc(const us id, switch (dbc.type()) { case pb::PressureBc: { - return new PressureBc(id,sys,dbc); - break; - } - case pb::AdiabaticWall: { - return new AdiabaticWall(id,sys,dbc); + return new PressureBc(sys,id,dbc); break; } + // case pb::AdiabaticWall: { + // return new AdiabaticWall(id,sys,dbc); + // break; + // } default: tasmet_assert(false,"Not implemented DuctBc"); break; diff --git a/src/ductbc/ductbc.h b/src/ductbc/ductbc.h index d700f87..41567ed 100644 --- a/src/ductbc/ductbc.h +++ b/src/ductbc/ductbc.h @@ -17,18 +17,21 @@ class Duct; class DuctBc :public Segment { pb::DuctBc _dbc; public: - DuctBc(const us id, + DuctBc(const TaSystem& sys, + const us id, const pb::DuctBc& dbc): - Segment(id,dbc.name()), + Segment(sys,id,dbc.name()), _dbc(dbc) {} - DuctBc(const DuctBc& o): Segment(o),_dbc(o._dbc) {} + DuctBc(const TaSystem&sys,const DuctBc& o): + Segment(sys,o), + _dbc(o._dbc) {} static DuctBc* newDuctBc(const us id, - const TaSystem& sys, - const pb::DuctBc&); + const TaSystem& sys, + const pb::DuctBc&); - const Duct& getDuct(const TaSystem&) const; + const Duct& getDuct() const; }; diff --git a/src/ductbc/pressurebc.cpp b/src/ductbc/pressurebc.cpp index 3ab6cf8..0ea7ab8 100644 --- a/src/ductbc/pressurebc.cpp +++ b/src/ductbc/pressurebc.cpp @@ -15,10 +15,10 @@ #include "perfectgas.h" #include "adiabatictemp.h" -PressureBc::PressureBc(const us id, - const TaSystem& sys, +PressureBc::PressureBc(const TaSystem& sys, + const us id, const pb::DuctBc& dbc): - DuctBc(id,dbc) + DuctBc(sys,id,dbc) { TRACE(15,"PressureBc(id,sys,dbc)"); vd time = sys.timeInstances(); @@ -54,50 +54,51 @@ PressureBc::PressureBc(const us id, } } -PressureBc::PressureBc(const PressureBc& o): - DuctBc(o),_p(o._p),_T(o._T) { +PressureBc::PressureBc(const TaSystem& sys, + const PressureBc& o): + DuctBc(sys,o), + _p(o._p), + _T(o._T) { TRACE(15,"PressureBc(o)"); - } PressureBc::~PressureBc() { } -PressureBc* PressureBc::copy() const { - return new PressureBc(*this); +PressureBc* PressureBc::copy(const TaSystem& sys) const { + return new PressureBc(sys,*this); } -vd PressureBc::initialSolution(const TaSystem& sys) const { +vd PressureBc::initialSolution() const { return vd(); } -void PressureBc::residual(const TaSystem& sys, - arma::subview_col&& residual +void PressureBc::residualJac(arma::subview_col&& residual ) const { TRACE(15,"PressureBc::residual()"); - const pb::Duct& dpb = getDuct(sys).getDuctPb(); + const pb::Duct& dpb = getDuct().getDuctPb(); us Ns = sys.Ns(); - const Duct& duct = getDuct(sys); + const Duct& duct = getDuct(); if(_side == pb::left) { - residual.subvec(0,Ns-1) = duct.pt(sys,0) - _p; + residual.subvec(0,Ns-1) = duct.pt(0) - _p; if(dpb.htmodel() != pb::Isentropic) { - residual.subvec(Ns,2*Ns-1) = duct.Tt(sys,0) - _T; + residual.subvec(Ns,2*Ns-1) = duct.Tt(0) - _T; } } else { - residual.subvec(0,Ns-1) = duct.pt(sys,-1) - _p; + residual.subvec(0,Ns-1) = duct.pt(-1) - _p; if(dpb.htmodel() != pb::Isentropic) { - residual.subvec(Ns,2*Ns-1) = duct.Tt(sys,-1) - _T; + residual.subvec(Ns,2*Ns-1) = duct.Tt(-1) - _T; } } } -us PressureBc::getNEqs(const TaSystem& sys) const { +us PressureBc::getNEqs() const { TRACE(15,"PressureBc::getNEqs()"); // We provide equations for @@ -105,7 +106,7 @@ us PressureBc::getNEqs(const TaSystem& sys) const { // T = prescribed, if htmodel of duct is not Isentropic // Ts = prescribed, if stempmodel of duct is not Prescribed - const pb::Duct& dpb = getDuct(sys).getDuctPb(); + const pb::Duct& dpb = getDuct().getDuctPb(); bool has_solideq = dpb.stempmodel() != pb::Prescribed; @@ -115,13 +116,8 @@ us PressureBc::getNEqs(const TaSystem& sys) const { VARTRACE(15,neqs); return neqs; } -void PressureBc::show(const TaSystem&,us verbosity_level) const { +void PressureBc::show(us verbosity_level) const { TRACE(15,"PressureBc::show()"); } -void PressureBc::jac(const TaSystem&,Jacobian&, - us dof_start,us eq_start) const { - TRACE(15,"PressureBc::jac()"); - -} ////////////////////////////////////////////////////////////////////// diff --git a/src/ductbc/pressurebc.h b/src/ductbc/pressurebc.h index 838649b..15fe435 100644 --- a/src/ductbc/pressurebc.h +++ b/src/ductbc/pressurebc.h @@ -27,35 +27,33 @@ class PressureBc: public DuctBc { temperature and solid temperature */ pb::DuctSide _side; /**< Duct side at which this b.c. works */ protected: - PressureBc(const PressureBc&); + PressureBc(const TaSystem&,const PressureBc&); public: - PressureBc(const us id, - const TaSystem& sys, + + PressureBc(const TaSystem& sys, + const us id, const pb::DuctBc&); + ~PressureBc(); - PressureBc* copy() const; + + PressureBc* copy(const TaSystem&) const; - vd initialSolution(const TaSystem&) const; + vd initialSolution() const; - virtual void residual(const TaSystem&, - arma::subview_col&& residual + virtual void residualJac(arma::subview_col&& residual ) const; // Return the total number of equations in this segment - virtual us getNEqs(const TaSystem&) const; + virtual us getNEqs() const; // Return the current mass in this segment - virtual d getMass(const TaSystem&) const { return 0;}; + virtual d getMass() const { return 0;}; - virtual void show(const TaSystem&,us verbosity_level) const; + virtual void show(us verbosity_level) const; // 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; - }; diff --git a/src/gui/add_duct_dialog.cpp b/src/gui/add_duct_dialog.cpp index d6eac3d..2d4ce90 100644 --- a/src/gui/add_duct_dialog.cpp +++ b/src/gui/add_duct_dialog.cpp @@ -5,20 +5,22 @@ // Description: // ////////////////////////////////////////////////////////////////////// -#include "ui_add_duct_dialog.h" #include "add_duct_dialog.h" #include "ui_add_duct_dialog.h" +#include "ui_add_duct_dialog.h" + +#include -#include #include #include "tasmet_constants.h" #include "tasmet_enum.h" #include "tasmet_tracer.h" #include "tasmet_assert.h" #include "tasmet_exception.h" +#include "sys/tasystem.h" #include "duct/duct.h" #include "tasmet_qt.h" -#include + DECLARE_ENUM(PreviewShow,CSArea,Porosity,HydraulicRadius,SolidTemperatureFunction) @@ -138,7 +140,9 @@ AddDuctDialog::~AddDuctDialog(){ void AddDuctDialog::accept(){ try { - std::unique_ptr duct (new Duct(0,_duct)); + // If duct can be built without exceptions, everything is OK + TaSystem sys(TaSystem::testSystem()); + Duct duct (sys,0,_duct); } catch(TaSMETError& e) { @@ -216,9 +220,10 @@ void AddDuctDialog::changed(){ PreviewShow pshow = (PreviewShow) _dialog->previewshow->currentIndex(); + TaSystem sys = TaSystem::testSystem(); std::unique_ptr duct; try { - duct = std::unique_ptr(new Duct(0,_duct)); + duct = std::unique_ptr(new Duct(sys,0,_duct)); } catch(TaSMETError& e) { return; diff --git a/src/solver/newton_raphson.cpp b/src/solver/newton_raphson.cpp index d0adcf0..0d1a0a4 100644 --- a/src/solver/newton_raphson.cpp +++ b/src/solver/newton_raphson.cpp @@ -23,43 +23,47 @@ void NewtonRaphson::start_implementation(GradientNonlinearSystem& system, vd guess = system.getSolution(); - vd residual=system.residual(); - #ifdef DEBUG_TASMET_SYSTEM cout << "Initial solution: " << endl; cout << guess << endl; #endif // DEBUG_TASMET_SYSTEM + SolverProgress progress; + SolverAction action; + + ResidualJac resjac; + system.residualJac(resjac); + + const sdmat& jac = resjac.jacobian; + const vd& residual = resjac.residual; + vd dx; + #ifdef DEBUG_TASMET_SYSTEM cout << "Initial residual: " << endl; cout << residual << endl; #endif // DEBUG_TASMET_SYSTEM + assert(jac.n_cols==residual.size()); + assert(jac.n_rows==residual.size()); - SolverProgress progress; - SolverAction action; - while (true) { - - sdmat jac=system.jacobian(); - assert(jac.n_cols==residual.size()); - assert(jac.n_rows==residual.size()); + #ifdef DEBUG_TASMET_SYSTEM + cout << "Residual: "; + cout << residual << endl; + #endif // DEBUG_TASMET_SYSTEM + + dx = -1*_dampfac*arma::spsolve(jac,residual,"superlu"); - vd dx = -1*_dampfac*arma::spsolve(jac,residual,"superlu"); + progress.rel_err = norm(dx); guess += dx; system.updateSolution(guess); - residual = system.residual(); - #ifdef DEBUG_TASMET_SYSTEM - cout << "Residual: "; - cout << residual << endl; - #endif // DEBUG_TASMET_SYSTEM + // Obtain a new residual and Jacobian matrix + system.residualJac(resjac); - - progress.rel_err = norm(dx); progress.fun_err = norm(residual); progress.iteration++; diff --git a/src/solver/system.h b/src/solver/system.h index 6d128e9..fe7ecbe 100644 --- a/src/solver/system.h +++ b/src/solver/system.h @@ -62,13 +62,28 @@ public: }; -class GradientNonlinearSystem: public NoGradientNonlinearSystem { +struct ResidualJac { + vd residual; // Store residual here + sdmat jacobian; // Store Jacobian here +}; + +class GradientNonlinearSystem { public: - - virtual sdmat jacobian() const=0; + /** + * Computes both the residual as well as the Jacobian, and stores + * them in the input parameter + * + */ + virtual void residualJac(ResidualJac&) const=0; + + virtual GradientNonlinearSystem* copy() const=0; + + // Obtain an initial guess of the solution + virtual vd getSolution() const=0; + + virtual void updateSolution(const vd& new_guess)=0; - GradientNonlinearSystem* copy() const=0; }; diff --git a/src/sys/segment.h b/src/sys/segment.h index 0280d09..441545b 100644 --- a/src/sys/segment.h +++ b/src/sys/segment.h @@ -24,21 +24,37 @@ class TaSystem; class Segment{ protected: + + const TaSystem& sys; + // ID us _id; // Name std::string _name; - - Segment(const us id,const std::string& name): _id(id),_name(name) {} - Segment(const Segment& o): Segment(o._id,o._name){} + Segment(const TaSystem& sys, + const us id, + const std::string& name): sys(sys),_id(id),_name(name) {} + Segment(const TaSystem& sys,const Segment& o): + Segment(sys,o._id,o._name){} Segment& operator=(const Segment&)=delete; public: virtual ~Segment(){} - virtual Segment* copy() const = 0; + /** + * Returns a copy of the current segment, coupled to the TaSystem + * as given function parameter. + * + * @param sys: The new TaSystem to which this segment should be + * coupled. + * @return the Segment copy + */ + virtual Segment* copy(const TaSystem&) const = 0; - virtual vd initialSolution(const TaSystem&) const = 0; + virtual vd initialSolution() const = 0; + + virtual void residualJac(arma::subview_col&& residual // Here we store the residual + ) const=0; // Get and set name const std::string& getName() const{return _name;} // This one is just the name @@ -49,30 +65,22 @@ public: // 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&, - arma::subview_col&& residual // Here we store the residual - ) const=0; - // Return the total number of equations in this segment - virtual us getNEqs(const TaSystem&) const { return 0;} + virtual us getNEqs() const { return 0;} // Return the total number of DOFS in this segment - virtual us getNDofs(const TaSystem&) const { return 0;} + virtual us getNDofs() 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 d getMass() const = 0; + virtual void dmtotdx(vd& dmtotdx,us dof_start) const {} - virtual void show(const TaSystem&,us verbosity_level) const=0; + virtual void show(us verbosity_level) const=0; // 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=0; - }; diff --git a/src/sys/tasystem.cpp b/src/sys/tasystem.cpp index f3d96a1..0922e18 100644 --- a/src/sys/tasystem.cpp +++ b/src/sys/tasystem.cpp @@ -34,7 +34,7 @@ TaSystem::TaSystem(const pb::System& sys): // d.first: id // d.second: duct description try { - _segs[d.first] = new Duct(d.first,d.second); + _segs[d.first] = new Duct(*this,d.first,d.second); if(!_segs[d.first]) throw TaSMETBadAlloc(); } catch(TaSMETError e) { @@ -69,6 +69,9 @@ TaSystem::TaSystem(const pb::System& sys): vd solution = vd(total_dofs); us i=0; + + + const Segment* seg; for(auto& seg_: _segs) { seg = seg_.second; @@ -76,11 +79,11 @@ TaSystem::TaSystem(const pb::System& sys): if(ndofs(i)>0) { if(i==0) { solution.subvec(0,ndofs(0)-1) = - seg->initialSolution(*this); + seg->initialSolution(); } else { solution.subvec(dofend(i-1),dofend(i)-1) = - seg->initialSolution(*this); + seg->initialSolution(); } i++; } @@ -111,7 +114,7 @@ TaSystem::TaSystem(const TaSystem& o): _segs = o._segs; for(auto& seg: _segs){ - seg.second = seg.second->copy(); + seg.second = seg.second->copy(*this); if(!seg.second) throw TaSMETBadAlloc(); } } @@ -150,7 +153,7 @@ int TaSystem::getArbitrateMassEq(const vus& neqs) const { } return arbitrateMassEq; } -vd TaSystem::residual() const { +void TaSystem::residualJac(ResidualJac& resjac) const { TRACE(15,"TaSystem::residual()"); vus neqs = getNEqs(); @@ -190,15 +193,15 @@ vd TaSystem::residual() const { if(i==0) { // Put the residual of the segment in the over-all residual - seg->residual(*this,residual.subvec(0,eqsend(0)-1)); + seg->residualJac(residual.subvec(0,eqsend(0)-1)); } else { - seg->residual(*this,residual.subvec(eqsend(i-1),eqsend(i)-1)); + seg->residualJac(residual.subvec(eqsend(i-1),eqsend(i)-1)); } // Count the mass, add it if(arbitrateMassEq!=-1) { - mass += seg->getMass(*this); + mass += seg->getMass(); } i++; @@ -215,7 +218,7 @@ vd TaSystem::residual() const { residual(arbitrateMassEq)=mass - _mass; } - return residual; + resjac.residual = residual; } vd TaSystem::getSolution() const { @@ -240,7 +243,7 @@ vus TaSystem::getNDofs() const { vus Ndofs(_segs.size()); us i=0; for (auto seg : _segs) { - Ndofs(i)=seg.second->getNDofs(*this); + Ndofs(i)=seg.second->getNDofs(); i++; } return Ndofs; @@ -250,7 +253,7 @@ vus TaSystem::getNEqs() const { vus Neqs(_segs.size()); us i=0; for (auto seg :_segs) { - Neqs(i)=seg.second->getNEqs(*this); + Neqs(i)=seg.second->getNEqs(); VARTRACE(15,Neqs(i)); i++; } @@ -266,78 +269,73 @@ void TaSystem::show(us detailnr){ if(detailnr>0){ for(auto seg:_segs){ cout << "Showing segment with ID " << seg.first << "\n"; - seg.second->show(*this,detailnr); + seg.second->show(detailnr); } } // detailnr>0 } -TripletList TaSystem::jacTriplets() const { +// 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); +// 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); +// 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++; - } +// 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"); - } +// if(dofnr!=eqnr){ +// throw TaSMETError("System of equations is over/underdetermined"); +// } - // Convert to tripletlist - TripletList jac=j; +// // Convert to tripletlist +// TripletList jac=j; - assert(arbitrateMassEq< (int) total_neqs); +// 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); +// // 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(){ +// dmat TaSystem::showJac(){ - TRACE(15,"TaSystem::showJac()"); - - return dmat(jacobian()); -} +// TRACE(15,"TaSystem::showJac()"); +// return dmat(jacobian()); +// } TaSystem::~TaSystem() { TRACE(25,"~TaSystem()"); cleanup(); diff --git a/src/sys/tasystem.h b/src/sys/tasystem.h index 204b7c5..fe06374 100644 --- a/src/sys/tasystem.h +++ b/src/sys/tasystem.h @@ -28,14 +28,24 @@ protected: std::unique_ptr _gas; - vd _solution; + vd _solution; /**< This column vector stores the + current solution. */ + + std::vector _dofs; /**< This vector of column vectors + stores columns of doubles that do + not own the memory they point + to. Instead, they use the memory of + _solution.*/ + vus _startdof; // Store the start DOFS for each segment TaSystem& operator=(const TaSystem& other)=delete; - TaSystem(const TaSystem& o); public: + TaSystem(const TaSystem& o); TaSystem(const pb::System&); - + static TaSystem testSystem() { + return TaSystem(pb::System::default_instance()); + } const Gas& gas() const {return *_gas;} // Set and get the mass in the system. If the mass is not set @@ -52,7 +62,7 @@ public: dmat showJac(); virtual void show(us detailnr=0); - vd residual() const; + void residualJac(ResidualJac&) const; vd getSolution() const; @@ -61,10 +71,6 @@ public: virtual void updateSolution(const vd& sol) {_solution = sol; } // 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); @@ -81,7 +87,6 @@ public: protected: virtual int getArbitrateMassEq(const vus& neqs) const; - virtual TripletList jacTriplets() const; private: void cleanup(); }; // class System