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::ExpListSharedPtrUpdateVessels (void)
 
- 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...
 
SOLVER_UTILS_EXPORT void SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
- 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, NekDoubleErrorExtraPoints (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::FieldMetaDataMapUpdateFieldMetaDataMap ()
 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, Array< OneD, NekDouble > &alpha)
 Riemann Problem for Bifurcation. More...
 
void MergingRiemann (Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0, Array< OneD, NekDouble > &alpha)
 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, Array< OneD, NekDouble > &alpha)
 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 ()
 
virtual SOLVER_UTILS_EXPORT void v_SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
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...
 
virtual SOLVER_UTILS_EXPORT bool UpdateTimeStepCheck ()
 
- 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::ExpListSharedPtrm_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_gamma
 
Array< OneD, Array< OneD, NekDouble > > m_alpha
 
Array< OneD, Array< OneD, NekDouble > > m_alpha_trace
 
Array< OneD, Array< OneD, NekDouble > > m_trace_fwd_normal
 
Array< OneD, Array< OneD, NekDouble > > m_pressure
 
PulseWavePressureAreaSharedPtr m_pressureArea
 
bool extraFields = false
 
Array< OneD, Array< OneD, NekDouble > > m_PWV
 
Array< OneD, Array< OneD, NekDouble > > m_W1
 
Array< OneD, Array< OneD, NekDouble > > m_W2
 
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::TimeIntegrationSchemeSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
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...
 
NekDouble m_steadyStateRes = 1.0
 
NekDouble m_steadyStateRes0 = 1.0
 
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...
 
NekDouble m_TimeIntegLambda =0.0
 coefff of spacial derivatives(rhs or m_F in GLM) in calculating the residual of the whole equation(used in unsteady time integrations) More...
 
bool m_flagImplicitItsStatistics
 
bool m_flagImplicitSolver = false
 
Array< OneD, NekDoublem_magnitdEstimat
 estimate the magnitude of each conserved varibles More...
 
Array< OneD, NekDoublem_locTimeStep
 local time step(notice only for jfnk other see m_cflSafetyFactor) More...
 
NekDouble m_inArrayNorm =-1.0
 
int m_TotLinItePerStep =0
 
int m_StagesPerStep =1
 
bool m_flagUpdatePreconMat
 
int m_maxLinItePerNewton
 
int m_TotNewtonIts =0
 
int m_TotLinIts =0
 
int m_TotImpStages =0
 
bool m_CalcPhysicalAV = true
 flag to update artificial viscosity More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
bool m_verbose
 
bool m_root
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtrm_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_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_timestepMax = -1.0
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
NekDouble m_lastCheckTime
 
NekDouble m_TimeIncrementFactor
 
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, NekDouble &alpha)
 

Additional Inherited Members

- Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 
NekDouble m_cflNonAcoustic
 
NekDouble m_CFLGrowth
 CFL growth rate. More...
 
NekDouble m_CFLEnd
 maximun cfl in cfl growth 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 87 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 }
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.

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,
Array< OneD, NekDouble > &  alpha 
)
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 741 of file PulseWaveSystem.cpp.

