Nektar++
Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | List of all members
Nektar::PulseWaveSystem Class Reference

Base class for unsteady solvers. More...

#include <PulseWaveSystem.h>

Inheritance diagram for Nektar::PulseWaveSystem:
[legend]

Public Member Functions

virtual ~PulseWaveSystem ()
 Destructor. More...
 
int GetNdomains ()
 
Array< OneD, MultiRegions::ExpListSharedPtr > UpdateVessels (void)
 
void CalcCharacteristicVariables (int omega)
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
virtual SOLVER_UTILS_EXPORT ~UnsteadySystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble >> &inarray)
 Calculate the larger time-step mantaining the problem stable. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void SetUpTraceNormals (void)
 
SOLVER_UTILS_EXPORT void InitObject ()
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise ()
 Perform any initialisation necessary before solving the problem. More...
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem. More...
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space. More...
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space. More...
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve. More...
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name. More...
 
template<class T >
std::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name. More...
 
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession ()
 Get Session name. More...
 
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure ()
 Get pressure field if available. More...
 
SOLVER_UTILS_EXPORT void ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 
SOLVER_UTILS_EXPORT void PrintSummary (std::ostream &out)
 Print a summary of parameters and solver characteristics. More...
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda. More...
 
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction (std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
 Get a SessionFunction by name. More...
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields. More...
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the L2 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > ErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write checkpoint file of custom data fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write input fields to the given filename. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
 Input field data from the given file to multiple domains. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, std::vector< std::string > &fieldStr, Array< OneD, Array< OneD, NekDouble > > &coeffs)
 Output a field. Input field data into array from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, MultiRegions::ExpListSharedPtr &pField, std::string &pFieldName)
 Output a field. Input field data into ExpList from the given file. More...
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary. More...
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMap & UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetFinalTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp ()
 Virtual function to identify if operator is negated in DoSolve. More...
 

Protected Member Functions

 PulseWaveSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises PulseWaveSystem class members. More...
 
virtual void v_InitObject ()
 
virtual void v_DoInitialise ()
 Sets up initial conditions. More...
 
virtual void v_DoSolve ()
 Solves an unsteady problem. More...
 
void LinkSubdomains (Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &fields)
 Links the subdomains. More...
 
void BifurcationRiemann (Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0)
 Riemann Problem for Bifurcation. More...
 
void MergingRiemann (Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0)
 Riemann Problem for Merging Flow. More...
 
void JunctionRiemann (Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0)
 Riemann Problem for Junction. More...
 
virtual void v_Output (void)
 
void CheckPoint_Output (const int n)
 
NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Compute the L2 error between fields and a given exact solution. More...
 
NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Compute the L_inf error between fields and a given exact solution. More...
 
void WriteVessels (const std::string &outname)
 Write input fields to the given filename. More...
 
void EnforceInterfaceConditions (const Array< OneD, const Array< OneD, NekDouble > > &fields)
 
- Protected Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises UnsteadySystem class members. More...
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. More...
 
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &s)
 Print a summary of time stepping parameters. More...
 
virtual SOLVER_UTILS_EXPORT void v_AppendOutput1D (Array< OneD, Array< OneD, NekDouble >> &solution1D)
 Print the solution at each solution point in a txt file. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_GetTimeStep (const Array< OneD, const Array< OneD, NekDouble >> &inarray)
 Return the timestep to be used for the next step in the time-marching loop. More...
 
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_RequireFwdTrans ()
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time, int &nchk)
 
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff (const Array< OneD, Array< OneD, NekDouble >> vel, StdRegions::VarCoeffMap &varCoeffMap)
 Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h time velocity. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys ()
 Virtual function for transformation to physical space. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff ()
 Virtual function for transformation to coefficient space. More...
 
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution (unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure (void)
 
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 

Protected Attributes

Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
 
int m_nDomains
 
int m_currentDomain
 
int m_nVariables
 
UpwindTypePulse m_upwindTypePulse
 
Array< OneD, int > m_fieldPhysOffset
 
NekDouble m_rho
 
NekDouble m_pext
 
NekDouble m_C
 
NekDouble m_RT
 
NekDouble m_pout
 
Array< OneD, Array< OneD, NekDouble > > m_A_0
 
Array< OneD, Array< OneD, NekDouble > > m_A_0_trace
 
Array< OneD, Array< OneD, NekDouble > > m_beta
 
Array< OneD, Array< OneD, NekDouble > > m_beta_trace
 
Array< OneD, Array< OneD, NekDouble > > m_trace_fwd_normal
 
std::vector< std::vector< InterfacePointShPtr > > m_vesselJcts
 
std::vector< std::vector< InterfacePointShPtr > > m_bifurcations
 
std::vector< std::vector< InterfacePointShPtr > > m_mergingJcts
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
int m_infosteps
 Number of time steps between outputting status information. More...
 
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
int m_nanSteps
 
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
 
NekDouble m_epsilon
 
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used. More...
 
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used. More...
 
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used. More...
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check. More...
 
std::ofstream m_errFile
 
std::vector< int > m_intVariables
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
NekDouble m_filterTimeWarning
 Number of time steps between outputting status information. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtr > m_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error. More...
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous) More...
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous) More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
int m_HomoDirec
 number of homogenous directions More...
 

Private Member Functions

void SetUpDomainInterfaces (void)
 
void FillDataFromInterfacePoint (InterfacePointShPtr &I, const Array< OneD, const Array< OneD, NekDouble > > &field, NekDouble &A, NekDouble &u, NekDouble &beta, NekDouble &A_0)
 

Additional Inherited Members

- Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 
- Protected Types inherited from Nektar::SolverUtils::EquationSystem
enum  HomogeneousType { eHomogeneous1D, eHomogeneous2D, eHomogeneous3D, eNotHomogeneous }
 Parameter for homogeneous expansions. More...
 
- Static Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
static std::string equationSystemTypeLookupIds []
 

Detailed Description

Base class for unsteady solvers.

Initialises the arterial subdomains in m_vessels and sets up all domain-linking conditions (bifurcations, junctions, merging flows). Detects the network structure and assigns boundary conditons. Also provides the underlying timestepping framework for pulse wave solvers including the general timestepping routines.

Definition at line 85 of file PulseWaveSystem.h.

Constructor & Destructor Documentation

◆ ~PulseWaveSystem()

Nektar::PulseWaveSystem::~PulseWaveSystem ( )
virtual

Destructor.

Destructor

Definition at line 74 of file PulseWaveSystem.cpp.

75  {
76  }

◆ PulseWaveSystem()

Nektar::PulseWaveSystem::PulseWaveSystem ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const SpatialDomains::MeshGraphSharedPtr &  pGraph 
)
protected

Initialises PulseWaveSystem class members.

Processes SolverInfo parameters from the session file and sets up timestepping-specific code.

Parameters
m_SessionSession object to read parameters from.

Definition at line 64 of file PulseWaveSystem.cpp.

67  : UnsteadySystem(pSession, pGraph)
68  {
69  }
Nektar::SolverUtils::UnsteadySystem::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.
Definition: UnsteadySystem.cpp:68

Member Function Documentation

◆ BifurcationRiemann()

void Nektar::PulseWaveSystem::BifurcationRiemann ( Array< OneD, NekDouble > &  Au,
Array< OneD, NekDouble > &  uu,
Array< OneD, NekDouble > &  beta,
Array< OneD, NekDouble > &  A_0 
)
protected

Riemann Problem for Bifurcation.

Solves the Riemann problem at a bifurcation by assuming subsonic flow at both sides of the boundary and by applying conservation of mass and continuity of the total pressure \( \frac{p}{rho} + \frac{u^{2}}{2}. \) The other 3 missing equations come from the characteristic variables. For further information see "Pulse WavePropagation in the human vascular system" Section 3.4.4

Definition at line 644 of file PulseWaveSystem.cpp.

References ASSERTL0, and m_rho.

Referenced by EnforceInterfaceConditions().

