From c81e17ceccceaf0fc513f885b3a0feb25f093e5b Mon Sep 17 00:00:00 2001 From: Anne de Jong Date: Thu, 26 Jan 2017 15:03:34 +0100 Subject: [PATCH] Working! New way of obtaining Jacobian succesfully implemented. Newton-Raphson solver working --- src/CMakeLists.txt | 3 +- src/duct/duct.cpp | 171 ++++++++++++++-------- src/duct/duct.h | 12 +- src/ductbc/adiabaticwall.cpp | 19 +-- src/ductbc/adiabaticwall.h | 11 +- src/ductbc/ductbc.cpp | 8 +- src/ductbc/ductbc.h | 2 + src/ductbc/pressurebc.cpp | 21 ++- src/ductbc/pressurebc.h | 6 +- src/material/gas.h | 8 +- src/solver/newton_raphson.cpp | 24 ++-- src/sys/jaccol.cpp | 13 -- src/sys/jaccol.h | 50 ------- src/sys/jacobian.cpp | 8 -- src/sys/jacobian.h | 54 +++++-- src/sys/jacrow.cpp | 40 ------ src/sys/jacrow.h | 55 -------- src/sys/segment.h | 8 +- src/sys/tasystem.cpp | 259 +++++++++++++++++++++------------- src/sys/tasystem.h | 43 ++++-- 20 files changed, 416 insertions(+), 399 deletions(-) delete mode 100644 src/sys/jaccol.cpp delete mode 100644 src/sys/jaccol.h delete mode 100644 src/sys/jacrow.cpp delete mode 100644 src/sys/jacrow.h diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 632f63a..7e9e488 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -65,9 +65,8 @@ add_library(tasmet_src # sys/enginesystem.cpp sys/globalconf.cpp sys/tasmet_variable.cpp - sys/jaccol.cpp + sys/jacobian.cpp - sys/jacrow.cpp sys/tasystem.cpp sys/triplets.cpp ) diff --git a/src/duct/duct.cpp b/src/duct/duct.cpp index ee4884b..e02ceee 100644 --- a/src/duct/duct.cpp +++ b/src/duct/duct.cpp @@ -60,12 +60,11 @@ Duct::~Duct() { // delete eq; // } } -void Duct::residualJac(arma::subview_col && residual) const { +void Duct::residualJac(SegPositionMapper& residual, + SegJacRows& jacrows) const { - TRACE(15,"Duct::residual()"); + TRACE(15,"Duct::residualJac()"); - const arma::subview_col sol = sys.getSolution(_id); - vd rho,u,T,p,Ts; // Solution at this gp vd rhop,up,Tp,pp,Tsp; // Solution at next gp @@ -78,28 +77,26 @@ void Duct::residualJac(arma::subview_col && residual) const { us number_eqs = 4; number_eqs += (has_solideq) ? 1 : 0; + // rho,u,T,p + // Ts maybe + us nvars_per_gp = (_ductpb.stempmodel() + == pb::HeatBalance ? 5 : 4); + VARTRACE(15,number_eqs); us Ns = sys.Ns(); - us eq_offset = 0; // Equation offset for current gp - us res_offset = 0; // Residual offset for current gp - us res_offsetp = 0; // Residual offset for next gp - us gp_jump = number_eqs * Ns; // The jump per gp - rhop = getvart(constants::rho,0); up = getvart(constants::u,0); Tp = getvart(constants::T,0); - pp = getvart(constants::p,0); + pp = getvart(constants::p,0); const Gas& gas = sys.gas(); + us i=0; + for(us gp=0;gp && residual) const { pp = getvart(constants::p,gp+1); cont = ((rhop%up)-(rho%u))/dx; + *residual.at(i) = cont; + + jacrows.at(i)[{_id,gp*nvars_per_gp+constants::u}] = + -diagmat(rho)/dx; + jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] = + -diagmat(u)/dx; + jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::rho}] = + diagmat(up)/dx; + jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::u}] = + diagmat(rhop)/dx; + TRACE(15,"SFSG"); + i++; mom = (rhop%up%up - rho%u%u + pp - p)/dx; + jacrows.at(i)[{_id,gp*nvars_per_gp+constants::u}] = + -diagmat(rho%u)/dx; + jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] = + -diagmat(u%u)/dx; + jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::rho}] = + diagmat(u%u)/dx; + jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::u}] = + diagmat(up%rhop)/dx; + jacrows.at(i)[{_id,(gp)*nvars_per_gp+constants::p}] = + -eye(Ns,Ns)/dx; + jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::p}] = + eye(Ns,Ns)/dx; + *residual.at(i) = mom; + TRACE(15,"SFSG"); + i++; switch (_ductpb.htmodel()) { case pb::Isentropic: { @@ -125,33 +149,53 @@ void Duct::residualJac(arma::subview_col && residual) const { en = p/p0 - pow(rho/rho0,gamma0); + *residual.at(i) = en; + + jacrows.at(i)[{_id,gp*nvars_per_gp+constants::p}] = + eye(Ns,Ns)/p0; + jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] = + -(gamma0/rho0)*diagmat(pow(rho/rho0,gamma0-1)); + TRACE(15,"SFSG"); } break; default: tasmet_assert(false,"Not implemented htmodel"); } + i++; + TRACE(15,"SFSG"); + // Equation of state st = gas.rho(T,p) - rho; - - residual.subvec(eq_offset+0*Ns,eq_offset+1*Ns-1) = cont; - residual.subvec(eq_offset+1*Ns,eq_offset+2*Ns-1) = mom; - residual.subvec(eq_offset+2*Ns,eq_offset+3*Ns-1) = en; - residual.subvec(eq_offset+3*Ns,eq_offset+4*Ns-1) = st; - + TRACE(15,"SFSG"); + jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] = + -eye(Ns,Ns); + jacrows.at(i)[{_id,gp*nvars_per_gp+constants::p}] = + diagmat(gas.drhodp(T,p)); + jacrows.at(i)[{_id,gp*nvars_per_gp+constants::T}] = + diagmat(gas.drhodT(T,p)); + + *residual.at(i) = st; + TRACE(15,"SFSG"); + i++; } - eq_offset += number_eqs*Ns; - // Equation of state for the last node st = gas.rho(Tp,pp) - rhop; - residual.subvec(eq_offset,eq_offset+Ns-1) = st; + *residual.at(i) = st; + + jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::rho}] = + -eye(Ns,Ns); + jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::p}] = + diagmat(gas.drhodp(T,p)); + jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::T}] = + diagmat(gas.drhodT(T,p)); + + i++; // Two more equations for the last grid point in case // the heat transfer model is not a transport equation. if(_ductpb.htmodel() == pb::Isentropic) { - eq_offset += Ns; - d T0 = gas.T0(); d p0 = gas.p0(); d rho0 = gas.rho0(); @@ -159,24 +203,43 @@ void Duct::residualJac(arma::subview_col && residual) const { en = p/p0 - pow(rho/rho0,gamma0); - residual.subvec(eq_offset,eq_offset+Ns-1) = en; + jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::p}] = + eye(Ns,Ns)/p0; + jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::rho}] = + -(gamma0/rho0)*diagmat(pow(rho/rho0,gamma0-1)); + + *residual.at(i) = en; i++; } } +PosId Duct::getDof(int varnr,int gp) const { + // Wraparound + if(gp<0) gp+=ngp(); + + // rho,u,T,p + // Ts maybe + us nvars_per_gp = (_ductpb.stempmodel() + == pb::HeatBalance ? 5 : 4); + + + return {_id,nvars_per_gp*gp+varnr}; +} + vd Duct::getvart(int varnr,int gp) const { - TRACE(15,"Duct::getvart()"); - const arma::subview_col sol = sys.getSolution(_id); + TRACE(14,"Duct::getvart()"); + const SegPositionMapper& sol = sys.getSolution(_id); us Ns = sys.Ns(); // Wraparound if(gp<0) gp+=ngp(); - us vars_per_gp = 4; - vars_per_gp+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0); + // rho,u,T,p + // Ts maybe + us nvars_per_gp = (_ductpb.stempmodel() + == pb::HeatBalance ? 5 : 4); - return sol.subvec((gp*vars_per_gp+varnr)*Ns, - (gp*vars_per_gp+varnr+1)*Ns-1); + return *sol.at((gp*nvars_per_gp+varnr)); } vd Duct::getvarx(int varnr,int t) const { vd res(ngp()); @@ -185,48 +248,44 @@ vd Duct::getvarx(int varnr,int t) const { } return res; } -vd Duct::initialSolution() const { +void Duct::initialSolution(SegPositionMapper& sol) const { TRACE(15,"Duct::initialSolution()"); - vd initsol(getNDofs()); - VARTRACE(15,initsol.size()); + VARTRACE(15,sol.size()); - us vars_per_gp = 4; - vars_per_gp+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0); - - us Ns = sys.Ns(); + us nvars_per_gp = 4; + nvars_per_gp+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0); const Gas& gas = sys.gas(); - + + us segdof = 0; + for(us i=0;i&& residual) const; - vd initialSolution() const; + // Solving + virtual void residualJac(SegPositionMapper& residual, + SegJacRows& jac) const; // virtual void getSolution(const TaSystem&,const us insertion_start) const; diff --git a/src/ductbc/adiabaticwall.cpp b/src/ductbc/adiabaticwall.cpp index b49d2e1..50a0c5b 100644 --- a/src/ductbc/adiabaticwall.cpp +++ b/src/ductbc/adiabaticwall.cpp @@ -35,20 +35,22 @@ AdiabaticWall::~AdiabaticWall() { AdiabaticWall* AdiabaticWall::copy(const TaSystem& sys) const { return new AdiabaticWall(sys,*this); } -vd AdiabaticWall::initialSolution() const { - return vd(); -} -void AdiabaticWall::residualJac(arma::subview_col&& residual - ) const { - TRACE(15,"AdiabaticWall::residual()"); +void AdiabaticWall::residualJac(SegPositionMapper& residual, + SegJacRows& jac) const { + + TRACE(15,"AdiabaticWall::residualJac()"); const Duct& duct = getDuct(); const pb::Duct& dpb = duct.getDuctPb(); us Ns = sys.Ns(); + us i=0; + if(_side == pb::left) { - residual.subvec(0,Ns-1) = duct.ut(0); + *residual.at(i) = duct.ut(0); + jac.at(i)[duct.getDof(constants::u,0)].eye(Ns,Ns); + if(dpb.htmodel() != pb::Isentropic) { // TODO: Put the boundary condition of zero heat flux here // residual.subvec(Ns,2*Ns-1) = @@ -56,7 +58,8 @@ void AdiabaticWall::residualJac(arma::subview_col&& residual } } else { - residual.subvec(0,Ns-1) = duct.ut(-1); + *residual.at(i) = duct.ut(-1); + jac.at(i)[duct.getDof(constants::u,-1)].eye(Ns,Ns); if(dpb.htmodel() != pb::Isentropic) { // TODO: Put the boundary condition of zero heat flux here // residual.subvec(Ns,2*Ns-1) = duct.Tt(sys,-1) - _T; diff --git a/src/ductbc/adiabaticwall.h b/src/ductbc/adiabaticwall.h index 648ae97..fb19848 100644 --- a/src/ductbc/adiabaticwall.h +++ b/src/ductbc/adiabaticwall.h @@ -34,10 +34,8 @@ public: ~AdiabaticWall(); AdiabaticWall* copy(const TaSystem& sys) const; - vd initialSolution() const; - - virtual void residualJac(arma::subview_col&& residual - ) const; + virtual void residualJac(SegPositionMapper& residual, + SegJacRows& jac) const; // Return the total number of equations in this segment virtual us getNEqs() const; @@ -50,11 +48,6 @@ public: // 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; - - }; #endif // ADIABATICWALL_H diff --git a/src/ductbc/ductbc.cpp b/src/ductbc/ductbc.cpp index deaa74f..44c53a4 100644 --- a/src/ductbc/ductbc.cpp +++ b/src/ductbc/ductbc.cpp @@ -31,10 +31,10 @@ DuctBc* DuctBc::newDuctBc(const us id, return new PressureBc(sys,id,dbc); break; } - // case pb::AdiabaticWall: { - // return new AdiabaticWall(id,sys,dbc); - // break; - // } + case pb::AdiabaticWall: { + return new AdiabaticWall(sys,id,dbc); + break; + } default: tasmet_assert(false,"Not implemented DuctBc"); break; diff --git a/src/ductbc/ductbc.h b/src/ductbc/ductbc.h index 41567ed..8ec6550 100644 --- a/src/ductbc/ductbc.h +++ b/src/ductbc/ductbc.h @@ -33,6 +33,8 @@ public: const Duct& getDuct() const; + virtual void initialSolution(SegPositionMapper&) const {} + }; diff --git a/src/ductbc/pressurebc.cpp b/src/ductbc/pressurebc.cpp index 25db0e6..6d527bf 100644 --- a/src/ductbc/pressurebc.cpp +++ b/src/ductbc/pressurebc.cpp @@ -69,31 +69,28 @@ PressureBc::~PressureBc() { PressureBc* PressureBc::copy(const TaSystem& sys) const { return new PressureBc(sys,*this); } -vd PressureBc::initialSolution() const { - return vd(); -} -void PressureBc::residualJac(arma::subview_col&& residual - ) const { +void PressureBc::residualJac(SegPositionMapper& residual, + SegJacRows& jacrows) const { - TRACE(15,"PressureBc::residual()"); + TRACE(15,"PressureBc::residualJac()"); const pb::Duct& dpb = getDuct().getDuctPb(); us Ns = sys.Ns(); const Duct& duct = getDuct(); if(_side == pb::left) { - residual.subvec(0,Ns-1) = duct.pt(0) - _p; - + *residual.at(0) = duct.pt(0) - _p; + jacrows.at(0)[duct.getDof(constants::p,0)] = eye(Ns,Ns); if(dpb.htmodel() != pb::Isentropic) { - residual.subvec(Ns,2*Ns-1) = duct.Tt(0) - _T; + // residual.subvec(Ns,2*Ns-1) = duct.Tt(0) - _T; } } else { - residual.subvec(0,Ns-1) = duct.pt(-1) - _p; - + *residual.at(0) = duct.pt(-1) - _p; + jacrows.at(0)[duct.getDof(constants::p,-1)] = eye(Ns,Ns); if(dpb.htmodel() != pb::Isentropic) { - residual.subvec(Ns,2*Ns-1) = duct.Tt(-1) - _T; + // residual.subvec(Ns,2*Ns-1) = duct.Tt(-1) - _T; } } } diff --git a/src/ductbc/pressurebc.h b/src/ductbc/pressurebc.h index 15fe435..6cbea2b 100644 --- a/src/ductbc/pressurebc.h +++ b/src/ductbc/pressurebc.h @@ -38,10 +38,8 @@ public: PressureBc* copy(const TaSystem&) const; - vd initialSolution() const; - - virtual void residualJac(arma::subview_col&& residual - ) const; + virtual void residualJac(SegPositionMapper& residual, + SegJacRows& jac) const; // Return the total number of equations in this segment virtual us getNEqs() const; diff --git a/src/material/gas.h b/src/material/gas.h index 60ce8b6..8c68c0e 100644 --- a/src/material/gas.h +++ b/src/material/gas.h @@ -21,11 +21,11 @@ using pb::air; using pb::helium; using pb::nitrogen; -#define element_wise(varname) \ - vd varname_(T.size()); \ +#define element_wise(function) \ + vd function_(T.size()); \ for(us i=0;i +#include class TripletList; +/** + * A PosId is a specifier for the position in the global Jacobian + * matrix. Both considering the column as well as the row + * + */ +struct PosId { + us segno; /**< This is the segment number */ + us segoffset; /**< This is the offset in the segment */ -class Jacobian{ - us _ndofs; -public: - Jacobian(us ndofs):_ndofs(ndofs){} - vector jacrows; + bool operator==(const PosId& other) const { + return (segno==other.segno) && + (segoffset == other.segoffset); + } - void operator+=(const Jacobian&); - void operator+=(const JacRow&); + bool operator<(const PosId& other) const { + if(other.segno < segno) return true; + else if(other.segno > segno) return false; + else return (other.segoffset < segoffset); - operator TripletList() const; + } +}; +/** + * To use the PosId in an unordered_map, it requires a custom hash + * function. This struct has a custom function to generate a hash from + * the given two integers. + * + */ + +struct hash_PosId +{ + std::size_t operator() ( const PosId& m ) const + { + const us a = m.segno+1 ; + const us b = m.segoffset+1 ; + return std::hash()( a*b + (a+b) ) ; + } }; - +typedef std::unordered_map JacRow; + +typedef std::vector SegJacRows; + +typedef std::map Jacobian; + +typedef std::vector SegPositionMapper; +typedef std::map GlobalPositionMapper; #endif // JACOBIAN_H ////////////////////////////////////////////////////////////////////// diff --git a/src/sys/jacrow.cpp b/src/sys/jacrow.cpp deleted file mode 100644 index 3b89247..0000000 --- a/src/sys/jacrow.cpp +++ /dev/null @@ -1,40 +0,0 @@ -#include "jacobian.h" -#include "tasmet_variable.h" -#include "globalconf.h" - - -JacRow& JacRow::operator+=(const JacCol& j){ - TRACE(10,"JacRow::operator+=(const JacCol& j)"); - jaccols.push_back(j); - return *this; -} - -JacRow& JacRow::operator*=(const d val){ - TRACE(2,"Jacobian::operator*=()"); - for(JacCol& col: jaccols) - col*=val; - return *this; -} -JacRow& JacRow::operator+=(const JacRow& other){ - // Reserve some space - this->jaccols.reserve(this->jaccols.capacity()+other.jaccols.size()-this->jaccols.size()); - - for(const JacCol& col :other.jaccols) - (*this)+=col; - - return *this; -} -JacRow JacRow::operator-() const{ - TRACE(15,"JacRow::operator-()"); - - JacRow result(*this); - for (JacCol& jaccol : result.jaccols) - jaccol*=-1; - - return result; -} -void JacRow::prePostMultiply(const dmat& pre,const dmat& post) { - for(JacCol& jaccol: jaccols) - jaccol.prePostMultiply(pre,post); -} - diff --git a/src/sys/jacrow.h b/src/sys/jacrow.h deleted file mode 100644 index 2716956..0000000 --- a/src/sys/jacrow.h +++ /dev/null @@ -1,55 +0,0 @@ -// jacrow.h -// -// Author: J.A. de Jong -// -// Description: -// -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef JACROW_H -#define JACROW_H -#include -#include "jaccol.h" - -class JacRow{ // Row in Jacobian matrix - -public: - - PosId id; /**< Specifies the row for this - JacRow */ - std::vector jaccols; /**< The columns in this row */ - - /** - * Initialize a Jacobian row with a single column block - * - * @param id The equation ID - * @param col The added column - */ - JacRow(PosId id,const JacCol& col):JacRow(id){jaccols.push_back(col);} - - JacRow(PosId id): id(id){} - - JacRow operator-() const; - - JacRow& operator+=(const JacCol&); - JacRow& operator+=(const JacRow& jacrow); - JacRow& operator*=(const d val); // Multiply all terms with constant value - - /** - * Pre- and postmultiply all Column blocks with the given - * matrices. Hence the new columns will be col_new = pre*col_old*post - * - * @param pre : The pre-multiply matrix - * @param post: The post-multiply matrix - */ - void prePostMultiply(const dmat& pre,const dmat& post); - - -}; - - - - - -#endif // JACROW_H -////////////////////////////////////////////////////////////////////// diff --git a/src/sys/segment.h b/src/sys/segment.h index 441545b..8ffd1df 100644 --- a/src/sys/segment.h +++ b/src/sys/segment.h @@ -15,9 +15,10 @@ #include "tasmet_types.h" #include "tasmet_exception.h" +#include "jacobian.h" // #include "phaseconstraint.h" -class Jacobian; + class GlobalConf; class TaSystem; @@ -51,9 +52,10 @@ public: */ virtual Segment* copy(const TaSystem&) const = 0; - virtual vd initialSolution() const = 0; + virtual void initialSolution(SegPositionMapper&) const = 0; - virtual void residualJac(arma::subview_col&& residual // Here we store the residual + virtual void residualJac(SegPositionMapper& residual, + SegJacRows& jacrows// Here we store the residual ) const=0; // Get and set name diff --git a/src/sys/tasystem.cpp b/src/sys/tasystem.cpp index 0922e18..6685b0a 100644 --- a/src/sys/tasystem.cpp +++ b/src/sys/tasystem.cpp @@ -60,47 +60,8 @@ TaSystem::TaSystem(const pb::System& sys): throw e; } } - - // Create the initial solution from the segments and store it - // here. Please be careful not to call any virtual functions! - vus ndofs = getNDofs(); - vus dofend = arma::cumsum(ndofs); - us total_dofs = arma::sum(ndofs); - vd solution = vd(total_dofs); - us i=0; - - - - const Segment* seg; - for(auto& seg_: _segs) { - seg = seg_.second; - - if(ndofs(i)>0) { - if(i==0) { - solution.subvec(0,ndofs(0)-1) = - seg->initialSolution(); - } - else { - solution.subvec(dofend(i-1),dofend(i)-1) = - seg->initialSolution(); - } - i++; - } - } - // TODO: work directly on the final solution array - _solution = solution; - - // Copy solution vector, if valid - // const auto& sol = sys.solution(); - // us size = sol.size(), i=0; - // if(size>0) { - // _solution = vd(size); - // for(auto& val: sol) { - // _solution(i) = val; - // i++; - // } - // } + initSolRes(); } TaSystem::TaSystem(const TaSystem& o): @@ -117,7 +78,89 @@ TaSystem::TaSystem(const TaSystem& o): seg.second = seg.second->copy(*this); if(!seg.second) throw TaSMETBadAlloc(); } + + initSolRes(); + } + void TaSystem::initSolRes() { + // Create the initial solution from the segments and store it + // here. Please be careful not to call any virtual functions! + vus ndofs = getNDofs(); + vus neqs = getNEqs(); + + us Ns = this->Ns(); + + us total_dofs = arma::sum(ndofs); + us total_eqs = arma::sum(neqs); + + if(total_dofs != total_eqs) { + throw TaSMETError("Number of equations does not match number of DOFS. " + "This is probably due to missing, or too much boundary conditions"); + } + + if(total_dofs*Ns > constants::maxndofs) { + stringstream error; + error << "Too many DOFS required. Problem too large. " + "Number of equations computed: "; + error << total_dofs*Ns; + throw TaSMETError(error); + } + + + _solution = zeros(total_dofs*Ns); + _residual = zeros(total_dofs*Ns)+arma::datum::nan; + + d* solptr = _solution.memptr(); + d* resptr = _residual.memptr(); + + us globalposdof=0; + us globalposeq=0; + + // Loop over the segments, and for each dof in the segment, we + // allocate a vd, with its memory pointed to from the global + // solution pointer. + + for(auto& seg_: _segs) { + + us segid = seg_.first; + VARTRACE(15,segid); + + for(us i=0;i< ndofs(segid); i++) { + + // false: own memory (means the vd does not own the memory + // true: strict true means bound to the memory for its + // lifetime. Does not allow memory adjustments. + _solution_dofs[segid].push_back(new vd(&solptr[globalposdof*Ns],Ns,false,true)); + globalposdof++; + VARTRACE(15,i); + } + + // Set the initial solution + if(ndofs(segid)>0) + seg_.second->initialSolution(_solution_dofs.at(segid)); + + + + for(us i=0;i< neqs(segid); i++) { + + // false: own memory (means the vd does not own the memory + // true: strict true means bound to the memory for its + // lifetime. Does not allow memory adjustments. + _residual_eqs[segid].push_back(new vd(&resptr[globalposeq*Ns],Ns,false,true)); + + _jac[segid].push_back(JacRow()); + + globalposeq++; + } + + } + + tasmet_assert(globalposdof == total_dofs,"Bug in assigning solution vectors"); + tasmet_assert(globalposeq == total_eqs,"Bug in assigning residual vectors"); + + } + + 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 @@ -156,56 +199,85 @@ int TaSystem::getArbitrateMassEq(const vus& neqs) const { void TaSystem::residualJac(ResidualJac& resjac) const { TRACE(15,"TaSystem::residual()"); - vus neqs = getNEqs(); - us total_neqs = arma::sum(neqs); - - if(total_neqs>constants::maxndofs) { - stringstream error; - error << "Too many DOFS required. Problem too large. Number of equations computed: "; - error << total_neqs; - throw TaSMETError(error); - } - - // This vector of indices stores the last equation number + 1 for - // each equation set in a Segment - vus eqsend = arma::cumsum(neqs); - - int arbitrateMassEq = getArbitrateMassEq(neqs); + int arbitrateMassEq = -1;//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,eqsend); - - vd residual(total_neqs); #ifdef TASMET_DEBUG - residual = arma::datum::nan*ones(total_neqs); + _residual *= arma::datum::nan; #endif - us i=0; + us segid; const Segment* seg; for(auto seg_: _segs) { - + segid = seg_.first; seg = seg_.second; - if(i==0) { // Put the residual of the segment in the over-all residual - seg->residualJac(residual.subvec(0,eqsend(0)-1)); - } - else { - seg->residualJac(residual.subvec(eqsend(i-1),eqsend(i)-1)); - } + seg->residualJac(_residual_eqs[segid],_jac[segid]); // Count the mass, add it if(arbitrateMassEq!=-1) { mass += seg->getMass(); } - i++; - } + } + + us Ns = this->Ns(); + vus neqs = getNEqs(); + vus neqs_sum = arma::cumsum(neqs); + us total_neqs = arma::sum(neqs); + + vus ndofs = getNDofs(); + vus ndofs_sum = arma::cumsum(ndofs); + us total_ndofs = arma::sum(ndofs); + + // Now we create the triplets + TripletList triplets(total_neqs*Ns); + + for(auto& segjacrows : _jac) { + + // The current segment number + us segno = segjacrows.first; + + // The row offset for this segment + us rowoffset = (segno==0) ? 0 : neqs_sum(segno-1)*Ns; + + for(auto& jacrow: segjacrows.second) { + + for(auto& jaccol : jacrow) { + const PosId& id = jaccol.first; + + us coloffset = (id.segno==0) ? 0 : ndofs_sum(segno-1)*Ns; + coloffset += id.segoffset; + coloffset*=Ns; + + dmat& colmat = jaccol.second; + + for(us i=0;i=0) { + _residual(arbitrateMassEq)=mass - _mass; } - resjac.residual = residual; + resjac.jacobian = sdmat(triplets); + resjac.residual = _residual; + } vd TaSystem::getSolution() const { - return _solution; } -const arma::subview_col TaSystem::getSolution(const us seg_id) const { - - vus ndofs = getNDofs(); - vus dofsend = arma::cumsum(ndofs); - VARTRACE(15,dofsend); - VARTRACE(15,seg_id); - if(seg_id == 0) { - VARTRACE(15,_solution.size()); - return _solution.subvec(0,dofsend(0)-1); - } - else { - return _solution.subvec(dofsend(seg_id-1),dofsend(seg_id)-1); - } +const SegPositionMapper& TaSystem::getSolution(const us seg_id) const { + return _solution_dofs.at(seg_id); } vus TaSystem::getNDofs() const { TRACE(0,"TaSystem::getNDofs()"); vus Ndofs(_segs.size()); - us i=0; for (auto seg : _segs) { - Ndofs(i)=seg.second->getNDofs(); - i++; + Ndofs(seg.first)=seg.second->getNDofs(); } return Ndofs; } vus TaSystem::getNEqs() const { TRACE(15,"TaSystem::getNEqs()"); vus Neqs(_segs.size()); - us i=0; for (auto seg :_segs) { - Neqs(i)=seg.second->getNEqs(); - VARTRACE(15,Neqs(i)); - i++; + Neqs(seg.first)=seg.second->getNEqs(); } return Neqs; } @@ -341,7 +398,19 @@ TaSystem::~TaSystem() { cleanup(); } void TaSystem::cleanup() { + purge(_segs); + for(auto& segdofs: _solution_dofs) { + for(vd* dof : segdofs.second){ + delete dof; + } + } + for(auto& segdofs: _residual_eqs) { + for(vd* dof : segdofs.second){ + delete dof; + } + } + } ////////////////////////////////////////////////////////////////////// diff --git a/src/sys/tasystem.h b/src/sys/tasystem.h index fe06374..d091905 100644 --- a/src/sys/tasystem.h +++ b/src/sys/tasystem.h @@ -11,9 +11,9 @@ #include "system.h" #include "globalconf.h" #include "triplets.h" +#include "jacobian.h" #include "gas.h" #include "protobuf/system.pb.h" -#include #include class Segment; @@ -28,21 +28,48 @@ protected: std::unique_ptr _gas; + mutable 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.*/ + mutable + GlobalPositionMapper _solution_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 + + mutable + vd _residual; + + mutable + GlobalPositionMapper _residual_eqs; + + mutable + Jacobian _jac; /**< Storing the Jacobian brings the + hugh benefit of only allocating the + memory once. Once this is done, all + values can be updated during each + iteration. */ TaSystem& operator=(const TaSystem& other)=delete; + +private: + /** + * Initializes the _solution and _solution_dofs vector. This + * should be done during construction of the object. The method is + * called from all constructors. Should be done after construction + * of the segments. + * + */ + void initSolRes(); public: TaSystem(const TaSystem& o); TaSystem(const pb::System&); + + + static TaSystem testSystem() { return TaSystem(pb::System::default_instance()); } @@ -67,7 +94,7 @@ public: vd getSolution() const; // Obtain the solution vector for the Segment with given id - const arma::subview_col getSolution(const us seg_id) const; + const SegPositionMapper& getSolution(const us seg_id) const; virtual void updateSolution(const vd& sol) {_solution = sol; } // Update the solution