diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9e601a8..36d37f1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -33,6 +33,7 @@ add_library(tasmet_src duct/grid.cpp duct/geom.cpp + duct/duct.cpp funcs/bessel.cpp funcs/cbessj.cpp diff --git a/src/duct/drag.cpp b/src/duct/drag.cpp new file mode 100644 index 0000000..6b67eed --- /dev/null +++ b/src/duct/drag.cpp @@ -0,0 +1,61 @@ +#include "drag.h" +#include "bccell.h" +#include "jacrow.h" + +namespace duct{ + namespace drag { + using tasystem::JacRow; + using std::make_tuple; + using std::tuple; + + vd DragResistance::drag(const Cell& v) const {return vd(v.gc->Ns(),fillwith::zeros);} + JacRow DragResistance::dDrag(const Cell& v) const {return JacRow(-1);} + // dmat DragResistance::drhoi(const Cell& v) const {return v.zero;} + // dmat DragResistance::dpi(const Cell& v) const {return v.zero;} + + d DragResistance::rho0(const Cell& v) { + TRACE(15,"DragResistance::rho0()"); + d T0=v.T()(0); + d p0=v.p()(0)+v.gc->p0(); + return v.gc->gas().rho(T0,p0); + } + d DragResistance::mu0(const Cell& v) { + TRACE(15,"DragResistance::mu0()"); + d T0=v.T()(0); + return v.gc->gas().mu(T0); + } + + vd DragResistance::shearWaveNumber(const Cell& v) const { + TRACE(15,"DragResistance::shearWaveNumber(const Cell& v)"); + const us Nf=v.gc->Nf(); + + const d omg=v.gc->getomg(); + const vd omgvec=linspace(0,Nf*omg,Nf+1); + d rh=v.rh(); + + return rh*sqrt((rho0(v)/mu0(v))*omgvec); + + } + } // namespace drag +} // namespace duct + + + + + + + + + + + + + + + + + + + + + diff --git a/src/duct/drag.h b/src/duct/drag.h new file mode 100644 index 0000000..15f0c0e --- /dev/null +++ b/src/duct/drag.h @@ -0,0 +1,43 @@ +#pragma once +#ifndef _DRAG_H_ +#define _DRAG_H_ +#include "tasmet_types.h" + +class DragResistance{ +public: + virtual ~DragResistance(){} + virtual vd drag(us gp) const; + virtual tasystem::JacRow dDrag(const Cell& v) const; + // virtual dmat drhoi(const Cell& v) const; + // virtual dmat dpi(const Cell& v) const; + + // Provide a vector of shear wave numbers for each frequency. The + // length of this vector is Nf+1, and obviously, the first // of this vector is always zero. + vd shearWaveNumber(const Cell& v) const; + + static d rho0(const Cell& v); + static d mu0(const Cell& v); +}; + +#endif /* _DRAG_H_ */ + + + + + + + + + + + + + + + + + + + + + diff --git a/src/duct/duct.cpp b/src/duct/duct.cpp new file mode 100644 index 0000000..ad64125 --- /dev/null +++ b/src/duct/duct.cpp @@ -0,0 +1,48 @@ +// duct.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#include "duct.h" +#include "tasystem.h" + +Duct::Duct(const pb::Duct& duct): + Segment(duct.name()), + Geom(duct) +{ + + + + +} +void Duct::updateSolution(const TaSystem&,const vd&) const { + +} +vd Duct::getSolution(const TaSystem& sys,vd& sol,const us insert_start) { + + + +} + + +Duct::Duct(const Duct& other): + Segment(other), + Geom(other) +{ + + +} +us Duct::getNEqs(const TaSystem& sys) const { + + us Ns = sys.Ns(); + us number_eqs = _eqs.size(); + + return Ns*number_eqs*ngp(); +} + +Duct* Duct::copy() const { + return new Duct(*this); +} +////////////////////////////////////////////////////////////////////// diff --git a/src/duct/duct.h b/src/duct/duct.h index 3e3d3dd..4c85bc1 100644 --- a/src/duct/duct.h +++ b/src/duct/duct.h @@ -10,23 +10,55 @@ #define DUCT_H #include "segment.h" #include "duct.pb.h" +#include "geom.h" class Equation; class Drag; class Heat; +class TaSystem; + +class Duct : public Segment, public Geom { + + Duct(const Duct& other); + + // Drag* _drag = nullptr; + // Heat* _heat = nullptr; + std::vector _eqs; + + std::vector _rho; + std::vector _u; + std::vector _T; + std::vector _p; + std::vector _Ts; -class Duct { - Duct(const string& name,const Geom& geom); public: - + Duct(const pb::Duct&); + virtual Duct* copy() const; const Geom& geom() const; - - - // Parsing a protobuf to generate a NEW Duct - static Duct* parseProto(const pb::Duct&); - - + // Solving + virtual void residual(const TaSystem&,vd&,const us insertion_start) const; + + virtual void updateSolution(const TaSystem&,const vd&); + virtual getSolution(const TaSystem&,vd& sol,const us insertion_start) const; + + // Return the total number of equations in this segment + virtual us getNEqs(const TaSystem&) const; + // Return the total number of DOFS in this segment + virtual us getNDofs(const TaSystem&) 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 void show(const TaSystem&,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; }; diff --git a/src/duct/eq.cpp b/src/duct/eq.cpp new file mode 100644 index 0000000..065b2e2 --- /dev/null +++ b/src/duct/eq.cpp @@ -0,0 +1,18 @@ +#include "ductequation.h" +#include "cell.h" + + +namespace duct{ + vd Equation::getp0t() const { + TRACE(0,"Equation::getp0t()"); + return v.gc->p0()*vd(v.gc->Ns(),fillwith::ones); + } + dmat Equation::eye() const { + TRACE(15,"Equation::eye()"); + return arma::eye(v.gc->Ns(),v.gc->Ns());} + dmat Equation::eye(const Cell& v) { + TRACE(15,"Equation::eye()"); + return arma::eye(v.gc->Ns(),v.gc->Ns());} + vd Equation::zeros() const {return vd(v.gc->Ns(),fillwith::zeros);} + +} // namespace duct diff --git a/src/duct/eq.h b/src/duct/eq.h new file mode 100644 index 0000000..7c416ea --- /dev/null +++ b/src/duct/eq.h @@ -0,0 +1,46 @@ +// eq.h +// +// Author: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef EQ_H +#define EQ_H + +#include "constants.h" + + +class JacRow; +class JacCol; + +class Duct; +class WeightFactors; + +class Equation{ +protected: + us dofnr; +public: + Equation(const Cell& v):v(v){TRACE(15,"Equation(v)");} + virtual enum EqType getType() const=0; + + // Return an eye of the right size: + dmat eye() const; + static dmat eye(const Cell&); + vd zeros() const; + virtual tasystem::JacRow jac() const=0; // Returns the local Jacobian of this equation + virtual vd error() const=0; + virtual void show() const=0; + virtual void domg(vd&) const {/* Placeholder */} + + vd getp0t() const; + vd getT0t() const; + virtual ~Equation(){} +private: + vd nu(const Cell&) const; // Function of d^2p/dx^2 +}; // class Equation + + +#endif // EQ_H +////////////////////////////////////////////////////////////////////// diff --git a/src/duct/laminardrag.cpp b/src/duct/laminardrag.cpp new file mode 100644 index 0000000..5c938c0 --- /dev/null +++ b/src/duct/laminardrag.cpp @@ -0,0 +1,105 @@ +#include "duct.h" +#include "bccell.h" +#include "laminardrag.h" +#include "geom.h" +#include "math_constants.h" +#include "jacrow.h" + +namespace duct{ + namespace drag { + using math_common::sq2; + using rottfuncs::RottFuncs; + using tasystem::JacRow; + using tasystem::JacCol; + + // Resistance force for laminar flow for the zero-frequency. + d zerodrag_vert(const Cell& v){ + TRACE(5,"zerodrag_vert"); + return 3*DragResistance::mu0(v)/(DragResistance::rho0(v)*pow(v.rhl,2)); + } + d zerodrag_circ(const Cell& v){ + TRACE(5,"zerodrag_circ"); + return 2*DragResistance::mu0(v)/(DragResistance::rho0(v)*pow(v.rhl,2)); + } + d zerodrag_inviscid(const Cell& v){ + TRACE(5,"zerodrag_inviscid"); + return 0; + } + + class ZeroFreqDragCoef { + d (*zerodrag_funptr)(const Cell&); + public: + ZeroFreqDragCoef(const Duct& t) { + TRACE(0,"ZeroFreqDragCoef::ZeroFreqDragCoef()"); + if(t.geom().shape().compare("vert")==0) + zerodrag_funptr=&zerodrag_vert; + else if(t.geom().shape().compare("circ")==0){ + TRACE(20,"Circular pore chosen"); + zerodrag_funptr=&zerodrag_circ; + } + else if(t.geom().shape().compare("inviscid")==0) + zerodrag_funptr=&zerodrag_inviscid; + else { + WARN("Warning: duct.geom.shape unknown for ZeroFreqDrag. Aborting..."); + abort(); + } + } + // We implement a cheap variant of polymorphism + d operator()(const Cell& v) const { + return (*zerodrag_funptr)(v); + } + }; + + + LaminarDragResistance::LaminarDragResistance(const Duct& t) + { + TRACE(10,"LaminarDragResistanc::LaminarDragResistance()"); + TRACE(11,"Entering redefinition of Rottfuncs"); + if(t.geom().isBlApprox()) + rf=RottFuncs("blapprox"); + else + rf=RottFuncs(t.geom().shape()); // Reinitialize thermoviscous functions with right shape + TRACE(11,"Exiting redefinition of Rottfuncs"); + zfd=new ZeroFreqDragCoef(t); + } + LaminarDragResistance::~LaminarDragResistance(){ + delete zfd; + } + vd LaminarDragResistance::drag(const Cell& v) const { + TRACE(10,"LaminarDragResistance::drag(v)"); + vd drag=dm(v)*v.ml()(); + return drag; // No momentum scale here, since this is already done in dUi!!!! + } + JacRow LaminarDragResistance::dDrag(const Cell& v) const { + TRACE(15,"LaminarDragResistance::dDrag()"); + return JacRow(-1,JacCol(v.ml(),dm(v))); + } + dmat LaminarDragResistance::dm(const Cell& v) const { // Derivative of drag resistance to velocity + TRACE(10,"LaminarDragResistance::dUi()"); + vc CResistance=ComplexResistancecoef(v); + tasystem::var resistance(*v.gc); + resistance.setadata(CResistance); + return resistance.freqMultiplyMat(); + } + vc LaminarDragResistance::ComplexResistancecoef(const Cell& v) const { + TRACE(0,"LaminarDragResistance::ComplexResistancecoef()"); + const us& Nf=v.gc->Nf(); + const us& i=v.geti(); + + const d& rh=v.rhl; + + vc rescoef(Nf+1); + rescoef(0)=(*zfd)(v); // Zero frequency drag divided by zero-frequency velocity + if(Nf>0){ + // Divided by sq2, see Eq. 5.23 of my thesis + const d omg=v.gc->getomg(); + const vd omgvec=omg*linspace(1,Nf,Nf); + vd rh_over_deltanu=(shearWaveNumber(v).subvec(1,Nf))/sq2; + vc fnu=rf.fx(rh_over_deltanu); // Viscous rott function + rescoef.subvec(1,Nf)=I*(omgvec%(fnu/(1.0-fnu))); + } + return rescoef; + } + } // namespace drag +} // namespace duct + diff --git a/src/duct/laminardrag.h b/src/duct/laminardrag.h new file mode 100644 index 0000000..d74f7fb --- /dev/null +++ b/src/duct/laminardrag.h @@ -0,0 +1,36 @@ +#include "drag.h" +#include "rottfuncs.h" + +namespace duct{ + SPOILNAMESPACE + namespace drag { + + // The Drag coefficient for "frequency" zero. + class ZeroFreqDragCoef; + + // Laminar drag resistance + class LaminarDragResistance:public DragResistance + { + ZeroFreqDragCoef* zfd; + rottfuncs::RottFuncs rf; + public: + LaminarDragResistance(const Duct& t); + LaminarDragResistance(const LaminarDragResistance&)=delete; + LaminarDragResistance& operator=(const LaminarDragResistance&)=delete; + ~LaminarDragResistance(); + + // Overloaded virtuals + vd drag(const Cell& cell) const; + tasystem::JacRow dDrag(const Cell&) const; + + private: + dmat dm(const Cell&) const; // Derivative of drag resistance + // to volume flow + + // Returns a complex vector of size Ns with drag resistance + // coefficients for every nonzero frequency (1..Nf) + vc ComplexResistancecoef(const Cell&) const; + + }; + } // namespace drag +} // namespace duct diff --git a/src/sys/segment.h b/src/sys/segment.h index cb95920..d2e8128 100644 --- a/src/sys/segment.h +++ b/src/sys/segment.h @@ -46,8 +46,8 @@ public: 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; + virtual void updateSolution(const TaSystem&,const vd&) = 0; + virtual 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;} @@ -61,7 +61,7 @@ public: virtual void show(const TaSystem&,us verbosity_level) const=0; // Reset amplitude data in higher harmonics - void resetHarmonics(); + // virtual void resetHarmonics() = 0; // Fill Jacobian with values from the equations in this // segment/connector. diff --git a/src/sys/tasystem.cpp b/src/sys/tasystem.cpp index 2156a69..cbd341a 100644 --- a/src/sys/tasystem.cpp +++ b/src/sys/tasystem.cpp @@ -14,13 +14,13 @@ #include "tasmet_constants.h" TaSystem::TaSystem(const GlobalConf& gc,const Gas& g): - _gc(new GlobalConf(gc)), + GlobalConf(gc)), _gas(g.copy()) { TRACE(14,"TaSystem::TaSystem(gc,gastype)"); } TaSystem::TaSystem(const TaSystem& o): - _gc(o._gc), // Share a ptr to the Global conf + GlobalConf(o), // Share a ptr to the Global conf _gas(o._gas->copy()) { TRACE(25,"TaSystem::TaSystem(TaSystem&) copy"); @@ -206,7 +206,7 @@ void TaSystem::show(us detailnr){ cout << "########################## Showing TaSystem...\n"; cout << "Showing Global configuration...\n"; - _gc->show(); + GlobalConf::show(); if(detailnr>0){ for(auto seg:_segs){ @@ -272,11 +272,11 @@ sdmat TaSystem::jacobian() const { TRACE(14,"TaSystem::Jac()"); return jacTriplets(); // Implicitly converts to sdmat } -void TaSystem::resetHarmonics(){ - for(auto seg: _segs) { - seg.second->resetHarmonics(); - } -} +// void TaSystem::resetHarmonics(){ +// for(auto seg: _segs) { +// seg.second->resetHarmonics(); +// } +// } dmat TaSystem::showJac(){ TRACE(15,"TaSystem::showJac()"); diff --git a/src/sys/tasystem.h b/src/sys/tasystem.h index f20a737..fa0533b 100644 --- a/src/sys/tasystem.h +++ b/src/sys/tasystem.h @@ -16,15 +16,13 @@ #include // Inherit all global configuration members -class TaSystem : public GradientNonlinearSystem { +class TaSystem : public GradientNonlinearSystem, public GlobalConf { protected: d _mass = -1; std::map _segs; - gc_ptr _gc; - Gas* _gas = nullptr; TaSystem& operator=(const TaSystem& other)=delete; @@ -68,7 +66,7 @@ public: void updateNf(us); // Reset amplitude data in higher harmonics - void resetHarmonics(); + // void resetHarmonics(); // void delseg(us n); // Not yet implemented. Delete a segment // from the system (we have to determine how elaborated the API