746 {
747  NekDouble rho = m_rho;
748  Array<OneD, NekDouble> W(3);
749  Array<OneD, NekDouble> P_Au(3);
750  Array<OneD, NekDouble> W_Au(3);
751  NekMatrix<NekDouble> invJ(6, 6);
752  NekVector<NekDouble> f(6);
753  NekVector<NekDouble> dx(6);
754 
755  int proceed = 1;
756  int iter = 0;
757  int MAX_ITER = 100;
758 
759  // Forward characteristic
760  m_pressureArea->GetW1(W[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
761 
762  // Backward characteristics
763  for (int i = 1; i < 3; ++i)
764  {
765  m_pressureArea->GetW2(W[i], uu[i], beta[i], Au[i], A_0[i], alpha[i]);
766  }
767 
768  // Tolerances for the algorithm
769  NekDouble Tol = 1.0E-10;
770 
771  // Newton Iteration
772  while ((proceed) && (iter < MAX_ITER))
773  {
774  iter += 1;
775 
776  /*
777  * We solve the six constraint equations via a multivariate Newton
778  * iteration. Equations are:
779  * 1. Forward characteristic: W1(A_L, U_L) = W1(Au_L, Uu_L)
780  * 2. Backward characteristic 1: W2(A_R1, U_R1) = W2(Au_R1, Uu_R1)
781  * 3. Backward characteristic 2: W2(A_R2, U_R2) = W2(Au_R2, Uu_R2)
782  * 4. Conservation of mass: Au_L * Uu_L = Au_R1 * Uu_R1 +
783  * Au_R2 * Uu_R2
784  * 5. Continuity of total pressure 1: rho * Uu_L * Uu_L / 2 + p(Au_L) =
785  * rho * Uu_R1 * Uu_R1 / 2 + p(Au_R1)
786  * 6. Continuity of total pressure 2: rho * Uu_L * Uu_L / 2 + p(Au_L) =
787  * rho * Uu_R2 * Uu_R2 / 2 + p(Au_R2)
788  */
789  for (int i = 0; i < 3; ++i)
790  {
791  m_pressureArea->GetPressure(P_Au[i], beta[i], Au[i], A_0[i], 0, 0,
792  alpha[i]);
793  }
794 
795  m_pressureArea->GetW1(W_Au[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
796  m_pressureArea->GetW2(W_Au[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
797  m_pressureArea->GetW2(W_Au[2], uu[2], beta[2], Au[2], A_0[2], alpha[2]);
798 
799  // Constraint equations set to zero
800  f[0] = W_Au[0] - W[0];
801  f[1] = W_Au[1] - W[1];
802  f[2] = W_Au[2] - W[2];
803  f[3] = Au[0] * uu[0] - Au[1] * uu[1] - Au[2] * uu[2];
804  f[4] = uu[0] * uu[0] + 2 * P_Au[0] / rho -
805  uu[1] * uu[1] - 2 * P_Au[1] / rho;
806  f[5] = uu[0] * uu[0] + 2 * P_Au[0] / rho -
807  uu[2] * uu[2] - 2 * P_Au[2] / rho;
808 
809  // Inverse Jacobian for x + dx = x - J^(-1) * f(x)
810  m_pressureArea->GetJacobianInverse(invJ, Au, uu, beta, A_0, alpha,
811  "Bifurcation");
812 
813  Multiply(dx, invJ, f);
814 
815  // Update the solution: x_new = x_old - dx
816  for (int i = 0; i < 3; ++i)
817  {
818  uu[i] -= dx[i];
819  Au[i] -= dx[i + 3];
820  }
821 
822  // Check if the error of the solution is smaller than Tol
823  if (Dot(dx, dx) < Tol)
824  {
825  proceed = 0;
826  }
827 
828  // Check if solver converges
829  if (iter >= MAX_ITER)
830  {
831  ASSERTL0(false, "Riemann solver for Bifurcation did not converge.");
832  }
833  }
834 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
PulseWavePressureAreaSharedPtr m_pressureArea
void Multiply(NekMatrix< ResultDataType, StandardMatrixTag > &result, const NekMatrix< LhsDataType, LhsMatrixType > &lhs, const ResultDataType &rhs)
double NekDouble
DataType Dot(const NekVector< DataType > &lhs, const NekVector< DataType > &rhs)
Definition: NekVector.cpp:1220

References ASSERTL0, Nektar::Dot(), m_pressureArea, m_rho, and Nektar::Multiply().

Referenced by EnforceInterfaceConditions().

◆ 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 1046 of file PulseWaveSystem.cpp.

1047 {
1048  std::stringstream outname;
1049  outname << m_sessionName << "_" << n << ".chk";
1050 
1051  WriteVessels(outname.str());
1052 }
void WriteVessels(const std::string &outname)
Write input fields to the given filename.
std::string m_sessionName
Name of the session.

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

Referenced by v_DoSolve().

◆ EnforceInterfaceConditions()

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

Definition at line 644 of file PulseWaveSystem.cpp.

646 {
647  int dom, bcpos;
648  Array<OneD, NekDouble> Au(3);
649  Array<OneD, NekDouble> uu(3);
650  Array<OneD, NekDouble> beta(3);
651  Array<OneD, NekDouble> A_0(3);
652  Array<OneD, NekDouble> alpha(3);
653 
654  // Enforce Bifurcations:
655  for (int n = 0; n < m_bifurcations.size(); ++n)
656  {
657  for (int i = 0; i < 3; ++i)
658  {
659  FillDataFromInterfacePoint(m_bifurcations[n][i], fields, Au[i],
660  uu[i], beta[i], A_0[i], alpha[i]);
661  }
662  // Solve the Riemann problem for a bifurcation
663  BifurcationRiemann(Au, uu, beta, A_0, alpha);
664 
665  // Store the values into the right positions:
666  for (int i = 0; i < 3; ++i)
667  {
668  dom = m_bifurcations[n][i]->m_domain;
669  bcpos = m_bifurcations[n][i]->m_bcPosition;
670  m_vessels[dom * m_nVariables]
671  ->UpdateBndCondExpansion(bcpos)
672  ->UpdatePhys()[0] = Au[i];
673  m_vessels[dom * m_nVariables + 1]
674  ->UpdateBndCondExpansion(bcpos)
675  ->UpdatePhys()[0] = uu[i];
676  }
677  }
678 
679  // Enforce Bifurcations;
680  for (int n = 0; n < m_mergingJcts.size(); ++n)
681  {
682  // Merged vessel
683  for (int i = 0; i < 3; ++i)
684  {
685  FillDataFromInterfacePoint(m_mergingJcts[n][i], fields, Au[i],
686  uu[i], beta[i], A_0[i], alpha[i]);
687  }
688 
689  // Solve the Riemann problem for a merging vessel
690  MergingRiemann(Au, uu, beta, A_0, alpha);
691 
692  // Store the values into the right positions:
693  for (int i = 0; i < 3; ++i)
694  {
695  int dom = m_mergingJcts[n][i]->m_domain;
696  int bcpos = m_mergingJcts[n][i]->m_bcPosition;
697  m_vessels[dom * m_nVariables]
698  ->UpdateBndCondExpansion(bcpos)
699  ->UpdatePhys()[0] = Au[i];
700  m_vessels[dom * m_nVariables + 1]
701  ->UpdateBndCondExpansion(bcpos)
702  ->UpdatePhys()[0] = uu[i];
703  }
704  }
705 
706  for (int n = 0; n < m_vesselJcts.size(); ++n)
707  {
708  for (int i = 0; i < 2; ++i)
709  {
710  FillDataFromInterfacePoint(m_vesselJcts[n][i], fields, Au[i], uu[i],
711  beta[i], A_0[i], alpha[i]);
712  }
713 
714  JunctionRiemann(Au, uu, beta, A_0, alpha);
715 
716  // Store the values into the right positions:
717  for (int i = 0; i < 2; ++i)
718  {
719  int dom = m_vesselJcts[n][i]->m_domain;
720  int bcpos = m_vesselJcts[n][i]->m_bcPosition;
721  m_vessels[dom * m_nVariables]
722  ->UpdateBndCondExpansion(bcpos)
723  ->UpdatePhys()[0] = Au[i];
724  m_vessels[dom * m_nVariables + 1]
725  ->UpdateBndCondExpansion(bcpos)
726  ->UpdatePhys()[0] = uu[i];
727  }
728  }
729 }
void JunctionRiemann(Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0, Array< OneD, NekDouble > &alpha)
Riemann Problem for Junction.
void MergingRiemann(Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0, Array< OneD, NekDouble > &alpha)
Riemann Problem for Merging Flow.
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
void FillDataFromInterfacePoint(InterfacePointShPtr &I, const Array< OneD, const Array< OneD, NekDouble > > &field, NekDouble &A, NekDouble &u, NekDouble &beta, NekDouble &A_0, NekDouble &alpha)
std::vector< std::vector< InterfacePointShPtr > > m_bifurcations
std::vector< std::vector< InterfacePointShPtr > > m_vesselJcts
std::vector< std::vector< InterfacePointShPtr > > m_mergingJcts
void BifurcationRiemann(Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0, Array< OneD, NekDouble > &alpha)
Riemann Problem for Bifurcation.

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

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

◆ FillDataFromInterfacePoint()

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

Definition at line 618 of file PulseWaveSystem.cpp.

622 {
623  int omega = I->m_domain;
624  int traceId = I->m_traceId;
625  int eid = I->m_elmt;
626  int vert = I->m_elmtVert;
627  int vesselID = omega * m_nVariables;
628  int phys_offset = m_vessels[vesselID]->GetPhys_Offset(eid);
630  Array<OneD, NekDouble> Tmp(1);
631 
632  m_vessels[vesselID]->GetExp(eid)->GetTracePhysVals(
633  vert, dumExp, fields[0] + m_fieldPhysOffset[omega] + phys_offset, Tmp);
634  A = Tmp[0];
635  m_vessels[vesselID]->GetExp(eid)->GetTracePhysVals(
636  vert, dumExp, fields[1] + m_fieldPhysOffset[omega] + phys_offset, Tmp);
637  u = Tmp[0];
638 
639  beta = m_beta_trace[omega][traceId];
640  A_0 = m_A_0_trace[omega][traceId];
641  alpha = m_alpha_trace[omega][traceId];
642 }
Array< OneD, int > m_fieldPhysOffset
Array< OneD, Array< OneD, NekDouble > > m_beta_trace
Array< OneD, Array< OneD, NekDouble > > m_alpha_trace
Array< OneD, Array< OneD, NekDouble > > m_A_0_trace
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68

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

Referenced by EnforceInterfaceConditions().

◆ GetNdomains()

int Nektar::PulseWaveSystem::GetNdomains ( )
inline

Definition at line 93 of file PulseWaveSystem.h.

94  {
95  return m_nDomains;
96  }

◆ JunctionRiemann()

void Nektar::PulseWaveSystem::JunctionRiemann ( Array< OneD, NekDouble > &  Au,
Array< OneD, NekDouble > &  uu,
Array< OneD, NekDouble > &  beta,
Array< OneD, NekDouble > &  A_0,
Array< OneD, NekDouble > &  alpha 
)
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 949 of file PulseWaveSystem.cpp.

954 {
955  NekDouble rho = m_rho;
956  Array<OneD, NekDouble> W(2);
957  Array<OneD, NekDouble> W_Au(2);
958  Array<OneD, NekDouble> P_Au(2);
959  NekMatrix<NekDouble> invJ(4, 4);
960  NekVector<NekDouble> f(4);
961  NekVector<NekDouble> dx(4);
962 
963  int proceed = 1;
964  int iter = 0;
965  int MAX_ITER = 15;
966  NekDouble Tol = 1.0E-10;
967 
968  // Forward and backward characteristics
969  m_pressureArea->GetW1(W[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
970  m_pressureArea->GetW2(W[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
971 
972  while ((proceed) && (iter < MAX_ITER))
973  {
974  iter += 1;
975 
976  /*
977  * We solve the four constraint equations via a multivariate Newton
978  * iteration. Equations are:
979  * 1. Forward characteristic: W1(A_L, U_L) = W1(Au_L, Uu_L)
980  * 2. Backward characteristic: W2(A_R, U_R) = W2(Au_R, Uu_R)
981  * 3. Conservation of mass: Au_L * Uu_L = Au_R * Uu_R
982  * 4. Continuity of total pressure: rho * Uu_L * Uu_L / 2 + p(Au_L) =
983  * rho * Uu_R * Uu_R / 2 + p(Au_R)
984  */
985  for (int i = 0; i < 2; ++i)
986  {
987  m_pressureArea->GetPressure(P_Au[i], beta[i], Au[i], A_0[i], 0, 0,
988  alpha[i]);
989  }
990 
991  m_pressureArea->GetW1(W_Au[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
992  m_pressureArea->GetW2(W_Au[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
993 
994  // Constraint equations set to zero
995  f[0] = W_Au[0] - W[0];
996  f[1] = W_Au[1] - W[1];
997  f[2] = Au[0] * uu[0] - Au[1] * uu[1];
998  f[3] = uu[0] * uu[0] + 2 * P_Au[0] / rho -
999  uu[1] * uu[1] - 2 * P_Au[1] / rho;
1000 
1001  // Inverse Jacobian for x + dx = x - J^(-1) * f(x)
1002  m_pressureArea->GetJacobianInverse(invJ, Au, uu, beta, A_0, alpha,
1003  "Junction");
1004 
1005  Multiply(dx, invJ, f);
1006 
1007  // Update solution: x_new = x_old - dx
1008  for (int i = 0; i < 2; ++i)
1009  {
1010  uu[i] -= dx[i];
1011  Au[i] -= dx[i + 2];
1012  }
1013 
1014  // Check if the error of the solution is smaller than Tol.
1015  if (Dot(dx, dx) < Tol)
1016  {
1017  proceed = 0;
1018  }
1019  }
1020 
1021  if (iter >= MAX_ITER)
1022  {
1023  ASSERTL0(false, "Riemann solver for Junction did not converge");
1024  }
1025 }

References ASSERTL0, Nektar::Dot(), m_pressureArea, m_rho, and Nektar::Multiply().

Referenced by EnforceInterfaceConditions().

◆ 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,
Array< OneD, NekDouble > &  alpha 
)
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 845 of file PulseWaveSystem.cpp.

850 {
851  NekDouble rho = m_rho;
852  Array<OneD, NekDouble> W(3);
853  Array<OneD, NekDouble> W_Au(3);
854  Array<OneD, NekDouble> P_Au(3);
855  NekMatrix<NekDouble> invJ(6, 6);
856  NekVector<NekDouble> f(6);
857  NekVector<NekDouble> dx(6);
858 
859  int proceed = 1;
860  int iter = 0;
861  int MAX_ITER = 15;
862 
863  // Backward characteristic
864  m_pressureArea->GetW2(W[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
865 
866  // Forward characteristics
867  for (int i = 1; i < 3; ++i)
868  {
869  m_pressureArea->GetW1(W[i], uu[i], beta[i], Au[i], A_0[i], alpha[i]);
870  }
871 
872  // Tolerances for the algorithm
873  NekDouble Tol = 1.0E-10;
874 
875  // Newton Iteration
876  while ((proceed) && (iter < MAX_ITER))
877  {
878  iter += 1;
879 
880  /*
881  * We solve the six constraint equations via a multivariate Newton
882  * iteration. Equations are:
883  * 1. Backward characteristic: W2(A_R, U_R) = W1(Au_R, Uu_R)
884  * 2. Forward characteristic 1: W1(A_L1, U_L1) = W1(Au_L1, Uu_R1)
885  * 3. Forward characteristic 2: W1(A_L2, U_L2) = W1(Au_L2, Uu_L2)
886  * 4. Conservation of mass: Au_R * Uu_R = Au_L1 * Uu_L1 +
887  * Au_L2 * Uu_L2
888  * 5. Continuity of total pressure 1: rho * Uu_R * Uu_R / 2 + p(Au_R) =
889  * rho * Uu_L1 * Uu_L1 / 2 + p(Au_L1)
890  * 6. Continuity of total pressure 2: rho * Uu_R * Uu_R / 2 + p(Au_R) =
891  * rho * Uu_L2 * Uu_L2 / 2 + p(Au_L2)
892  */
893  for (int i = 0; i < 3; ++i)
894  {
895  m_pressureArea->GetPressure(P_Au[i], beta[i], Au[i], A_0[i], 0, 0,
896  alpha[i]);
897  }
898 
899  m_pressureArea->GetW2(W_Au[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
900  m_pressureArea->GetW1(W_Au[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
901  m_pressureArea->GetW1(W_Au[2], uu[2], beta[2], Au[2], A_0[2], alpha[2]);
902 
903  // Constraint equations set to zero
904  f[0] = W_Au[0] - W[0];
905  f[1] = W_Au[1] - W[1];
906  f[2] = W_Au[2] - W[2];
907  f[3] = Au[0] * uu[0] - Au[1] * uu[1] - Au[2] * uu[2];
908  f[4] = uu[0] * uu[0] + 2 * P_Au[0] / rho -
909  uu[1] * uu[1] - 2 * P_Au[1] / rho;
910  f[5] = uu[0] * uu[0] + 2 * P_Au[0] / rho -
911  uu[2] * uu[2] - 2 * P_Au[2] / rho;
912 
913  // Inverse Jacobian for x + dx = x - J^(-1) * f(x)
914  m_pressureArea->GetJacobianInverse(invJ, Au, uu, beta, A_0, alpha,
915  "Merge");
916 
917  Multiply(dx, invJ, f);
918 
919  // Update the solution: x_new = x_old - dx
920  for (int i = 0; i < 3; ++i)
921  {
922  uu[i] -= dx[i];
923  Au[i] -= dx[i + 3];
924  }
925 
926  // Check if the error of the solution is smaller than Tol
927  if (Dot(dx, dx) < Tol)
928  {
929  proceed = 0;
930  }
931 
932  // Check if solver converges
933  if (iter >= MAX_ITER)
934  {
935  ASSERTL0(false, "Riemann solver for Merging Flow did not converge");
936  }
937  }
938 }

References ASSERTL0, Nektar::Dot(), m_pressureArea, m_rho, and Nektar::Multiply().

Referenced by EnforceInterfaceConditions().

◆ SetUpDomainInterfaces()

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

Definition at line 320 of file PulseWaveSystem.cpp.

321 {
322  map<int, std::vector<InterfacePointShPtr> > VidToDomain;
323 
324  /* Loop over domain and find out if we have any undefined
325  * boundary conditions representing interfaces. If so make a
326  * map based around vid and storing the domains that are
327  * part of interfaces.
328  */
329  for (int omega = 0; omega < m_nDomains; ++omega)
330  {
331  int vesselID = omega * m_nVariables;
332 
333  for (int i = 0; i < 2; ++i)
334  {
335  if (m_vessels[vesselID]
336  ->GetBndConditions()[i]
337  ->GetBoundaryConditionType() == SpatialDomains::eNotDefined)
338  {
339  // Get Vid of interface
340  int vid = m_vessels[vesselID]
341  ->UpdateBndCondExpansion(i)
342  ->GetExp(0)
343  ->GetGeom()
344  ->GetVid(0);
345  cout << "Shared vertex ID: " << vid << endl;
346  MultiRegions::ExpListSharedPtr trace = m_vessels[vesselID]->GetTrace();
348 
349  bool finish = false;
350  // find which elmt, the local vertex and the data offset of point
351  for (int n = 0; n < m_vessels[vesselID]->GetExpSize(); ++n)
352  {
353  for (int p = 0; p < 2; ++p)
354  {
355  if (m_vessels[vesselID]
356  ->GetTraceMap()
357  ->GetElmtToTrace()[n][p]
358  ->as<LocalRegions::Expansion>()
359  ->GetGeom()
360  ->GetVid(0) == vid)
361  {
362  int eid = m_vessels[vesselID]
363  ->GetTraceMap()
364  ->GetElmtToTrace()[n][p]
365  ->GetElmtId();
366 
367  int tid =
368  m_vessels[vesselID]->GetTrace()->GetCoeff_Offset(
369  eid);
370  Ipt =
372  vid, omega, n, p, tid, i);
373 
374  cout << "Global VID of interface point: " << vid << endl;
375  cout << "Domain interface point belongs to: " << omega << endl;
376  cout << "Element ID of vertex: " << n << endl;
377  cout << "Vertex ID within local element: " << p << endl;
378  cout << "Element ID within the trace: " << tid << endl;
379  cout << "Position of boundary condition in region: " << i << endl;
380 
381  finish = true;
382  break;
383  }
384  }
385  if (finish == true)
386  {
387  break;
388  }
389  }
390 
391  VidToDomain[vid].push_back(Ipt);
392 
393  // finally reset boundary condition to Dirichlet
394  m_vessels[vesselID]->GetBndConditions()[i]->SetBoundaryConditionType(
396 
397  m_vessels[vesselID + 1]
398  ->GetBndConditions()[i]
399  ->SetBoundaryConditionType(SpatialDomains::eDirichlet);
400  }
401  }
402  }
403 
404  // loop over map and set up Interface information;
405  for (auto &iter : VidToDomain)
406  {
407  if (iter.second.size() == 2) // Vessel jump interface
408  {
409  m_vesselJcts.push_back(iter.second);
410  }
411  else if (iter.second.size() == 3) // Bifurcation or Merging junction.
412  {
413  int nbeg = 0;
414  int nend = 0;
415 
416  // Determine if bifurcation or merging junction
417  // through number of element vertices that meet at
418  // junction. Only one vertex using a m_elmtVert=1
419  // indicates a bifurcation
420  for (int i = 0; i < 3; ++i)
421  {
422  if (iter.second[i]->m_elmtVert == 0)
423  {
424  nbeg += 1;
425  }
426  else
427  {
428  nend += 1;
429  }
430  }
431 
432  // Set up Bifurcation information
433  if (nbeg == 2)
434  {
435  // ensure first InterfacePoint is parent
436  if (iter.second[0]->m_elmtVert ==
437  1) // m_elmtVert: Vertex id in local element
438  {
439  m_bifurcations.push_back(iter.second);
440  }
441  else
442  {
443  // order points according to Riemann solver convention
445  // find merging vessel
446  if (iter.second[1]->m_elmtVert == 1)
447  {
448  I = iter.second[0];
449  iter.second[0] = iter.second[1];
450  iter.second[1] = I;
451  }
452  else if (iter.second[2]->m_elmtVert == 1)
453  {
454  I = iter.second[0];
455  iter.second[0] = iter.second[2];
456  iter.second[2] = I;
457  }
459  "This routine has not been checked");
460  }
461  }
462  else
463  {
464  // ensure last InterfacePoint is merged vessel
465  if (iter.second[0]->m_elmtVert == 0)
466  {
467  m_mergingJcts.push_back(iter.second);
468  }
469  else
470  {
471  // order points according to Riemann solver convention
473  // find merging vessel
474  if (iter.second[1]->m_elmtVert == 0)
475  {
476  I = iter.second[0];
477  iter.second[0] = iter.second[1];
478  iter.second[1] = I;
479  }
480  else if (iter.second[2]->m_elmtVert == 0)
481  {
482  I = iter.second[0];
483  iter.second[0] = iter.second[2];
484  iter.second[2] = I;
485  }
487  "This routine has not been checked");
488  }
489  }
490 
491  }
492  else
493  {
494  ASSERTL0(false, "Unknown junction type");
495  }
496  }
497 }
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< InterfacePoint > InterfacePointShPtr

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().

◆ UpdateVessels()

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

Definition at line 98 of file PulseWaveSystem.h.

99  {
100  return m_vessels;
101  }

◆ 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 506 of file PulseWaveSystem.cpp.

507 {
508  if (m_session->GetComm()->GetRank() == 0)
509  {
510  cout << "Initial Conditions: " << endl;
511  }
512 
513  /* Loop over all subdomains to initialize all with the Initial
514  * Conditions read from the inputfile*/
515  for (int omega = 0; omega < m_nDomains; omega++)
516  {
517  for (int i = 0; i < 2; ++i)
518  {
519  m_fields[i] = m_vessels[m_nVariables * omega + i];
520  }
521 
522  if (m_session->GetComm()->GetRank() == 0)
523  {
524  cout << "Subdomain: " << omega << endl;
525  }
526 
527  SetInitialConditions(0.0, 0, omega);
528  }
529 
530  // Reset to first definition
531  for (int i = 0; i < 2; ++i)
532  {
533  m_fields[i] = m_vessels[i];
534  }
535 }
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.

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

◆ 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 independently 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 554 of file PulseWaveSystem.cpp.

555 {
556  NekDouble IntegrationTime = 0.0;
557  int i;
558  int n;
559  int nchk = 1;
560 
561  Array<OneD, Array<OneD, NekDouble> > fields(m_nVariables);
562 
563  for (int i = 0; i < m_nVariables; ++i)
564  {
565  fields[i] = m_vessels[i]->UpdatePhys();
566  m_fields[i]->SetPhysState(false);
567  }
568 
569  m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
570 
571  // Time loop
572  for (n = 0; n < m_steps; ++n)
573  {
574  LibUtilities::Timer timer;
575  timer.Start();
576  fields = m_intScheme->TimeIntegrate(n, m_timestep, m_ode);
577  m_time += m_timestep;
578  timer.Stop();
579  IntegrationTime += timer.TimePerTest(1);
580 
581  // Write out status information.
582  if (m_session->GetComm()->GetRank() == 0 && !((n + 1) % m_infosteps))
583  {
584  cout << "Steps: " << n + 1 << "\t Time: " << m_time
585  << "\t Time-step: " << m_timestep << "\t" << endl;
586  }
587 
588  // Transform data if needed
589  if (!((n + 1) % m_checksteps))
590  {
591  for (i = 0; i < m_nVariables; ++i)
592  {
593  int cnt = 0;
594  for (int omega = 0; omega < m_nDomains; omega++)
595  {
596  m_vessels[omega * m_nVariables + i]->FwdTrans(
597  fields[i] + cnt,
598  m_vessels[omega * m_nVariables + i]->UpdateCoeffs());
599  cnt += m_vessels[omega * m_nVariables + i]->GetTotPoints();
600  }
601  }
602  CheckPoint_Output(nchk++);
603  }
604 
605  } // end of timeintegration
606 
607  // Copy Array To Vessel Phys Fields
608  for (int i = 0; i < m_nVariables; ++i)
609  {
610  Vmath::Vcopy(fields[i].size(), fields[i], 1,
611  m_vessels[i]->UpdatePhys(), 1);
612  }
613 
614  cout << "Time-integration timing: " << IntegrationTime << " s" << endl
615  << endl;
616 }
void CheckPoint_Output(const int n)
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
int m_steps
Number of steps to take.
int m_checksteps
Number of steps between checkpoints.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
int m_infosteps
Number of time steps between outputting status information.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199

References CheckPoint_Output(), Nektar::SolverUtils::EquationSystem::m_checksteps, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::UnsteadySystem::m_infosteps, Nektar::SolverUtils::UnsteadySystem::m_intScheme, 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().

◆ 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

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

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

Definition at line 87 of file PulseWaveSystem.cpp.

88 {
89  // Initialise base class
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"),
97  "Pulse solver only set up for Discontinuous projections");
99  ASSERTL0(m_graph->GetMeshDimension() == 1,
100  "Pulse wave solver only set up for expansion dimension equal to 1");
101 
102  int i;
103  m_nVariables = m_session->GetVariables().size();
104 
105  m_fields = Array<OneD, MultiRegions::ExpListSharedPtr>(m_nVariables);
106  m_vessels =
107  Array<OneD, MultiRegions::ExpListSharedPtr>(m_nVariables * m_nDomains);
108 
109  /* If the PressureArea property is not specified, default to the Beta law;
110  * it's the most well known and this way previous examples that did not
111  * specify the tube law will still work.
112  */
113  if (m_session->DefinesSolverInfo("PressureArea"))
114  {
116  m_session->GetSolverInfo("PressureArea"), m_vessels, m_session);
117  }
118  else
119  {
122  }
123 
124  m_pressure = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
125 
126  m_PWV = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
127  m_W1 = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
128  m_W2 = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
129 
130  const std::vector<SpatialDomains::CompositeMap> domain = m_graph->GetDomain();
131 
132  SpatialDomains::BoundaryConditions Allbcs(m_session, m_graph);
133 
134  // Set up domains and put geometry to be only one space dimension.
135  int cnt = 0;
136  bool SetToOneSpaceDimension = true;
137 
138  if (m_session->DefinesCmdLineArgument("SetToOneSpaceDimension"))
139  {
140  std::string cmdline =
141  m_session->GetCmdLineArgument<std::string>("SetToOneSpaceDimension");
142  if (boost::to_upper_copy(cmdline) == "FALSE")
143  {
144  SetToOneSpaceDimension = false;
145  }
146  }
147 
148  for (i = 0; i < m_nDomains; ++i)
149  {
150  for (int j = 0; j < m_nVariables; ++j)
151  {
152  m_vessels[cnt++] =
154  m_session, m_graph, domain[i], Allbcs,
155  m_session->GetVariable(j), SetToOneSpaceDimension);
156  }
157  }
158 
159  // Reset coeff and phys space to be continuous over all domains
160  int totcoeffs = 0;
161  int totphys = 0;
162  m_fieldPhysOffset = Array<OneD, int>(m_nDomains + 1, 0);
163  for (i = 0; i < m_nDomains; ++i)
164  {
165  totcoeffs += m_vessels[i * m_nVariables]->GetNcoeffs();
166 
167  m_fieldPhysOffset[i] = totphys;
168  totphys += m_vessels[i * m_nVariables]->GetTotPoints();
169  }
170  m_fieldPhysOffset[m_nDomains] = totphys;
171 
172  for (int n = 0; n < m_nVariables; ++n)
173  {
174  Array<OneD, NekDouble> coeffs(totcoeffs, 0.0);
175  Array<OneD, NekDouble> phys(totphys, 0.0);
176  Array<OneD, NekDouble> tmpcoeffs, tmpphys;
177 
178  m_vessels[n]->SetCoeffsArray(coeffs);
179  m_vessels[n]->SetPhysArray(phys);
180 
181  int cnt = m_vessels[n]->GetNcoeffs();
182  int cnt1 = m_vessels[n]->GetTotPoints();
183 
184  for (i = 1; i < m_nDomains; ++i)
185  {
186  m_vessels[i * m_nVariables + n]->SetCoeffsArray(tmpcoeffs = coeffs + cnt);
187  m_vessels[i * m_nVariables + n]->SetPhysArray(tmpphys = phys + cnt1);
188  cnt += m_vessels[i * m_nVariables + n]->GetNcoeffs();
189  cnt1 += m_vessels[i * m_nVariables + n]->GetTotPoints();
190  }
191  }
192 
193  for (int i = 0; i < 2; ++i)
194  {
195  m_fields[i] = m_vessels[i];
196  }
197 
198  // Zero all physical fields initially.
199  ZeroPhysFields();
200 
201  // If Discontinuous Galerkin determine upwinding method to use
202  for (int i = 0; i < (int)SIZE_UpwindTypePulse; ++i)
203  {
204  bool match;
205  m_session->MatchSolverInfo("UPWINDTYPEPULSE", UpwindTypeMapPulse[i],
206  match, false);
207  if (match)
208  {
210  break;
211  }
212  }
213 
214  // Load blood density and external pressure
215  m_session->LoadParameter("rho", m_rho, 0.5);
216  m_session->LoadParameter("pext", m_pext, 0.0);
217 
218  int nq = 0;
219  /*
220  * Gets the Material Properties of each arterial segment
221  * specified in the inputfile from section MaterialProperties
222  * Also gets the Area at static equilibrium A_0 specified in the
223  * inputfile.
224  *
225  * Having found these points also extract the values at the
226  * trace points and the normal direction consistent with the
227  * left adjacent definition of Fwd and Bwd
228  */
229  m_beta = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
230  m_beta_trace = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
231  m_gamma = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
232  m_alpha = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
233  m_alpha_trace = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
234  m_A_0 = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
235  m_A_0_trace = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
236  m_trace_fwd_normal = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
237 
238  for (int omega = 0; omega < m_nDomains; omega++)
239  {
240  nq = m_vessels[2 * omega]->GetNpoints();
241  m_fields[0] = m_vessels[2 * omega];
242 
243  m_beta[omega] = Array<OneD, NekDouble>(nq);
244  GetFunction("MaterialProperties")
245  ->Evaluate("beta", m_beta[omega], m_time, omega);
246 
247  // If the input file doesn't specify viscoelasticity, set it to zero.
248  m_gamma[omega] = Array<OneD, NekDouble>(nq);
249 
250  if (m_session->DefinesFunction("Viscoelasticity"))
251  {
252  GetFunction("Viscoelasticity")
253  ->Evaluate("gamma", m_gamma[omega], m_time, omega);
254  }
255  else
256  {
257  for (int j = 0; j < nq; ++j)
258  {
259  m_gamma[omega][j] = 0;
260  }
261  }
262 
263  /* If the input file doesn't specify strain-stiffening, set it to
264  * 0.5 for elastic behaviour.
265  */
266  m_alpha[omega] = Array<OneD, NekDouble>(nq);
267 
268  if (m_session->DefinesFunction("StrainStiffening"))
269  {
270  GetFunction("StrainStiffening")
271  ->Evaluate("alpha", m_alpha[omega], m_time, omega);
272  }
273  else
274  {
275  for (int j = 0; j < nq; ++j)
276  {
277  m_alpha[omega][j] = 0.5;
278  }
279  }
280 
281  m_A_0[omega] = Array<OneD, NekDouble>(nq);
282  GetFunction("A_0")->Evaluate("A_0", m_A_0[omega], m_time, omega);
283 
284  int nqTrace = GetTraceTotPoints();
285 
286  m_beta_trace[omega] = Array<OneD, NekDouble>(nqTrace);
287  m_fields[0]->ExtractTracePhys(m_beta[omega], m_beta_trace[omega]);
288 
289  m_A_0_trace[omega] = Array<OneD, NekDouble>(nqTrace);
290  m_fields[0]->ExtractTracePhys(m_A_0[omega], m_A_0_trace[omega]);
291 
292  m_alpha_trace[omega] = Array<OneD, NekDouble>(nqTrace);
293  m_fields[0]->ExtractTracePhys(m_alpha[omega], m_alpha_trace[omega]);
294 
295  if (SetToOneSpaceDimension)
296  {
297  m_trace_fwd_normal[omega] = Array<OneD, NekDouble>(nqTrace, 0.0);
298 
299  MultiRegions::ExpListSharedPtr trace = m_fields[0]->GetTrace();
300  int nelmt_trace = trace->GetExpSize();
301 
302  Array<OneD, Array<OneD, NekDouble> > normals(nelmt_trace);
303 
304  for (int i = 0; i < nelmt_trace; ++i)
305  {
306  normals[i] = m_trace_fwd_normal[omega] + i;
307  }
308 
309  // need to set to 1 for consistency since boundary
310  // conditions may not have coordim = 1
311  trace->GetExp(0)->GetGeom()->SetCoordim(1);
312 
313  trace->GetNormals(normals);
314  }
315  }
316 
318 }
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:145
Array< OneD, Array< OneD, NekDouble > > m_A_0
Array< OneD, Array< OneD, NekDouble > > m_W2
Array< OneD, Array< OneD, NekDouble > > m_trace_fwd_normal
Array< OneD, Array< OneD, NekDouble > > m_W1
UpwindTypePulse m_upwindTypePulse
Array< OneD, Array< OneD, NekDouble > > m_pressure
Array< OneD, Array< OneD, NekDouble > > m_gamma
Array< OneD, Array< OneD, NekDouble > > m_alpha
Array< OneD, Array< OneD, NekDouble > > m_PWV
Array< OneD, Array< OneD, NekDouble > > m_beta
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SOLVER_UTILS_EXPORT void ZeroPhysFields()
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
PressureAreaFactory & GetPressureAreaFactory()
@ SIZE_UpwindTypePulse
Length of enum list.
const char *const UpwindTypeMapPulse[]

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::MultiRegions::eDiscontinuous, Nektar::SolverUtils::EquationSystem::GetFunction(), Nektar::GetPressureAreaFactory(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), m_A_0, m_A_0_trace, m_alpha, m_alpha_trace, m_beta, m_beta_trace, m_fieldPhysOffset, Nektar::SolverUtils::EquationSystem::m_fields, m_gamma, Nektar::SolverUtils::EquationSystem::m_graph, m_nDomains, m_nVariables, m_pext, m_pressure, m_pressureArea, Nektar::SolverUtils::EquationSystem::m_projectionType, m_PWV, m_rho, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_time, m_trace_fwd_normal, m_upwindTypePulse, m_vessels, m_W1, m_W2, SetUpDomainInterfaces(), Nektar::SIZE_UpwindTypePulse, Nektar::UpwindTypeMapPulse, Nektar::SolverUtils::UnsteadySystem::v_InitObject(), and Nektar::SolverUtils::EquationSystem::ZeroPhysFields().

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

◆ 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 1152 of file PulseWaveSystem.cpp.

1155 {
1156  NekDouble L2error = 0.0;
1157  NekDouble L2error_dom;
1158  NekDouble Vol = 0.0;
1159 
1160  if (m_NumQuadPointsError == 0)
1161  {
1162  for (int omega = 0; omega < m_nDomains; omega++)
1163  {
1164  int vesselid = field + omega * m_nVariables;
1165 
1166  if (m_vessels[vesselid]->GetPhysState() == false)
1167  {
1168  m_vessels[vesselid]->BwdTrans(m_vessels[vesselid]->GetCoeffs(),
1169  m_vessels[vesselid]->UpdatePhys());
1170  }
1171 
1172  if (exactsoln.size())
1173  {
1174  L2error_dom = m_vessels[vesselid]->L2(
1175  m_vessels[vesselid]->GetPhys(), exactsoln);
1176  }
1177  else if (m_session->DefinesFunction("ExactSolution"))
1178  {
1179  Array<OneD, NekDouble> exactsoln(
1180  m_vessels[vesselid]->GetNpoints());
1181 
1183  m_session->GetFunction("ExactSolution", field, omega);
1184  GetFunction("ExactSolution")
1185  ->Evaluate(m_session->GetVariable(field), exactsoln,
1186  m_time);
1187 
1188  L2error_dom = m_vessels[vesselid]->L2(
1189  m_vessels[vesselid]->GetPhys(), exactsoln);
1190  }
1191  else
1192  {
1193  L2error_dom =
1194  m_vessels[vesselid]->L2(m_vessels[vesselid]->GetPhys());
1195  }
1196 
1197  L2error += L2error_dom * L2error_dom;
1198 
1199  if (Normalised == true)
1200  {
1201  Array<OneD, NekDouble> one(m_vessels[vesselid]->GetNpoints(),
1202  1.0);
1203 
1204  Vol += m_vessels[vesselid]->Integral(one);
1205  }
1206  }
1207  }
1208  else
1209  {
1210  ASSERTL0(false, "Not set up");
1211  }
1212 
1213  if (Normalised == true)
1214  {
1215  m_comm->AllReduce(Vol, LibUtilities::ReduceSum);
1216 
1217  L2error = sqrt(L2error / Vol);
1218  }
1219  else
1220  {
1221  L2error = sqrt(L2error);
1222  }
1223 
1224  return L2error;
1225 }
LibUtilities::CommSharedPtr m_comm
Communicator.
SOLVER_UTILS_EXPORT int GetNpoints()
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:131
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267

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, Nektar::LibUtilities::ReduceSum, and tinysimd::sqrt().

◆ 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 1233 of file PulseWaveSystem.cpp.

1235 {
1236  NekDouble LinferrorDom, Linferror = -1.0;
1237 
1238  for (int omega = 0; omega < m_nDomains; ++omega)
1239  {
1240  int vesselid = field + omega * m_nVariables;
1241 
1242  if (m_NumQuadPointsError == 0)
1243  {
1244  if (m_vessels[vesselid]->GetPhysState() == false)
1245  {
1246  m_vessels[vesselid]->BwdTrans(m_vessels[vesselid]->GetCoeffs(),
1247  m_vessels[vesselid]->UpdatePhys());
1248  }
1249 
1250  if (exactsoln.size())
1251  {
1252  LinferrorDom =
1253  m_vessels[vesselid]->Linf(m_vessels[vesselid]->GetPhys(),
1254  exactsoln);
1255  }
1256  else if (m_session->DefinesFunction("ExactSolution"))
1257  {
1258  Array<OneD, NekDouble>
1259  exactsoln(m_vessels[vesselid]->GetNpoints());
1260 
1261  GetFunction("ExactSolution")
1262  ->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
1263 
1264  LinferrorDom =
1265  m_vessels[vesselid]->Linf(m_vessels[vesselid]->GetPhys(),
1266  exactsoln);
1267  }
1268  else
1269  {
1270  LinferrorDom = 0.0;
1271  }
1272 
1273  Linferror = (Linferror > LinferrorDom) ? Linferror : LinferrorDom;
1274  }
1275  else
1276  {
1277  ASSERTL0(false, "ErrorExtraPoints not allowed for this solver");
1278  }
1279  }
1280  return Linferror;
1281 }

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.

◆ 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 1031 of file PulseWaveSystem.cpp.

1032 {
1033  /**
1034  * Write the field data to file. The file is named according to the session
1035  * name with the extension .fld appended.
1036  */
1037  std::string outname = m_sessionName + ".fld";
1038 
1039  WriteVessels(outname);
1040 }

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

◆ 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 1058 of file PulseWaveSystem.cpp.

1059 {
1060  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1061  std::vector<std::string> variables = m_session->GetVariables();
1062 
1063  for (int n = 0; n < m_nDomains; ++n)
1064  {
1065  m_vessels[n * m_nVariables]->GetFieldDefinitions(FieldDef);
1066  }
1067 
1068  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
1069 
1070  int nFieldDefPerDomain = FieldDef.size() / m_nDomains;
1071  int cnt;
1072 
1073  // Copy Data into FieldData and set variable
1074  for (int n = 0; n < m_nDomains; ++n)
1075  {
1076  // Outputs area and velocity
1077  for (int j = 0; j < m_nVariables; ++j)
1078  {
1079  for (int i = 0; i < nFieldDefPerDomain; ++i)
1080  {
1081  cnt = n * nFieldDefPerDomain + i;
1082 
1083  FieldDef[cnt]->m_fields.push_back(variables[j]);
1084 
1085  m_vessels[n * m_nVariables]->AppendFieldData(
1086  FieldDef[cnt], FieldData[cnt],
1087  m_vessels[n * m_nVariables + j]->UpdateCoeffs());
1088  }
1089  }
1090 
1091  // Outputs pressure
1092  Array<OneD, NekDouble> PFwd(m_vessels[n * m_nVariables]->GetNcoeffs());
1093 
1094  m_vessels[n * m_nVariables]->FwdTrans_IterPerExp(m_pressure[n], PFwd);
1095 
1096  for (int i = 0; i < nFieldDefPerDomain; ++i)
1097  {
1098  cnt = n * nFieldDefPerDomain + i;
1099 
1100  FieldDef[cnt]->m_fields.push_back("P");
1101 
1102  m_vessels[n * m_nVariables]->AppendFieldData(FieldDef[cnt],
1103  FieldData[cnt], PFwd);
1104  }
1105 
1106  if (extraFields)
1107  {
1108  Array<OneD, NekDouble> PWVFwd(m_vessels[n *
1109  m_nVariables]->GetNcoeffs());
1110  Array<OneD, NekDouble> W1Fwd(m_vessels[n *
1111  m_nVariables]->GetNcoeffs());
1112  Array<OneD, NekDouble> W2Fwd(m_vessels[n *
1113  m_nVariables]->GetNcoeffs());
1114 
1115  m_vessels[n * m_nVariables]->FwdTrans_IterPerExp(m_PWV[n], PWVFwd);
1116  m_vessels[n * m_nVariables]->FwdTrans_IterPerExp(m_W1[n], W1Fwd);
1117  m_vessels[n * m_nVariables]->FwdTrans_IterPerExp(m_W2[n], W2Fwd);
1118 
1119  for (int i = 0; i < nFieldDefPerDomain; ++i)
1120  {
1121  cnt = n * nFieldDefPerDomain + i;
1122 
1123  FieldDef[cnt]->m_fields.push_back("c");
1124  FieldDef[cnt]->m_fields.push_back("W1");
1125  FieldDef[cnt]->m_fields.push_back("W2");
1126 
1127  m_vessels[n * m_nVariables]->AppendFieldData(FieldDef[cnt],
1128  FieldData[cnt], PWVFwd);
1129  m_vessels[n * m_nVariables]->AppendFieldData(FieldDef[cnt],
1130  FieldData[cnt], W1Fwd);
1131  m_vessels[n * m_nVariables]->AppendFieldData(FieldDef[cnt],
1132  FieldData[cnt], W2Fwd);
1133  }
1134  }
1135  }
1136 
1137  // Update time in field info if required
1138  if (m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
1139  {
1140  m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
1141  }
1142 
1143  LibUtilities::Write(outname, FieldDef, FieldData, m_fieldMetaDataMap);
1144 }
SOLVER_UTILS_EXPORT int GetNcoeffs()
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
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

References extraFields, Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::m_fieldMetaDataMap, m_nDomains, m_nVariables, m_pressure, m_PWV, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_time, m_vessels, m_W1, m_W2, and Nektar::LibUtilities::Write().

Referenced by CheckPoint_Output(), and v_Output().

Member Data Documentation

◆ extraFields

bool Nektar::PulseWaveSystem::extraFields = false
protected

Definition at line 131 of file PulseWaveSystem.h.

Referenced by Nektar::PulseWavePropagation::GetFluxVector(), and WriteVessels().

◆ m_A_0

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

◆ m_A_0_trace

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

◆ m_alpha

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

◆ m_alpha_trace

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

◆ m_beta

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

◆ m_beta_trace

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

◆ m_bifurcations

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

Definition at line 137 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

◆ 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_gamma

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

Definition at line 122 of file PulseWaveSystem.h.

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

◆ m_mergingJcts

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

Definition at line 138 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_nDomains

int Nektar::PulseWaveSystem::m_nDomains
protected

◆ m_nVariables

int Nektar::PulseWaveSystem::m_nVariables
protected

◆ m_pext

NekDouble Nektar::PulseWaveSystem::m_pext
protected

Definition at line 112 of file PulseWaveSystem.h.

Referenced by v_InitObject().

◆ m_pout

NekDouble Nektar::PulseWaveSystem::m_pout
protected

Definition at line 116 of file PulseWaveSystem.h.

◆ m_pressure

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

◆ m_pressureArea

PulseWavePressureAreaSharedPtr Nektar::PulseWaveSystem::m_pressureArea
protected

◆ m_PWV

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

◆ m_rho

NekDouble Nektar::PulseWaveSystem::m_rho
protected

◆ 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 125 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 136 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_vessels

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

◆ m_W1

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

◆ m_W2

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