646  {
647  NekDouble rho = m_rho;
648  Array<OneD, NekDouble> W(3);
649  Array<OneD, NekDouble> W_Au(3);
650  Array<OneD, NekDouble> P_Au(3);
651  Array<OneD, NekDouble> f(6);
652  Array<OneD, NekDouble> g(6);
653  Array<OneD, NekDouble> tmp(6);
654  Array<OneD, Array<OneD, NekDouble> > inv_J(6);
655  for (int i=0; i<6; i++)
656  {
657  inv_J[i] = Array<OneD, NekDouble> (6);
658  }
659  NekDouble k = 0.0;
660  NekDouble k1 = 0.0;
661  NekDouble k2 = 0.0;
662  NekDouble k3 = 0.0;
663 
664  int proceed = 1;
665  int iter = 0;
666  int MAX_ITER = 7;
667 
668  // Calculated from input
669  W[0] = uu[0] + 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
670  W[1] = uu[1] - 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
671  W[2] = uu[2] - 4*sqrt(beta[2]/(2*rho))*(sqrt(sqrt(Au[2])) - sqrt(sqrt(A_0[2])));
672 
673  // Tolerances for the algorithm
674  NekDouble Tol = 1.0e-10;
675 
676  // Newton Iteration
677  while ((proceed) && (iter < MAX_ITER))
678  {
679  iter = iter+1;
680 
681  // Calculate the constraint vector, six equations:
682  // 3 characteristic variables, mass conservation,
683  // total pressure
684  W_Au[0] = 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
685  W_Au[1] = 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
686  W_Au[2] = 4*sqrt(beta[2]/(2*rho))*(sqrt(sqrt(Au[2])) - sqrt(sqrt(A_0[2])));
687 
688  P_Au[0] = beta[0]*(sqrt(Au[0]) - sqrt(A_0[0]));
689  P_Au[1] = beta[1]*(sqrt(Au[1]) - sqrt(A_0[1]));
690  P_Au[2] = beta[2]*(sqrt(Au[2]) - sqrt(A_0[2]));
691 
692  f[0] = uu[0] + W_Au[0] - W[0];
693  f[1] = uu[1] - W_Au[1] - W[1];
694  f[2] = uu[2] - W_Au[2] - W[2];
695  f[3] = Au[0]*uu[0] - Au[1]*uu[1] - Au[2]*uu[2];
696  f[4] = uu[0]*uu[0] + 2.0/rho*P_Au[0] - uu[1]*uu[1] - 2.0/rho*P_Au[1];
697  f[5] = uu[0]*uu[0] + 2.0/rho*P_Au[0] - uu[2]*uu[2] - 2.0/rho*P_Au[2];
698 
699  // Calculate the wave speed at each vessel
700  NekDouble c1 = sqrt(beta[0]/(2*rho))*sqrt(sqrt(Au[0]));
701  NekDouble c2 = sqrt(beta[1]/(2*rho))*sqrt(sqrt(Au[1]));
702  NekDouble c3 = sqrt(beta[2]/(2*rho))*sqrt(sqrt(Au[2]));
703 
704  // Inverse Jacobian matrix J(x[n])^(-1), is
705  // already inverted here analytically
706  k = c1*Au[1]*c3+Au[0]*c3*c2+Au[2]*c1*c2;
707  k1 = (c1-uu[0])*k;
708  inv_J[0][0] = (-c2*uu[0]*c3*Au[0]+Au[2]*c2*c1*c1+Au[1]*c1*c1*c3)/k1;
709  inv_J[0][1] = Au[1]*(c2-uu[1])*c1*c3/k1;
710  inv_J[0][2] = Au[2]*(c3-uu[2])*c1*c2/k1;
711  inv_J[0][3] = c1*c2*c3/k1;
712  inv_J[0][4] = -0.5*c1*Au[1]*c3/k1;
713  inv_J[0][5] = -0.5*Au[2]*c1*c2/k1;
714 
715  k2 = (c2+uu[1])*k;
716  inv_J[1][0] = Au[0]*(c1+uu[0])*c2*c3/k2;
717  inv_J[1][1] = (c1*uu[1]*c3*Au[1]+Au[2]*c1*c2*c2+c3*c2*c2*Au[0])/k2;
718  inv_J[1][2] = -Au[2]*(c3-uu[2])*c1*c2/k2;
719  inv_J[1][3] = -c1*c2*c3/k2;
720  inv_J[1][4] = -0.5*(c1*Au[2]+Au[0]*c3)*c2/k2;
721  inv_J[1][5] = 0.5*Au[2]*c1*c2/k2;
722 
723  k3 = (c3+uu[2])*k;
724  inv_J[2][0] = Au[0]*(c1+uu[0])*c2*c3/k3;
725  inv_J[2][1] = -Au[1]*(c2-uu[1])*c1*c3/k3;
726  inv_J[2][2] = (c1*c2*uu[2]*Au[2]+c1*Au[1]*c3*c3+c2*c3*c3*Au[0])/k3;
727  inv_J[2][3] = -c1*c2*c3/k3;
728  inv_J[2][4] = 0.5*c1*Au[1]*c3/k3;
729  inv_J[2][5] = -0.5*(Au[1]*c1+c2*Au[0])*c3/k3;
730 
731  inv_J[3][0] = Au[0]*(Au[0]*c3*c2-uu[0]*c3*Au[1]-uu[0]*c2*Au[2])/k1;
732  inv_J[3][1] = -Au[0]*Au[1]*(c2-uu[1])*c3/k1;
733  inv_J[3][2] = -Au[0]*Au[2]*(c3-uu[2])*c2/k1;
734  inv_J[3][3] = -Au[0]*c3*c2/k1;
735  inv_J[3][4] = 0.5*Au[0]*Au[1]*c3/k1;
736  inv_J[3][5] = 0.5*Au[0]*c2*Au[2]/k1;
737 
738  inv_J[4][0] = Au[0]*Au[1]*(c1+uu[0])*c3/k2;
739  inv_J[4][1] = -Au[1]*(c1*Au[1]*c3+c1*uu[1]*Au[2]+c3*uu[1]*Au[0])/k2;
740  inv_J[4][2] = -Au[2]*Au[1]*(c3-uu[2])*c1/k2;
741  inv_J[4][3] = -c1*Au[1]*c3/k2;
742  inv_J[4][4] = -0.5*Au[1]*(c1*Au[2]+Au[0]*c3)/k2;
743  inv_J[4][5] = 0.5*Au[2]*Au[1]*c1/k2;
744 
745  inv_J[5][0] = Au[0]*Au[2]*(c1+uu[0])*c2/k3;
746  inv_J[5][1] = -Au[2]*Au[1]*(c2-uu[1])*c1/k3;
747  inv_J[5][2] = -Au[2]*(Au[2]*c1*c2+c1*uu[2]*Au[1]+c2*uu[2]*Au[0])/k3;
748  inv_J[5][3] = -Au[2]*c1*c2/k3;
749  inv_J[5][4] = 0.5*Au[2]*Au[1]*c1/k3;
750  inv_J[5][5] = -0.5*Au[2]*(Au[1]*c1+c2*Au[0])/k3;
751 
752 
753  // Solve the system by multiplying the Jacobian with the vector f:
754  // g = (inv_J)*f
755  for (int j=0; j<6; j++)
756  {
757  tmp[j] =0.0;
758  g[j] = 0.0;
759  }
760 
761  for (int j=0; j<6; j++)
762  {
763  for (int i=0; i<6; i++)
764  {
765  tmp[j] = inv_J[j][i]*f[i];
766  g[j] += tmp[j];
767  }
768  }
769 
770  // Update the solution: x_new = x_old - dx
771  uu[0] = uu[0] - g[0];
772  uu[1] = uu[1] - g[1];
773  uu[2] = uu[2] - g[2];
774  Au[0] = Au[0] - g[3];
775  Au[1] = Au[1] - g[4];
776  Au[2] = Au[2] - g[5];
777 
778  // Check if the error of the solution is smaller than Tol
779  if ((g[0]*g[0] + g[1]*g[1] + g[2]*g[2] + g[3]*g[3]+ g[4]*g[4] + g[5]*g[5]) < Tol)
780  {
781  proceed = 0;
782  }
783 
784  // Check if solver converges
785  if (iter >= MAX_ITER)
786  {
787  ASSERTL0(false,"Riemann solver for Bifurcation did not converge");
788  }
789  }
790  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43

◆ CalcCharacteristicVariables()

void Nektar::PulseWaveSystem::CalcCharacteristicVariables ( int  omega)

Definition at line 1281 of file PulseWaveSystem.cpp.

References m_beta, m_rho, m_vessels, Vmath::Smul(), Vmath::Vadd(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vsub().

1282  {
1283  int nq = m_vessels[omega]->GetTotPoints();
1284  Array<OneD, NekDouble> A = m_vessels[omega]->UpdatePhys();
1285  Array<OneD, NekDouble> u = m_vessels[omega+1]->UpdatePhys();
1286  Array<OneD, NekDouble> c(nq);
1287 
1288  //Calc 4*c
1289  Vmath::Vsqrt(nq,A,1,c,1);
1290  Vmath::Vmul(nq,m_beta[omega],1,c,1,c,1);
1291  Vmath::Smul(nq,16.0/(2*m_rho),c,1,c,1);
1292  Vmath::Vsqrt(nq,c,1,c,1);
1293 
1294  // Characteristics
1295  Vmath::Vadd(nq,u,1,c,1,A,1);
1296  Vmath::Vsub(nq,u,1,c,1,u,1);
1297  }
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:411
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Nektar::PulseWaveSystem::m_beta
Array< OneD, Array< OneD, NekDouble > > m_beta
Definition: PulseWaveSystem.h:120
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ CheckPoint_Output()

void Nektar::PulseWaveSystem::CheckPoint_Output ( const int  n)
protected

Writes the .fld file at the end of the simulation. Similar to the normal v_Output however the Multidomain output has to be prepared.

Definition at line 1090 of file PulseWaveSystem.cpp.

References Nektar::SolverUtils::EquationSystem::m_sessionName, and WriteVessels().

Referenced by v_DoSolve().

1091  {
1092  std::stringstream outname;
1093  outname << m_sessionName << "_" << n << ".chk";
1094 
1095  WriteVessels(outname.str());
1096  }
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:345
Nektar::PulseWaveSystem::WriteVessels
void WriteVessels(const std::string &outname)
Write input fields to the given filename.
Definition: PulseWaveSystem.cpp:1103

◆ EnforceInterfaceConditions()

void Nektar::PulseWaveSystem::EnforceInterfaceConditions ( const Array< OneD, const Array< OneD, NekDouble > > &  fields)
protected

Definition at line 555 of file PulseWaveSystem.cpp.

References BifurcationRiemann(), FillDataFromInterfacePoint(), JunctionRiemann(), m_bifurcations, m_mergingJcts, m_nVariables, m_vesselJcts, m_vessels, and MergingRiemann().

Referenced by Nektar::PulseWavePropagation::DoOdeRhs().

556  {
557  int dom, bcpos;
558  Array<OneD, NekDouble> Au(3),uu(3),beta(3),A_0(3);
559 
560 
561  // Enfore Bifurcations;
562  for(int n = 0; n < m_bifurcations.size(); ++n)
563  {
564  FillDataFromInterfacePoint(m_bifurcations[n][0],fields,
565  Au[0],uu[0],beta[0],A_0[0]);
566  FillDataFromInterfacePoint(m_bifurcations[n][1],fields,
567  Au[1],uu[1],beta[1],A_0[1]);
568  FillDataFromInterfacePoint(m_bifurcations[n][2],fields,
569  Au[2],uu[2],beta[2],A_0[2]);
570 
571  // Solve the Riemann problem for a bifurcation
572  BifurcationRiemann(Au, uu, beta, A_0);
573 
574  // Store the values into the right positions:
575  for(int i = 0; i < 3; ++i)
576  {
577  dom = m_bifurcations[n][i]->m_domain;
578  bcpos = m_bifurcations[n][i]->m_bcPosition;
579  m_vessels[dom*m_nVariables] ->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = Au[i];
580  m_vessels[dom*m_nVariables+1]->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = uu[i];
581  }
582  }
583 
584  // Enfore Bifurcations;
585  for(int n = 0; n < m_mergingJcts.size(); ++n)
586  {
587  // Merged vessel
588  FillDataFromInterfacePoint(m_mergingJcts[n][0],fields,
589  Au[0],uu[0],beta[0],A_0[0]);
590  FillDataFromInterfacePoint(m_mergingJcts[n][1],fields,
591  Au[1],uu[1],beta[1],A_0[1]);
592  FillDataFromInterfacePoint(m_mergingJcts[n][2],fields,
593  Au[2],uu[2],beta[2],A_0[2]);
594 
595  // Solve the Riemann problem for a merging vessel
596  MergingRiemann(Au, uu, beta, A_0);
597 
598  // Store the values into the right positions:
599  for(int i = 0; i < 3; ++i)
600  {
601  int dom = m_mergingJcts[n][i]->m_domain;
602  int bcpos = m_mergingJcts[n][i]->m_bcPosition;
603  m_vessels[dom*m_nVariables] ->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = Au[i];
604  m_vessels[dom*m_nVariables+1]->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = uu[i];
605  }
606  }
607 
608  for(int n = 0; n < m_vesselJcts.size(); ++n)
609  {
610 
611  FillDataFromInterfacePoint(m_vesselJcts[n][0],fields,
612  Au[0],uu[0],beta[0],A_0[0]);
613  FillDataFromInterfacePoint(m_vesselJcts[n][1],fields,
614  Au[1],uu[1],beta[1],A_0[1]);
615 
616  JunctionRiemann(Au, uu, beta, A_0);
617 
618  // Store the values into the right positions:
619  for(int i = 0; i < 2; ++i)
620  {
621  int dom = m_vesselJcts[n][i]->m_domain;
622  int bcpos = m_vesselJcts[n][i]->m_bcPosition;
623  m_vessels[dom*m_nVariables] ->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = Au[i];
624  m_vessels[dom*m_nVariables+1]->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = uu[i];
625  //cout<<"Au: "<<Au[i]<<endl;
626  //cout<<"uu: "<<uu[i]<<endl;
627  //cout<<endl;
628  }
629 
630  }
631  }
Nektar::PulseWaveSystem::BifurcationRiemann
void BifurcationRiemann(Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0)
Riemann Problem for Bifurcation.
Definition: PulseWaveSystem.cpp:644
Nektar::PulseWaveSystem::m_vesselJcts
std::vector< std::vector< InterfacePointShPtr > > m_vesselJcts
Definition: PulseWaveSystem.h:125
Nektar::PulseWaveSystem::MergingRiemann
void MergingRiemann(Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0)
Riemann Problem for Merging Flow.
Definition: PulseWaveSystem.cpp:803
Nektar::PulseWaveSystem::m_bifurcations
std::vector< std::vector< InterfacePointShPtr > > m_bifurcations
Definition: PulseWaveSystem.h:126
Nektar::PulseWaveSystem::FillDataFromInterfacePoint
void FillDataFromInterfacePoint(InterfacePointShPtr &I, const Array< OneD, const Array< OneD, NekDouble > > &field, NekDouble &A, NekDouble &u, NekDouble &beta, NekDouble &A_0)
Definition: PulseWaveSystem.cpp:529
Nektar::PulseWaveSystem::JunctionRiemann
void JunctionRiemann(Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0)
Riemann Problem for Junction.
Definition: PulseWaveSystem.cpp:960
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104
Nektar::PulseWaveSystem::m_mergingJcts
std::vector< std::vector< InterfacePointShPtr > > m_mergingJcts
Definition: PulseWaveSystem.h:127

◆ FillDataFromInterfacePoint()

void Nektar::PulseWaveSystem::FillDataFromInterfacePoint ( InterfacePointShPtr &  I,
const Array< OneD, const Array< OneD, NekDouble > > &  field,
NekDouble &  A,
NekDouble &  u,
NekDouble &  beta,
NekDouble &  A_0 
)
private

Definition at line 529 of file PulseWaveSystem.cpp.

References m_A_0_trace, m_beta_trace, m_fieldPhysOffset, m_nVariables, and m_vessels.

Referenced by EnforceInterfaceConditions().

533  {
534  int omega, traceId, eid, vert, phys_offset, vesselID;
535 
536 
537  // Parent vessel
538  omega = I->m_domain;
539  traceId = I->m_traceId;
540  eid = I->m_elmt;
541  vert = I->m_elmtVert;
542  vesselID = omega*m_nVariables;
543 
544  phys_offset = m_vessels[vesselID]->GetPhys_Offset(eid);
545 
546  m_vessels[vesselID]->GetExp(eid)->
547  GetVertexPhysVals(vert, fields[0]+m_fieldPhysOffset[omega]+phys_offset, A);
548  m_vessels[vesselID]->GetExp(eid)->
549  GetVertexPhysVals(vert, fields[1]+m_fieldPhysOffset[omega]+phys_offset, u);
550 
551  beta = m_beta_trace[omega][traceId];
552  A_0 = m_A_0_trace [omega][traceId];
553  }
Nektar::PulseWaveSystem::m_fieldPhysOffset
Array< OneD, int > m_fieldPhysOffset
Definition: PulseWaveSystem.h:110
Nektar::PulseWaveSystem::m_beta_trace
Array< OneD, Array< OneD, NekDouble > > m_beta_trace
Definition: PulseWaveSystem.h:121
Nektar::PulseWaveSystem::m_A_0_trace
Array< OneD, Array< OneD, NekDouble > > m_A_0_trace
Definition: PulseWaveSystem.h:119
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104

◆ GetNdomains()

int Nektar::PulseWaveSystem::GetNdomains ( )
inline

Definition at line 91 of file PulseWaveSystem.h.

92  {
93  return m_nDomains;
94  }

◆ JunctionRiemann()

void Nektar::PulseWaveSystem::JunctionRiemann ( Array< OneD, NekDouble > &  Au,
Array< OneD, NekDouble > &  uu,
Array< OneD, NekDouble > &  beta,
Array< OneD, NekDouble > &  A_0 
)
protected

Riemann Problem for Junction.

Solves the Riemann problem at an interdomain junction by assuming subsonic flow at both sides of the boundary and by applying conservation of mass and continuity of the total pressure \( \frac{p}{rho} + \frac{u^{2}}{2}. \) The other 2 missing equations come from the characteristic variables. For further information see "Pulse WavePropagation in the human vascular system" Section 3.4.

Definition at line 960 of file PulseWaveSystem.cpp.

References ASSERTL0, and m_rho.

Referenced by EnforceInterfaceConditions().

962  {
963  NekDouble rho = m_rho;
964  Array<OneD, NekDouble> W(2);
965  Array<OneD, NekDouble> W_Au(2);
966  Array<OneD, NekDouble> P_Au(2);
967  Array<OneD, NekDouble> f(4);
968  Array<OneD, NekDouble> g(4);
969  Array<OneD, NekDouble> tmp(4);
970  Array<OneD, Array<OneD, NekDouble> > inv_J(4);
971  for (int i=0; i<4; i++)
972  {
973  inv_J[i] = Array<OneD, NekDouble> (4);
974  }
975  NekDouble k = 0.0;
976  NekDouble k1 = 0.0;
977  NekDouble k2 = 0.0;
978 
979  int proceed = 1;
980  int iter = 0;
981  int MAX_ITER = 7;
982  NekDouble Tol = 1.0e-10;
983 
984  // Calculated from input
985  W[0] = uu[0] + 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
986  W[1] = uu[1] - 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
987 
988  while((proceed) && (iter < MAX_ITER))
989  {
990  iter = iter+1;
991 
992  // Calculate the constraint vector, 4 equations:
993  // 2 characteristic variables, mass conservation,
994  // total pressure
995  W_Au[0] = 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
996  W_Au[1] = 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
997 
998  P_Au[0] = beta[0]*(sqrt(Au[0]) - sqrt(A_0[0]));
999  P_Au[1] = beta[1]*(sqrt(Au[1]) - sqrt(A_0[1]));
1000 
1001  f[0] = uu[0] + W_Au[0] - W[0];
1002  f[1] = uu[1] - W_Au[1] - W[1];
1003  f[2] = Au[0]*uu[0] - Au[1]*uu[1];
1004  f[3] = uu[0]*uu[0] + 2.0/rho*P_Au[0] - uu[1]*uu[1] - 2.0/rho*P_Au[1];
1005 
1006  // Calculate the wave speed at each vessel
1007  NekDouble cl = sqrt(beta[0]/(2*rho))*sqrt(sqrt(Au[0]));
1008  NekDouble cr = sqrt(beta[1]/(2*rho))*sqrt(sqrt(Au[1]));
1009 
1010  // Inverse Jacobian matrix J(x[n])^(-1), is already inverted here analytically
1011  k = (cl*Au[1]+Au[0]*cr);
1012  k1 = (cl-uu[0])*k;
1013  inv_J[0][0] = (Au[1]*cl*cl-cr*uu[0]*Au[0])/k1;
1014  inv_J[0][1] = Au[1]*(cr-uu[1])*cl/k1;
1015  inv_J[0][2] = cl*cr/k1;
1016  inv_J[0][3] = -0.5*cl*Au[1]/k1;
1017 
1018  k2 = (cr+uu[1])*k;
1019  inv_J[1][0] = Au[0]*(cl+uu[0])*cr/k2;
1020  inv_J[1][1] = (cl*uu[1]*Au[1]+cr*cr*Au[0])/k2;
1021  inv_J[1][2] = -cl*cr/k2;
1022  inv_J[1][3] = -0.5*Au[0]*cr/k2;
1023 
1024  inv_J[2][0] = Au[0]*(Au[0]*cr-uu[0]*Au[1])/k1;
1025  inv_J[2][1] = -Au[0]*Au[1]*(cr-uu[1])/k1;
1026  inv_J[2][2] = -Au[0]*cr/k1;
1027  inv_J[2][3] = 0.5*Au[1]*Au[0]/k1;
1028 
1029  inv_J[3][0] = Au[0]*Au[1]*(cl+uu[0])/k2;
1030  inv_J[3][1] = -Au[1]*(cl*Au[1]+uu[1]*Au[0])/k2;
1031  inv_J[3][2] = -cl*Au[1]/k2;
1032  inv_J[3][3] = -0.5*Au[1]*Au[0]/k2;
1033 
1034  // Solve the system by multiplying the Jacobian with the vector f:
1035  // g = (inv_J)*f
1036  for (int j=0; j<4; j++)
1037  {
1038  tmp[j] =0.0;
1039  g[j] = 0.0;
1040  }
1041 
1042  for (int j=0; j<4; j++)
1043  {
1044  for (int i=0; i<4; i++)
1045  {
1046  tmp[j] = inv_J[j][i]*f[i];
1047  g[j] += tmp[j];
1048  }
1049  }
1050 
1051  // Update solution: x_new = x_old - dx
1052  uu[0] -= g[0];
1053  uu[1] -= g[1];
1054  Au[0] -= g[2];
1055  Au[1] -= g[3];
1056 
1057  // Check if the error of the solution is smaller than Tol.
1058  if((g[0]*g[0] + g[1]*g[1] + g[2]*g[2] + g[3]*g[3]) < Tol)
1059  proceed = 0;
1060  }
1061 
1062  if(iter >= MAX_ITER)
1063  {
1064  ASSERTL0(false,"Riemann solver for Junction did not converge");
1065  }
1066  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43

◆ LinkSubdomains()

void Nektar::PulseWaveSystem::LinkSubdomains ( Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  fields)
protected

Links the subdomains.

◆ MergingRiemann()

void Nektar::PulseWaveSystem::MergingRiemann ( Array< OneD, NekDouble > &  Au,
Array< OneD, NekDouble > &  uu,
Array< OneD, NekDouble > &  beta,
Array< OneD, NekDouble > &  A_0 
)
protected

Riemann Problem for Merging Flow.

Solves the Riemann problem at an merging flow condition by assuming subsonic flow at both sides of the boundary and by applying conservation of mass and continuity of the total pressure \( \frac{p}{rho} + \frac{u^{2}}{2}. \) The other 3 missing equations come from the characteristic variables. For further information see "Pulse WavePropagation in the human vascular system" Section 3.4.4

Definition at line 803 of file PulseWaveSystem.cpp.

References ASSERTL0, and m_rho.

Referenced by EnforceInterfaceConditions().

804  {
805  NekDouble rho = m_rho;
806  Array<OneD, NekDouble> W(3);
807  Array<OneD, NekDouble> W_Au(3);
808  Array<OneD, NekDouble> P_Au(3);
809  Array<OneD, NekDouble> f(6);
810  Array<OneD, NekDouble> g(6);
811  Array<OneD, NekDouble> tmp(6);
812  Array<OneD, Array<OneD, NekDouble> > inv_J(6);
813 
814  for (int i=0; i<6; i++)
815  {
816  inv_J[i] = Array<OneD, NekDouble> (6);
817  }
818 
819  NekDouble k = 0.0;
820  NekDouble k1 = 0.0;
821  NekDouble k2 = 0.0;
822  NekDouble k3 = 0.0;
823 
824  int proceed = 1;
825  int iter = 0;
826  int MAX_ITER = 7;
827 
828  // Calculated from input
829  W[0] = uu[0] - 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
830  W[1] = uu[1] + 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
831  W[2] = uu[2] + 4*sqrt(beta[2]/(2*rho))*(sqrt(sqrt(Au[2])) - sqrt(sqrt(A_0[2])));
832 
833  // Tolerances for the algorithm
834  NekDouble Tol = 1.0e-10;
835 
836  // Newton Iteration
837  while ((proceed) && (iter < MAX_ITER))
838  {
839  iter = iter+1;
840 
841  // Calculate the constraint vector, six equations:
842  // 3 characteristic variables, mass conservation,
843  // total pressure
844  W_Au[0] = 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
845  W_Au[1] = 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
846  W_Au[2] = 4*sqrt(beta[2]/(2*rho))*(sqrt(sqrt(Au[2])) - sqrt(sqrt(A_0[2])));
847 
848  P_Au[0] = beta[0]*(sqrt(Au[0]) - sqrt(A_0[0]));
849  P_Au[1] = beta[1]*(sqrt(Au[1]) - sqrt(A_0[1]));
850  P_Au[2] = beta[2]*(sqrt(Au[2]) - sqrt(A_0[2]));
851 
852  f[0] = uu[0] - W_Au[0] - W[0];
853  f[1] = uu[1] + W_Au[1] - W[1];
854  f[2] = uu[2] + W_Au[2] - W[2];
855  f[3] = Au[0]*uu[0] - Au[1]*uu[1] - Au[2]*uu[2];
856  f[4] = uu[0]*uu[0] + 2.0/rho*P_Au[0] - uu[1]*uu[1] - 2.0/rho*P_Au[1];
857  f[5] = uu[0]*uu[0] + 2.0/rho*P_Au[0] - uu[2]*uu[2] - 2.0/rho*P_Au[2];
858 
859  // Calculate the wave speed at each vessel
860  NekDouble c1 = sqrt(beta[0]/(2*rho))*sqrt(sqrt(Au[0]));
861  NekDouble c2 = sqrt(beta[1]/(2*rho))*sqrt(sqrt(Au[1]));
862  NekDouble c3 = sqrt(beta[2]/(2*rho))*sqrt(sqrt(Au[2]));
863 
864  // Inverse Jacobian matrix J(x[n])^(-1), is already inverted here analytically
865  k = c1*Au[1]*c3+Au[0]*c3*c2+Au[2]*c1*c2;
866  k1 = (c1+uu[0])*k;
867  inv_J[0][0] = (c2*uu[0]*c3*Au[0]+Au[2]*c2*c1*c1+Au[1]*c1*c1*c3)/k1;
868  inv_J[0][1] = Au[1]*(c2+uu[1])*c1*c3/k1;
869  inv_J[0][2] = Au[2]*(c3+uu[2])*c1*c2/k1;
870  inv_J[0][3] = c1*c2*c3/k1;
871  inv_J[0][4] = 0.5*Au[1]*c1*c3/k1;
872  inv_J[0][5] = 0.5*Au[2]*c1*c2/k1;
873 
874  k2 = (c2-uu[1])*k;
875  inv_J[1][0] = Au[0]*(c1-uu[0])*c2*c3/k2;
876  inv_J[1][1] = (-c1*uu[1]*c3*Au[1]+Au[2]*c1*c2*c2+c3*c2*c2*Au[0])/k2;
877  inv_J[1][2] = -Au[2]*(c3+uu[2])*c1*c2/k2;
878  inv_J[1][3] = -c1*c2*c3/k2;
879  inv_J[1][4] = 0.5*(c1*Au[2]+Au[0]*c3)*c2/k2;
880  inv_J[1][5] = -0.5*Au[2]*c1*c2/k2;
881 
882  k3 = (c3-uu[2])*k;
883  inv_J[2][0] = Au[0]*(c1-uu[0])*c2*c3/k3;
884  inv_J[2][1] = -Au[1]*(c2+uu[1])*c1*c3/k3;
885  inv_J[2][2] = -(c1*uu[2]*c2*Au[2]-Au[1]*c1*c3*c3-c2*c3*c3*Au[0])/k3;
886  inv_J[2][3] = -c1*c2*c3/k3;
887  inv_J[2][4] = -0.5*Au[1]*c1*c3/k3;
888  inv_J[2][5] = 0.5*(Au[1]*c1+Au[0]*c2)*c3/k3;
889 
890  inv_J[3][0] = -Au[0]*(Au[0]*c3*c2+uu[0]*c3*Au[1]+uu[0]*c2*Au[2])/k1;
891  inv_J[3][1] = Au[0]*Au[1]*(c2+uu[1])*c3/k1;
892  inv_J[3][2] = Au[0]*Au[2]*(c3+uu[2])*c2/k1;
893  inv_J[3][3] = Au[0]*c3*c2/k1;
894  inv_J[3][4] = 0.5*Au[0]*Au[1]*c3/k1;
895  inv_J[3][5] = 0.5*Au[0]*c2*Au[2]/k1;
896 
897  inv_J[4][0] = -Au[0]*Au[1]*(c1-uu[0])*c3/k2;
898  inv_J[4][1] = Au[1]*(Au[1]*c1*c3-c1*uu[1]*Au[2]-c3*uu[1]*Au[0])/k2;
899  inv_J[4][2] = Au[2]*Au[1]*(c3+uu[2])*c1/k2;
900  inv_J[4][3] = Au[1]*c1*c3/k2;
901  inv_J[4][4] = -0.5*Au[1]*(c1*Au[2]+Au[0]*c3)/k2;
902  inv_J[4][5] = 0.5*Au[2]*Au[1]*c1/k2;
903 
904  inv_J[5][0] = -Au[0]*Au[2]*(c1-uu[0])*c2/k3;
905  inv_J[5][1] = Au[2]*Au[1]*(c2+uu[1])*c1/k3;
906  inv_J[5][2] = Au[2]*(Au[2]*c1*c2-c1*uu[2]*Au[1]-c2*uu[2]*Au[0])/k3;
907  inv_J[5][3] = Au[2]*c1*c2/k3;
908  inv_J[5][4] = 0.5*Au[2]*Au[1]*c1/k3;
909  inv_J[5][5] = -0.5*Au[2]*(Au[1]*c1+Au[0]*c2)/k3;
910 
911  // Solve the system by multiplying the Jacobian with the vector f:
912  // g = (inv_J)*f
913  for (int j=0; j<6; j++)
914  {
915  tmp[j] =0.0;
916  g[j] = 0.0;
917  }
918 
919  for (int j=0; j<6; j++)
920  {
921  for (int i=0; i<6; i++)
922  {
923  tmp[j] = inv_J[j][i]*f[i];
924  g[j] += tmp[j];
925  }
926  }
927 
928  // Update the solution: x_new = x_old - dx
929  uu[0] = uu[0] - g[0];
930  uu[1] = uu[1] - g[1];
931  uu[2] = uu[2] - g[2];
932  Au[0] = Au[0] - g[3];
933  Au[1] = Au[1] - g[4];
934  Au[2] = Au[2] - g[5];
935 
936  // Check if the error of the solution is smaller than Tol
937  if ((g[0]*g[0] + g[1]*g[1] + g[2]*g[2] + g[3]*g[3]+ g[4]*g[4] + g[5]*g[5]) < Tol)
938  {
939  proceed = 0;
940  }
941 
942  // Check if solver converges
943  if (iter >= MAX_ITER)
944  {
945  ASSERTL0(false,"Riemann solver for Merging Flow did not converge");
946  }
947 
948  }
949  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43

◆ SetUpDomainInterfaces()

void Nektar::PulseWaveSystem::SetUpDomainInterfaces ( void  )
private

Definition at line 255 of file PulseWaveSystem.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNotDefined, Nektar::ErrorUtil::ewarning, m_bifurcations, m_mergingJcts, m_nDomains, m_nVariables, m_vesselJcts, m_vessels, NEKERROR, and CellMLToNektar.cellml_metadata::p.

Referenced by v_InitObject().

256  {
257  map<int,std::vector<InterfacePointShPtr> > VidToDomain;
258 
259  // loop over domain and find out if we have any undefined
260  // boundary conditions representing interfaces. If so make a
261  // map based around vid and storing the domains that are
262  // part of interfaces.
263  for(int omega = 0; omega < m_nDomains; ++omega)
264  {
265  int vesselID = omega*m_nVariables;
266 
267  for(int i = 0; i < 2; ++i)
268  {
269  if(m_vessels[vesselID]->GetBndConditions()[i]->GetBoundaryConditionType() == SpatialDomains::eNotDefined)
270  {
271  // Get Vid of interface
272  int vid = m_vessels[vesselID]->UpdateBndCondExpansion(i)->GetExp(0)->GetGeom()->GetVid(0);
273  cout<<"Shared vertex id: "<<vid<<endl;
274  MultiRegions::ExpListSharedPtr trace = m_vessels[vesselID]->GetTrace();
275  InterfacePointShPtr Ipt;
276 
277  bool finish = false;
278  // find which elmt, the lcoal vertex and the data offset of point
279  for(int n = 0; n < m_vessels[vesselID]->GetExpSize(); ++n)
280  {
281  for(int p = 0; p < 2; ++p)
282  {
283  if(m_vessels[vesselID]->GetTraceMap()->GetElmtToTrace()[n][p]->as<LocalRegions::Expansion>()->GetGeom()->GetVid(0) == vid)
284  {
285  int eid = m_vessels[vesselID]->GetTraceMap()->GetElmtToTrace()[n][p]->GetElmtId();
286 
287  int tid = m_vessels[vesselID]->GetTrace()->GetCoeff_Offset(eid);
288  Ipt = MemoryManager<InterfacePoint>::AllocateSharedPtr(vid,omega,n,p,tid,i);
289 
290  cout<<"Global Vid of interface point: "<<vid<<endl;
291  cout<<"Domain interface point belongs to: "<<omega<<endl;
292  cout<<"Element id of vertex: "<<n<<endl;
293  cout<<"Vertex id within local element: "<<p<<endl;
294  cout<<"Element id within the trace: "<<tid<<endl;
295  cout<<"Position of boundary condition in region: "<<i<<endl;
296 
297  finish = true;
298  break;
299  }
300  }
301  if(finish == true)
302  {
303  break;
304  }
305  }
306 
307  VidToDomain[vid].push_back(Ipt);
308 
309  // finally reset boundary condition to Dirichlet
310  m_vessels[vesselID]->GetBndConditions()[i]
311  ->SetBoundaryConditionType(SpatialDomains::eDirichlet);
312 
313  m_vessels[vesselID+1]->GetBndConditions()[i]
314  ->SetBoundaryConditionType(SpatialDomains::eDirichlet);
315  }
316  }
317  }
318 
319  // loop over map and set up Interface information;
320  for(auto &iter : VidToDomain)
321  {
322  if(iter.second.size() == 2) // Vessel jump interface
323  {
324  m_vesselJcts.push_back(iter.second);
325  }
326  else if(iter.second.size() == 3) // Bifurcation or Merging junction.
327  {
328  int nbeg = 0;
329  int nend = 0;
330 
331  // Determine if bifurcation or merging junction
332  // through number of elemnt vertices that meet at
333  // junction. Only one vertex using a m_elmtVert=1
334  // indicates a bifurcation
335  for(int i = 0; i < 3; ++i)
336  {
337  if(iter.second[i]->m_elmtVert == 0)
338  {
339  nbeg += 1;
340  }
341  else
342  {
343  nend += 1;
344  }
345  }
346 
347  // Set up Bifurcation information
348  if(nbeg == 2)
349  {
350  // ensure first InterfacePoint is parent
351  if(iter.second[0]->m_elmtVert == 1) //m_elmtVert: Vertex id in local element
352  {
353  m_bifurcations.push_back(iter.second);
354  }
355  else
356  {
357  //order points according to Riemann solver convention
358  InterfacePointShPtr I;
359  //find merging vessel
360  if(iter.second[1]->m_elmtVert == 1)
361  {
362  I = iter.second[0];
363  iter.second[0] = iter.second[1];
364  iter.second[1] = I;
365  }
366  else if (iter.second[2]->m_elmtVert == 1)
367  {
368  I = iter.second[0];
369  iter.second[0] = iter.second[2];
370  iter.second[2] = I;
371  }
372  NEKERROR(ErrorUtil::ewarning,"This routine has not been checked");
373  }
374  }
375  else
376  {
377  // ensure last InterfacePoint is merged vessel
378  if(iter.second[0]->m_elmtVert == 0)
379  {
380  m_mergingJcts.push_back(iter.second);
381  }
382  else
383  {
384  //order points according to Riemann solver convention
385  InterfacePointShPtr I;
386  //find merging vessel
387  if(iter.second[1]->m_elmtVert == 0)
388  {
389  I = iter.second[0];
390  iter.second[0] = iter.second[1];
391  iter.second[1] = I;
392  }
393  else if (iter.second[2]->m_elmtVert == 0)
394  {
395  I = iter.second[0];
396  iter.second[0] = iter.second[2];
397  iter.second[2] = I;
398  }
399  NEKERROR(ErrorUtil::ewarning,"This routine has not been checked");
400  }
401  }
402 
403  }
404  else
405  {
406  ASSERTL0(false,"Unknown junction type");
407  }
408  }
409  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:209
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::PulseWaveSystem::m_vesselJcts
std::vector< std::vector< InterfacePointShPtr > > m_vesselJcts
Definition: PulseWaveSystem.h:125
Nektar::PulseWaveSystem::m_bifurcations
std::vector< std::vector< InterfacePointShPtr > > m_bifurcations
Definition: PulseWaveSystem.h:126
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::InterfacePointShPtr
std::shared_ptr< InterfacePoint > InterfacePointShPtr
Definition: PulseWaveSystem.h:82
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104
Nektar::PulseWaveSystem::m_mergingJcts
std::vector< std::vector< InterfacePointShPtr > > m_mergingJcts
Definition: PulseWaveSystem.h:127

◆ UpdateVessels()

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::PulseWaveSystem::UpdateVessels ( void  )
inline

Definition at line 96 of file PulseWaveSystem.h.

97  {
98  return m_vessels;
99  }
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104

◆ v_DoInitialise()

void Nektar::PulseWaveSystem::v_DoInitialise ( void  )
protectedvirtual

Sets up initial conditions.

Initialisation routine for multiple subdomain case. Sets the initial conditions for all arterial subdomains read from the inputfile. Sets the material properties and the A_0 area for all subdomains and fills the domain-linking boundary conditions with the initial values of their domain.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 419 of file PulseWaveSystem.cpp.

References Nektar::SolverUtils::EquationSystem::m_fields, m_nDomains, m_nVariables, Nektar::SolverUtils::EquationSystem::m_session, m_vessels, and Nektar::SolverUtils::EquationSystem::SetInitialConditions().

420  {
421 
422  if (m_session->GetComm()->GetRank() == 0)
423  {
424  cout << "Initial Conditions:" << endl;
425  }
426 
427  /* Loop over all subdomains to initialize all with the Initial
428  * Conditions read from the inputfile*/
429  for (int omega = 0; omega < m_nDomains; omega++)
430  {
431  m_fields[0] = m_vessels[m_nVariables*omega];
432  m_fields[1] = m_vessels[m_nVariables*omega+1];
433 
434  if (m_session->GetComm()->GetRank() == 0)
435  {
436  cout << "Subdomain = " <<omega<<endl;
437  }
438 
439  SetInitialConditions(0.0,0,omega);
440  }
441  // Reset to first definition
442  m_fields[0] = m_vessels[0];
443  m_fields[1] = m_vessels[1];
444 
445  }
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:333
Nektar::SolverUtils::EquationSystem::SetInitialConditions
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
Definition: EquationSystem.h:620
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104

◆ v_DoSolve()

void Nektar::PulseWaveSystem::v_DoSolve ( void  )
protectedvirtual

Solves an unsteady problem.

NEEDS Updating:

DoSolve routine for PulseWavePropagation with multiple subdomains taken from UnsteadySystem and modified for multidomain case. Initialises the time integration scheme (as specified in the session file), and perform the time integration. Within the timestepping loop the following is done: 1. Link all arterial segments according to the network structure, solve the Riemann problem between different arterial segments and assign the values to the boundary conditions (LinkSubdomains) 2. Every arterial segment is solved independentl for this timestep. This is done by handing the solution vector \( \mathbf{u} \) and the right hand side m_ode, which is the PulseWavePropagation class in this example over to the time integration scheme

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 464 of file PulseWaveSystem.cpp.

References CheckPoint_Output(), Nektar::SolverUtils::EquationSystem::m_checksteps, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::UnsteadySystem::m_infosteps, Nektar::SolverUtils::UnsteadySystem::m_intScheme, Nektar::SolverUtils::UnsteadySystem::m_intSoln, m_nDomains, m_nVariables, Nektar::SolverUtils::UnsteadySystem::m_ode, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_steps, Nektar::SolverUtils::EquationSystem::m_time, Nektar::SolverUtils::EquationSystem::m_timestep, m_vessels, Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), Nektar::LibUtilities::Timer::TimePerTest(), and Vmath::Vcopy().

465  {
466  NekDouble IntegrationTime = 0.0;
467  int i,n,nchk = 1;
468 
469  Array<OneD, Array<OneD,NekDouble> > fields(m_nVariables);
470 
471  for(int i = 0; i < m_nVariables; ++i)
472  {
473  fields[i] = m_vessels[i]->UpdatePhys();
474  m_fields[i]->SetPhysState(false);
475  }
476 
477  m_intSoln = m_intScheme->InitializeScheme(
478  m_timestep,fields,m_time,m_ode);
479 
480  // Time loop
481  for(n = 0; n < m_steps; ++n)
482  {
483  LibUtilities::Timer timer;
484  timer.Start();
485  fields = m_intScheme->TimeIntegrate(n,m_timestep,m_intSoln,m_ode);
486  //cout<<"integration: "<<fields[0][fields[0].num_elements()-1]<<endl;
487  m_time += m_timestep;
488  timer.Stop();
489  IntegrationTime += timer.TimePerTest(1);
490 
491  // Write out status information.
492  if(m_session->GetComm()->GetRank() == 0 && !((n+1)%m_infosteps))
493  {
494  cout << "Steps: " << n+1
495  << "\t Time: " << m_time
496  << "\t Time-step: " << m_timestep << "\t" << endl;
497  }
498 
499  // Transform data if needed
500  if(!((n+1)%m_checksteps))
501  {
502  for (i = 0; i < m_nVariables; ++i)
503  {
504  int cnt = 0;
505  for (int omega = 0; omega < m_nDomains; omega++)
506  {
507  m_vessels[omega*m_nVariables+i]->FwdTrans(fields[i]+cnt,
508  m_vessels[omega*m_nVariables+i]->UpdateCoeffs());
509  cnt += m_vessels[omega*m_nVariables+i]->GetTotPoints();
510  }
511  }
512  CheckPoint_Output(nchk++);
513  }
514 
515  }//end of timeintegration
516 
517  //Copy Array To Vessel Phys Fields
518  for(int i = 0; i < m_nVariables; ++i)
519  {
520  Vmath::Vcopy(fields[i].num_elements(), fields[i],1,m_vessels[i]->UpdatePhys(),1);
521  }
522 
523  cout <<"Time-integration timing : "
524  << IntegrationTime << " s" << endl << endl;
525  }
Nektar::PulseWaveSystem::CheckPoint_Output
void CheckPoint_Output(const int n)
Definition: PulseWaveSystem.cpp:1090
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:347
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:71
Nektar::SolverUtils::EquationSystem::m_timestep
NekDouble m_timestep
Time step size.
Definition: EquationSystem.h:353
Nektar::SolverUtils::EquationSystem::m_checksteps
int m_checksteps
Number of steps between checkpoints.
Definition: EquationSystem.h:363
Nektar::SolverUtils::EquationSystem::m_steps
int m_steps
Number of steps to take.
Definition: EquationSystem.h:361
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:333
Nektar::SolverUtils::UnsteadySystem::m_intScheme
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
Wrapper to the time integration scheme.
Definition: UnsteadySystem.h:69
Nektar::SolverUtils::UnsteadySystem::m_infosteps
int m_infosteps
Number of time steps between outputting status information.
Definition: UnsteadySystem.h:62
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104
Nektar::SolverUtils::UnsteadySystem::m_intSoln
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
Definition: UnsteadySystem.h:73

◆ v_InitObject()

void Nektar::PulseWaveSystem::v_InitObject ( )
protectedvirtual

Initialisation routine for multidomain solver. Sets up the expansions for every arterial segment (m_vessels) and for one complete field m_outfield which is needed to write the postprocessing output. Also determines which upwind strategy is used (currently only upwinding scheme available) and reads blodd flow specific parameters from the inputfile 

Gets the Material Properties of each arterial segment
specified in the inputfile from section MaterialProperties

Also gets the Area at static equilibrium A_0 specified in the inputfile.

Having found these points also extract the values at the trace points and the normal direction consistent with the left adjacent definition of Fwd and Bwd

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Reimplemented in Nektar::PulseWavePropagation, and Nektar::PulseWaveSystemOutput.

Definition at line 87 of file PulseWaveSystem.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::MultiRegions::eDiscontinuous, Nektar::SolverUtils::EquationSystem::GetFunction(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), m_A_0, m_A_0_trace, m_beta, m_beta_trace, m_fieldPhysOffset, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_graph, m_nDomains, m_nVariables, m_pext, Nektar::SolverUtils::EquationSystem::m_projectionType, m_rho, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_time, m_trace_fwd_normal, m_upwindTypePulse, m_vessels, SetUpDomainInterfaces(), Nektar::SIZE_UpwindTypePulse, Nektar::UpwindTypeMapPulse, Nektar::SolverUtils::UnsteadySystem::v_InitObject(), and Nektar::SolverUtils::EquationSystem::ZeroPhysFields().

Referenced by Nektar::PulseWaveSystemOutput::v_InitObject(), and Nektar::PulseWavePropagation::v_InitObject().

88  {
89  // Initialise base class
90  UnsteadySystem::v_InitObject();
91 
92  // Read the geometry and the expansion information
93  m_nDomains = m_graph->GetDomain().size();
94 
95  // Determine projectiontype
96  ASSERTL0(m_session->MatchSolverInfo("Projection","DisContinuous"),"Pulse solver only set up for Discontinuous projections");
97  m_projectionType = MultiRegions::eDiscontinuous;
98  ASSERTL0(m_graph->GetMeshDimension() == 1,"Pulse wave solver only set up for expansion dimension equal to 1");
99 
100  int i;
101  m_nVariables = m_session->GetVariables().size();
102 
103  m_fields = Array<OneD, MultiRegions::ExpListSharedPtr> (m_nVariables);
104  m_vessels = Array<OneD, MultiRegions::ExpListSharedPtr> (m_nVariables*m_nDomains);
105 
106  const std::vector<SpatialDomains::CompositeMap> domain = m_graph->GetDomain();
107 
108  SpatialDomains::BoundaryConditions Allbcs(m_session, m_graph);
109 
110  // Set up domains and put geometry to be only one space dimension.
111  int cnt = 0;
112  bool SetToOneSpaceDimension = true;
113 
114  if(m_session->DefinesCmdLineArgument("SetToOneSpaceDimension"))
115  {
116  std::string cmdline = m_session->GetCmdLineArgument<std::string>("SetToOneSpaceDimension");
117  if(boost::to_upper_copy(cmdline) == "FALSE")
118  {
119  SetToOneSpaceDimension = false;
120  }
121  }
122 
123  for(i = 0 ; i < m_nDomains; ++i)
124  {
125  for(int j = 0; j < m_nVariables; ++j)
126  {
127  m_vessels[cnt++] = MemoryManager<MultiRegions::DisContField1D>
128  ::AllocateSharedPtr(m_session, m_graph, domain[i],
129  Allbcs,
130  m_session->GetVariable(j),
131  SetToOneSpaceDimension);
132  }
133  }
134 
135  // Reset coeff and phys space to be continuous over all domains
136  int totcoeffs = 0;
137  int totphys = 0;
138  m_fieldPhysOffset = Array<OneD, int>(m_nDomains+1,0);
139  for(i = 0; i < m_nDomains; ++i)
140  {
141  totcoeffs += m_vessels[i*m_nVariables]->GetNcoeffs();
142 
143  m_fieldPhysOffset[i] = totphys;
144  totphys += m_vessels[i*m_nVariables]->GetTotPoints();
145  }
146  m_fieldPhysOffset[m_nDomains] = totphys;
147 
148  for(int n = 0; n < m_nVariables; ++n)
149  {
150  Array<OneD, NekDouble> coeffs(totcoeffs,0.0);
151  Array<OneD, NekDouble> phys(totphys,0.0);
152  Array<OneD, NekDouble> tmpcoeffs,tmpphys;
153 
154  m_vessels[n]->SetCoeffsArray(coeffs);
155  m_vessels[n]->SetPhysArray(phys);
156 
157  int cnt = m_vessels[n]->GetNcoeffs();
158  int cnt1 = m_vessels[n]->GetTotPoints();
159 
160  for(i = 1; i < m_nDomains; ++i)
161  {
162  m_vessels[i*m_nVariables+n]->SetCoeffsArray(tmpcoeffs = coeffs + cnt);
163  m_vessels[i*m_nVariables+n]->SetPhysArray (tmpphys = phys + cnt1);
164  cnt += m_vessels[i*m_nVariables+n]->GetNcoeffs();
165  cnt1 += m_vessels[i*m_nVariables+n]->GetTotPoints();
166  }
167  }
168 
169  m_fields[0] = m_vessels[0];
170  m_fields[1] = m_vessels[1];
171 
172  // Zero all physical fields initially.
173  ZeroPhysFields();
174 
175  // If Discontinuous Galerkin determine upwinding method to use
176  for (int i = 0; i < (int)SIZE_UpwindTypePulse; ++i)
177  {
178  bool match;
179  m_session->MatchSolverInfo("UPWINDTYPEPULSE", UpwindTypeMapPulse[i], match, false);
180  if (match)
181  {
182  m_upwindTypePulse = (UpwindTypePulse) i;
183  break;
184  }
185  }
186 
187  // Load blood density
188  m_session->LoadParameter("rho", m_rho, 0.5);
189  // Load external pressure
190  m_session->LoadParameter("pext", m_pext, 0.0);
191 
192  int nq = 0;
193  /**
194  * Gets the Material Properties of each arterial segment
195  * specified in the inputfile from section MaterialProperties
196  * Also gets the Area at static equilibrium A_0 specified in the
197  * inputfile.
198  *
199  * Having found these points also extract the values at the
200  * trace points and the normal direction consistent with the
201  * left adjacent definition of Fwd and Bwd
202  */
203  m_beta = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
204  m_beta_trace = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
205  m_A_0 = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
206  m_A_0_trace = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
207  m_trace_fwd_normal = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
208 
209  for (int omega = 0; omega < m_nDomains; omega++)
210  {
211  nq = m_vessels[2*omega]->GetNpoints();
212  m_fields[0] = m_vessels[2*omega];
213 
214  m_beta[omega] = Array<OneD, NekDouble>(nq);
215  GetFunction("MaterialProperties")->Evaluate("beta", m_beta[omega], m_time, omega);
216 
217  m_A_0[omega] = Array<OneD, NekDouble>(nq);
218  GetFunction("A_0")->Evaluate("A_0", m_A_0[omega], m_time, omega);
219 
220  int nqTrace = GetTraceTotPoints();
221 
222  m_beta_trace[omega] = Array<OneD, NekDouble>(nqTrace);
223  m_fields[0]->ExtractTracePhys(m_beta[omega],m_beta_trace[omega]);
224 
225  m_A_0_trace[omega] = Array<OneD, NekDouble>(nqTrace);
226  m_fields[0]->ExtractTracePhys(m_A_0[omega],m_A_0_trace[omega]);
227 
228 
229  if(SetToOneSpaceDimension)
230  {
231  m_trace_fwd_normal[omega] = Array<OneD, NekDouble>(nqTrace,0.0);
232 
233  MultiRegions::ExpListSharedPtr trace = m_fields[0]->GetTrace();
234  int nelmt_trace = trace->GetExpSize();
235 
236  Array<OneD, Array<OneD, NekDouble> > normals(nelmt_trace);
237 
238  for(int i = 0 ; i < nelmt_trace; ++i)
239  {
240  normals[i] = m_trace_fwd_normal[omega]+i;
241  }
242 
243  // need to set to 1 for consistency since boundary
244  // conditions may not have coordim=1
245  trace->GetExp(0)->GetGeom()->SetCoordim(1);
246 
247  trace->GetNormals(normals);
248  }
249  }
250 
251  SetUpDomainInterfaces();
252 
253  }
Nektar::PulseWaveSystem::m_fieldPhysOffset
Array< OneD, int > m_fieldPhysOffset
Definition: PulseWaveSystem.h:110
Nektar::PulseWaveSystem::m_beta_trace
Array< OneD, Array< OneD, NekDouble > > m_beta_trace
Definition: PulseWaveSystem.h:121
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::PulseWaveSystem::m_upwindTypePulse
UpwindTypePulse m_upwindTypePulse
Definition: PulseWaveSystem.h:108
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:347
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
Nektar::SolverUtils::EquationSystem::m_projectionType
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:387
Nektar::UpwindTypeMapPulse
const char *const UpwindTypeMapPulse[]
Definition: PulseWaveSystem.h:52
Nektar::PulseWaveSystem::m_trace_fwd_normal
Array< OneD, Array< OneD, NekDouble > > m_trace_fwd_normal
Definition: PulseWaveSystem.h:122
Nektar::PulseWaveSystem::m_A_0_trace
Array< OneD, Array< OneD, NekDouble > > m_A_0_trace
Definition: PulseWaveSystem.h:119
Nektar::PulseWaveSystem::m_A_0
Array< OneD, Array< OneD, NekDouble > > m_A_0
Definition: PulseWaveSystem.h:118
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::SIZE_UpwindTypePulse
Length of enum list.
Definition: PulseWaveSystem.h:49
Nektar::SolverUtils::UnsteadySystem::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:80
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:677
Nektar::SolverUtils::EquationSystem::GetFunction
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
Definition: EquationSystem.cpp:688
Nektar::PulseWaveSystem::m_beta
Array< OneD, Array< OneD, NekDouble > > m_beta
Definition: PulseWaveSystem.h:120
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:333
Nektar::SolverUtils::EquationSystem::ZeroPhysFields
SOLVER_UTILS_EXPORT void ZeroPhysFields()
Definition: EquationSystem.cpp:1060
Nektar::SolverUtils::EquationSystem::m_graph
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Definition: EquationSystem.h:343
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104

◆ v_L2Error()

NekDouble Nektar::PulseWaveSystem::v_L2Error ( unsigned int  field,
const Array< OneD, NekDouble > &  exactsoln = NullNekDouble1DArray,
bool  Normalised = false 
)
protectedvirtual

Compute the L2 error between fields and a given exact solution.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1149 of file PulseWaveSystem.cpp.

References ASSERTL0, Nektar::SolverUtils::EquationSystem::GetFunction(), Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::SolverUtils::EquationSystem::m_comm, m_nDomains, Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError, m_nVariables, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_time, m_vessels, and Nektar::LibUtilities::ReduceSum.

1152  {
1153  NekDouble L2error = 0.0;
1154  NekDouble L2error_dom;
1155  NekDouble Vol = 0.0;
1156 
1157  if(m_NumQuadPointsError == 0)
1158  {
1159  for (int omega = 0; omega < m_nDomains; omega++)
1160  {
1161  int vesselid = field + omega*m_nVariables;
1162 
1163  if(m_vessels[vesselid]->GetPhysState() == false)
1164  {
1165  m_vessels[vesselid]->BwdTrans(m_vessels[vesselid]->GetCoeffs(),
1166  m_vessels[vesselid]->UpdatePhys());
1167  }
1168 
1169  if(exactsoln.num_elements())
1170  {
1171  L2error_dom = m_vessels[vesselid]->L2(
1172  m_vessels[vesselid]->GetPhys(),
1173  exactsoln);
1174  }
1175  else if (m_session->DefinesFunction("ExactSolution"))
1176  {
1177  Array<OneD, NekDouble> exactsoln(m_vessels[vesselid]->GetNpoints());
1178 
1179  LibUtilities::EquationSharedPtr vEqu
1180  = m_session->GetFunction("ExactSolution",field,omega);
1181  GetFunction("ExactSolution")->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
1182 
1183  L2error_dom = m_vessels[vesselid]->L2(
1184  m_vessels[vesselid]->GetPhys(),
1185  exactsoln);
1186 
1187  }
1188  else
1189  {
1190  L2error_dom = m_vessels[vesselid]->L2(
1191  m_vessels[vesselid]->GetPhys());
1192  }
1193 
1194  L2error += L2error_dom*L2error_dom;
1195 
1196  if(Normalised == true)
1197  {
1198  Array<OneD, NekDouble> one(m_vessels[vesselid]->GetNpoints(), 1.0);
1199 
1200  Vol += m_vessels[vesselid]->PhysIntegral(one);
1201  }
1202 
1203  }
1204  }
1205  else
1206  {
1207  ASSERTL0(false,"Not set up");
1208  }
1209 
1210 
1211  if(Normalised == true)
1212  {
1213  m_comm->AllReduce(Vol, LibUtilities::ReduceSum);
1214 
1215  L2error = sqrt(L2error/Vol);
1216  }
1217  else
1218  {
1219  L2error = sqrt(L2error);
1220  }
1221 
1222  return L2error;
1223  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:347
Nektar::SolverUtils::EquationSystem::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: EquationSystem.h:331
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::LibUtilities::EquationSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:131
Nektar::SolverUtils::EquationSystem::GetFunction
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
Definition: EquationSystem.cpp:688
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:712
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:333
Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Definition: EquationSystem.h:396
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104

◆ v_LinfError()

NekDouble Nektar::PulseWaveSystem::v_LinfError ( unsigned int  field,
const Array< OneD, NekDouble > &  exactsoln = NullNekDouble1DArray 
)
protectedvirtual

Compute the L_inf error between fields and a given exact solution.

Compute the error in the L_inf-norm

Parameters
fieldThe field to compare.
exactsolnThe exact solution to compare with.
Returns
Error in the L_inft-norm.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1232 of file PulseWaveSystem.cpp.

References ASSERTL0, Nektar::SolverUtils::EquationSystem::GetFunction(), Nektar::SolverUtils::EquationSystem::GetNpoints(), m_nDomains, Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError, m_nVariables, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_time, and m_vessels.

1234  {
1235  NekDouble LinferrorDom, Linferror = -1.0;
1236 
1237  for (int omega = 0; omega < m_nDomains; omega++)
1238  {
1239  int vesselid = field + omega*m_nVariables;
1240 
1241  if(m_NumQuadPointsError == 0)
1242  {
1243  if(m_vessels[vesselid]->GetPhysState() == false)
1244  {
1245  m_vessels[vesselid]->BwdTrans(m_vessels[vesselid]->GetCoeffs(),
1246  m_vessels[vesselid]->UpdatePhys());
1247  }
1248 
1249  if(exactsoln.num_elements())
1250  {
1251  LinferrorDom = m_vessels[vesselid]->Linf(
1252  m_vessels[vesselid]->GetPhys(),
1253  exactsoln);
1254  }
1255  else if (m_session->DefinesFunction("ExactSolution"))
1256  {
1257  Array<OneD, NekDouble> exactsoln(m_vessels[vesselid]->GetNpoints());
1258 
1259  GetFunction("ExactSolution")->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
1260 
1261  LinferrorDom = m_vessels[vesselid]->Linf(
1262  m_vessels[vesselid]->GetPhys(),
1263  exactsoln);
1264  }
1265  else
1266  {
1267  LinferrorDom = 0.0;
1268  }
1269 
1270  Linferror = (Linferror > LinferrorDom)? Linferror:LinferrorDom;
1271 
1272  }
1273  else
1274  {
1275  ASSERTL0(false,"ErrorExtraPoints not allowed for this solver");
1276  }
1277  }
1278  return Linferror;
1279  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:347
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::EquationSystem::GetFunction
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
Definition: EquationSystem.cpp:688
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:712
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:333
Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Definition: EquationSystem.h:396
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104

◆ v_Output()

void Nektar::PulseWaveSystem::v_Output ( void  )
protectedvirtual

Writes the .fld file at the end of the simulation. Similar to the normal v_Output however the Multidomain output has to be prepared.

Write the field data to file. The file is named according to the session name with the extension .fld appended.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1074 of file PulseWaveSystem.cpp.

References Nektar::SolverUtils::EquationSystem::m_sessionName, and WriteVessels().

1075  {
1076  /**
1077  * Write the field data to file. The file is named according to the session
1078  * name with the extension .fld appended.
1079  */
1080  std::string outname = m_sessionName + ".fld";
1081 
1082  WriteVessels(outname);
1083  }
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:345
Nektar::PulseWaveSystem::WriteVessels
void WriteVessels(const std::string &outname)
Write input fields to the given filename.
Definition: PulseWaveSystem.cpp:1103

◆ WriteVessels()

void Nektar::PulseWaveSystem::WriteVessels ( const std::string &  outname)
protected

Write input fields to the given filename.

Writes the field data to a file with the given filename.

Parameters
outnameFilename to write to.

Definition at line 1103 of file PulseWaveSystem.cpp.

References Nektar::SolverUtils::EquationSystem::m_fieldMetaDataMap, m_nDomains, m_nVariables, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_time, m_vessels, and Nektar::LibUtilities::Write().

Referenced by CheckPoint_Output(), and v_Output().

1104  {
1105 
1106  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1107  std::vector<std::string> variables = m_session->GetVariables();
1108 
1109  for(int n = 0; n < m_nDomains; ++n)
1110  {
1111  m_vessels[n*m_nVariables]->GetFieldDefinitions(FieldDef);
1112  }
1113 
1114  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
1115 
1116  int nFieldDefPerDomain = FieldDef.size()/m_nDomains;
1117  int cnt;
1118  // Copy Data into FieldData and set variable
1119  for(int n = 0; n < m_nDomains; ++n)
1120  {
1121  for(int j = 0; j < m_nVariables; ++j)
1122  {
1123  for(int i = 0; i < nFieldDefPerDomain; ++i)
1124  {
1125  cnt = n*nFieldDefPerDomain+i;
1126  FieldDef[cnt]->m_fields.push_back(variables[j]);
1127  m_vessels[n*m_nVariables]->AppendFieldData(FieldDef[cnt], FieldData[cnt], m_vessels[n*m_nVariables+j]->UpdateCoeffs());
1128  }
1129  }
1130  }
1131 
1132  // Update time in field info if required
1133  if(m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
1134  {
1135  m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
1136  }
1137 
1138  //LibUtilities::CombineFields(FieldDef, FieldData);
1139 
1140  LibUtilities::Write(outname, FieldDef, FieldData, m_fieldMetaDataMap);
1141  }
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:347
Nektar::SolverUtils::EquationSystem::m_fieldMetaDataMap
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
Definition: EquationSystem.h:393
Nektar::LibUtilities::Write
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap, const bool backup)
This function allows for data to be written to an FLD file when a session and/or communicator is not ...
Definition: FieldIO.cpp:249
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:333
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:104

Member Data Documentation

◆ m_A_0

Array<OneD, Array<OneD, NekDouble> > Nektar::PulseWaveSystem::m_A_0
protected

Definition at line 118 of file PulseWaveSystem.h.

Referenced by Nektar::PulseWavePropagation::GetFluxVector(), Nektar::PulseWavePropagation::SetPulseWaveBoundaryConditions(), and v_InitObject().

◆ m_A_0_trace

Array<OneD, Array<OneD, NekDouble> > Nektar::PulseWaveSystem::m_A_0_trace
protected

Definition at line 119 of file PulseWaveSystem.h.

Referenced by FillDataFromInterfacePoint(), Nektar::PulseWavePropagation::GetA0(), and v_InitObject().

◆ m_beta

Array<OneD, Array<OneD, NekDouble> > Nektar::PulseWaveSystem::m_beta
protected

Definition at line 120 of file PulseWaveSystem.h.

Referenced by CalcCharacteristicVariables(), Nektar::PulseWavePropagation::GetFluxVector(), Nektar::PulseWavePropagation::SetPulseWaveBoundaryConditions(), and v_InitObject().

◆ m_beta_trace

Array<OneD, Array<OneD, NekDouble> > Nektar::PulseWaveSystem::m_beta_trace
protected

Definition at line 121 of file PulseWaveSystem.h.

Referenced by FillDataFromInterfacePoint(), Nektar::PulseWavePropagation::GetBeta(), and v_InitObject().

◆ m_bifurcations

std::vector<std::vector<InterfacePointShPtr> > Nektar::PulseWaveSystem::m_bifurcations
protected

Definition at line 126 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_C

NekDouble Nektar::PulseWaveSystem::m_C
protected

Definition at line 114 of file PulseWaveSystem.h.

◆ m_currentDomain

int Nektar::PulseWaveSystem::m_currentDomain
protected

Definition at line 106 of file PulseWaveSystem.h.

Referenced by Nektar::PulseWavePropagation::DoOdeRhs(), Nektar::PulseWavePropagation::GetA0(), Nektar::PulseWavePropagation::GetBeta(), Nektar::PulseWavePropagation::GetFluxVector(), and Nektar::PulseWavePropagation::GetN().

◆ m_fieldPhysOffset

Array<OneD, int> Nektar::PulseWaveSystem::m_fieldPhysOffset
protected

Definition at line 110 of file PulseWaveSystem.h.

Referenced by FillDataFromInterfacePoint(), and v_InitObject().

◆ m_mergingJcts

std::vector<std::vector<InterfacePointShPtr> > Nektar::PulseWaveSystem::m_mergingJcts
protected

Definition at line 127 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_nDomains

int Nektar::PulseWaveSystem::m_nDomains
protected

Definition at line 105 of file PulseWaveSystem.h.

Referenced by Nektar::PulseWavePropagation::DoOdeRhs(), Nektar::PulseWavePropagation::SetPulseWaveBoundaryConditions(), SetUpDomainInterfaces(), v_DoInitialise(), v_DoSolve(), v_InitObject(), v_L2Error(), v_LinfError(), and WriteVessels().

◆ m_nVariables

int Nektar::PulseWaveSystem::m_nVariables
protected

Definition at line 107 of file PulseWaveSystem.h.

Referenced by Nektar::PulseWavePropagation::DoOdeProjection(), Nektar::PulseWavePropagation::DoOdeRhs(), EnforceInterfaceConditions(), FillDataFromInterfacePoint(), Nektar::PulseWavePropagation::GetFluxVector(), SetUpDomainInterfaces(), v_DoInitialise(), v_DoSolve(), v_InitObject(), v_L2Error(), v_LinfError(), and WriteVessels().

◆ m_pext

NekDouble Nektar::PulseWaveSystem::m_pext
protected

Definition at line 112 of file PulseWaveSystem.h.

Referenced by Nektar::PulseWavePropagation::GetFluxVector(), Nektar::PulseWavePropagation::GetPext(), and v_InitObject().

◆ m_pout

NekDouble Nektar::PulseWaveSystem::m_pout
protected

Definition at line 116 of file PulseWaveSystem.h.

◆ m_rho

NekDouble Nektar::PulseWaveSystem::m_rho
protected

Definition at line 111 of file PulseWaveSystem.h.

Referenced by BifurcationRiemann(), CalcCharacteristicVariables(), Nektar::PulseWavePropagation::GetFluxVector(), Nektar::PulseWavePropagation::GetRho(), JunctionRiemann(), MergingRiemann(), and v_InitObject().

◆ m_RT

NekDouble Nektar::PulseWaveSystem::m_RT
protected

Definition at line 115 of file PulseWaveSystem.h.

◆ m_trace_fwd_normal

Array<OneD, Array<OneD, NekDouble> > Nektar::PulseWaveSystem::m_trace_fwd_normal
protected

Definition at line 122 of file PulseWaveSystem.h.

Referenced by Nektar::PulseWavePropagation::GetN(), and v_InitObject().

◆ m_upwindTypePulse

UpwindTypePulse Nektar::PulseWaveSystem::m_upwindTypePulse
protected

Definition at line 108 of file PulseWaveSystem.h.

Referenced by Nektar::PulseWavePropagation::v_InitObject(), and v_InitObject().

◆ m_vesselJcts

std::vector<std::vector<InterfacePointShPtr> > Nektar::PulseWaveSystem::m_vesselJcts
protected

Definition at line 125 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_vessels

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::PulseWaveSystem::m_vessels
protected

Definition at line 104 of file PulseWaveSystem.h.

Referenced by CalcCharacteristicVariables(), Nektar::PulseWavePropagation::DoOdeRhs(), EnforceInterfaceConditions(), FillDataFromInterfacePoint(), Nektar::PulseWavePropagation::GetFluxVector(), Nektar::PulseWavePropagation::SetPulseWaveBoundaryConditions(), SetUpDomainInterfaces(), v_DoInitialise(), v_DoSolve(), Nektar::PulseWavePropagation::v_InitObject(), v_InitObject(), v_L2Error(), v_LinfError(), and WriteVessels().

Generated on Wed Jan 20 2021 03:33:55 for Nektar++ by   doxygen 1.8.13