anne/tasmet
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

In progress. Ready to add Jac to residualJac

Browse Source
This commit is contained in:
Anne de Jong 2017-01-22 10:34:46 +01:00
parent b7bd1df4aa
commit a71159cf7f
14 changed files with 283 additions and 258 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
src/duct/duct.cpp
View File

@ -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)
@ -48,10 +50,9 @@ Duct::Duct(const Duct& other):
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<d> && residual) const {
void Duct::residualJac(arma::subview_col<d> && residual) const {
TRACE(15,"Duct::residual()");
@ -86,10 +87,10 @@ void Duct::residual(const TaSystem& sys,arma::subview_col<d> && 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<d> && 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<d> && 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<d> 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<ngp();i++){
res(i) = getvart(sys,varnr,i)(t);
res(i) = getvart(varnr,i)(t);
}
return res;
}
vd Duct::initialSolution(const TaSystem& sys) const {
vd Duct::initialSolution() const {
TRACE(15,"Duct::initialSolution()");
vd initsol(getNDofs(sys));
vd initsol(getNDofs());
VARTRACE(15,initsol.size());
us vars_per_gp = 4;
@ -228,7 +229,7 @@ vd Duct::initialSolution(const TaSystem& sys) const {
}
us Duct::getNEqs(const TaSystem& sys) const {
us Duct::getNEqs() const {
TRACE(15,"Duct::getNEqs()");
us Ns = sys.Ns();
@ -253,7 +254,7 @@ us Duct::getNEqs(const TaSystem& sys) const {
VARTRACE(15,neqs);
return neqs;
}
us Duct::getNDofs(const TaSystem& sys) const {
us Duct::getNDofs() const {
TRACE(15,"Duct::getNDofs()");
us Ns = sys.Ns();
@ -264,14 +265,11 @@ us Duct::getNDofs(const TaSystem& sys) const {
return Ns*nvars_per_gp*ngp();
}
d Duct::getMass(const TaSystem& sys) const {
d Duct::getMass() const {
return 0;
}
void Duct::jac(const TaSystem&,Jacobian&,us dof_start,us eq_start) const {
}
void Duct::show(const TaSystem&,us verbosity_level) const {
void Duct::show(us verbosity_level) const {
}

64
src/duct/duct.h
View File

@ -22,16 +22,14 @@ class Duct : public Segment, public Geom {
// Drag* _drag = nullptr;
// Heat* _heat = nullptr;
std::vector<Equation*> _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<d>&& residual) const;
virtual void residualJac(arma::subview_col<d>&& 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

30
src/ductbc/adiabaticwall.cpp
View File

@ -12,19 +12,19 @@
#include "tasystem.h"
#include "duct.h"
AdiabaticWall::AdiabaticWall(const us id,
const TaSystem& sys,
AdiabaticWall::AdiabaticWall(const TaSystem& sys,
const us id,
const pb::DuctBc& dbc):
DuctBc(id,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<d>&& residual
void AdiabaticWall::residualJac(arma::subview_col<d>&& 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()");
}
//////////////////////////////////////////////////////////////////////

21
src/ductbc/adiabaticwall.h
View File

@ -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,
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<d>&& residual
virtual void residualJac(arma::subview_col<d>&& 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;

12
src/ductbc/ductbc.cpp
View File

@ -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<const Duct&>(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;

11
src/ductbc/ductbc.h
View File

@ -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 Duct& getDuct(const TaSystem&) const;
const Duct& getDuct() const;
};

46
src/ductbc/pressurebc.cpp
View File

@ -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<d>&& residual
void PressureBc::residualJac(arma::subview_col<d>&& 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()");
}
//////////////////////////////////////////////////////////////////////

26
src/ductbc/pressurebc.h
View File

@ -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;
vd initialSolution(const TaSystem&) const;
PressureBc* copy(const TaSystem&) const;
virtual void residual(const TaSystem&,
arma::subview_col<d>&& residual
vd initialSolution() const;
virtual void residualJac(arma::subview_col<d>&& 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;
};

15
src/gui/add_duct_dialog.cpp
View File

@ -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 <QMessageBox>
#include <QSignalBlocker>
#include <qcustomplot.h>
#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 <QMessageBox>
DECLARE_ENUM(PreviewShow,CSArea,Porosity,HydraulicRadius,SolidTemperatureFunction)
@ -138,7 +140,9 @@ AddDuctDialog::~AddDuctDialog(){
void AddDuctDialog::accept(){
try {
std::unique_ptr<Duct> 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> duct;
try {
duct = std::unique_ptr<Duct>(new Duct(0,_duct));
duct = std::unique_ptr<Duct>(new Duct(sys,0,_duct));
}
catch(TaSMETError& e) {
return;

36
src/solver/newton_raphson.cpp
View File

@ -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
SolverProgress progress;
SolverAction action;
while (true) {
sdmat jac=system.jacobian();
assert(jac.n_cols==residual.size());
assert(jac.n_rows==residual.size());
vd dx = -1*_dampfac*arma::spsolve(jac,residual,"superlu");
while (true) {
guess += dx;
system.updateSolution(guess);
residual = system.residual();
#ifdef DEBUG_TASMET_SYSTEM
cout << "Residual: ";
cout << residual << endl;
#endif // DEBUG_TASMET_SYSTEM
dx = -1*_dampfac*arma::spsolve(jac,residual,"superlu");
progress.rel_err = norm(dx);
guess += dx;
system.updateSolution(guess);
// Obtain a new residual and Jacobian matrix
system.residualJac(resjac);
progress.fun_err = norm(residual);
progress.iteration++;

21
src/solver/system.h
View File

@ -62,13 +62,28 @@ public:
};
class GradientNonlinearSystem: public NoGradientNonlinearSystem<vd> {
struct ResidualJac {
vd residual; // Store residual here
sdmat jacobian; // Store Jacobian here
};
class GradientNonlinearSystem {
public:
/**
* Computes both the residual as well as the Jacobian, and stores
* them in the input parameter
*
*/
virtual void residualJac(ResidualJac&) const=0;
virtual sdmat jacobian() 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;
};

44
src/sys/segment.h
View File

@ -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<d>&& 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<d>&& 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;
};

126
src/sys/tasystem.cpp
View File

@ -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;k<dmtotdxsize;k++)
if(dmtotdx(k)!=0){
// TRACE(20,"k: " << k);
// TRACE(20,"dmtotdx:"<< dmtotdx(k));
jac.addTriplet(Triplet(arbitrateMassEq,k,dmtotdx(k)));
}
}
// 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
}
// return jac;
// }
// void TaSystem::resetHarmonics(){
// for(auto seg: _segs) {
// seg.second->resetHarmonics();
// }
// }
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();

23
src/sys/tasystem.h
View File

@ -28,14 +28,24 @@ protected:
std::unique_ptr<Gas> _gas;
vd _solution;
vd _solution; /**< This column vector stores the
current solution. */
std::vector<vd> _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