Nektar++
Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | 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

int GetNdomains ()
 
Array< OneD, MultiRegions::ExpListSharedPtrUpdateVessels (void)
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT ~UnsteadySystem () override=default
 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 NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeSharedPtrGetTimeIntegrationScheme ()
 Returns the time integration scheme. More...
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeOperatorsGetTimeIntegrationSchemeOperators ()
 Returns the time integration scheme operators. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void InitObject (bool DeclareField=true)
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise (bool dumpInitialConditions=true)
 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 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 NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT NekDouble H1Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the H1 error between fields and a given exact solution. 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 GetTime ()
 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 void SetSteps (const int steps)
 
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, const Array< OneD, const NekDouble > &input)
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > & UpdatePhysField (const int i)
 
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 int GetInfoSteps ()
 
SOLVER_UTILS_EXPORT void SetInfoSteps (int num)
 
SOLVER_UTILS_EXPORT void SetIterationNumberPIT (int num)
 
SOLVER_UTILS_EXPORT void SetWindowNumberPIT (int num)
 
SOLVER_UTILS_EXPORT Array< OneD, const Array< OneD, NekDouble > > GetTraceNormals ()
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
SOLVER_UTILS_EXPORT bool NegatedOp ()
 Identify if operator is negated in DoSolve. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::ALEHelper
virtual ~ALEHelper ()=default
 
virtual SOLVER_UTILS_EXPORT void v_ALEInitObject (int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 
SOLVER_UTILS_EXPORT void InitObject (int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 
virtual SOLVER_UTILS_EXPORT void v_UpdateGridVelocity (const NekDouble &time)
 
virtual SOLVER_UTILS_EXPORT void v_ALEPreMultiplyMass (Array< OneD, Array< OneD, NekDouble > > &fields)
 
SOLVER_UTILS_EXPORT void ALEDoElmtInvMass (Array< OneD, Array< OneD, NekDouble > > &traceNormals, Array< OneD, Array< OneD, NekDouble > > &fields, NekDouble time)
 Update m_fields with u^n by multiplying by inverse mass matrix. That's then used in e.g. checkpoint output and L^2 error calculation. More...
 
SOLVER_UTILS_EXPORT void ALEDoElmtInvMassBwdTrans (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
SOLVER_UTILS_EXPORT void MoveMesh (const NekDouble &time, Array< OneD, Array< OneD, NekDouble > > &traceNormals)
 
const Array< OneD, const Array< OneD, NekDouble > > & GetGridVelocity ()
 
SOLVER_UTILS_EXPORT const Array< OneD, const Array< OneD, NekDouble > > & GetGridVelocityTrace ()
 
SOLVER_UTILS_EXPORT void ExtraFldOutputGridVelocity (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 

Protected Member Functions

 PulseWaveSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises PulseWaveSystem class members. More...
 
 ~PulseWaveSystem () override=default
 
void v_InitObject (bool DeclareField=false) override
 
void v_DoInitialise (bool dumpInitialConditions=false) override
 Sets up initial conditions. More...
 
void v_DoSolve () override
 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 InterfaceRiemann (Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0, Array< OneD, NekDouble > &alpha)
 Riemann Problem for Interface/Junction. More...
 
void v_Output (void) override
 
void CheckPoint_Output (const int n)
 
NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false) override
 Compute the L2 error between fields and a given exact solution. More...
 
NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray) override
 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 void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. More...
 
SOLVER_UTILS_EXPORT void v_DoSolve () override
 Solves an unsteady problem. More...
 
virtual SOLVER_UTILS_EXPORT void v_PrintStatusInformation (const int step, const NekDouble cpuTime)
 Print Status Information. More...
 
virtual SOLVER_UTILS_EXPORT void v_PrintSummaryStatistics (const NekDouble intTime)
 Print Summary Statistics. More...
 
SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true) override
 Sets up initial conditions. More...
 
SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &s) override
 Print a summary of time stepping parameters. 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)
 
virtual SOLVER_UTILS_EXPORT bool v_UpdateTimeStepCheck ()
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. More...
 
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...
 
SOLVER_UTILS_EXPORT void DoDummyProjection (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
 Perform dummy projection. 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_InitObject (bool DeclareFeld=true)
 Initialisation object for EquationSystem. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true)
 Virtual function for initialisation implementation. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoSolve ()
 Virtual function for solve implementation. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Virtual function for the L_inf error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the L_2 error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_H1Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the H_1 error computation between fields and a given exact solution. 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_GenerateSummary (SummaryList &l)
 Virtual function for generating summary information. 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 void v_Output (void)
 
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure (void)
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp (void)
 Virtual function to identify if operator is negated in DoSolve. More...
 
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
 
size_t m_nDomains
 
size_t m_currentDomain
 
size_t 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
 
std::map< int, SpatialDomains::CompositeMapm_domain
 
std::vector< int > m_domOrder
 
std::map< int, std::vector< int > > m_domainToFilterIDs
 
std::map< int, int > m_filterToVesselID
 
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_vesselIntfcs
 
std::vector< std::vector< InterfacePointShPtr > > m_bifurcations
 
std::vector< std::vector< InterfacePointShPtr > > m_mergingJcts
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check. More...
 
std::vector< int > m_intVariables
 
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 
NekDouble m_CFLGrowth
 CFL growth rate. More...
 
NekDouble m_CFLEnd
 Maximun cfl in cfl growth. More...
 
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
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...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
std::ofstream m_errFile
 
NekDouble m_epsilon
 Diffusion coefficient. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
bool m_verbose
 
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_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_infosteps
 Number of time steps between outputting status information. More...
 
int m_iterPIT = 0
 Number of parallel-in-time time iteration. More...
 
int m_windowPIT = 0
 Index of windows for parallel-in-time time iteration. 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...
 
Array< OneD, NekDoublem_movingFrameData
 Moving reference frame status in the inertial frame X, Y, Z, Theta_x, Theta_y, Theta_z, U, V, W, Omega_x, Omega_y, Omega_z, A_x, A_y, A_z, DOmega_x, DOmega_y, DOmega_z, pivot_x, pivot_y, pivot_z. More...
 
std::vector< std::string > m_strFrameData
 variable name in m_movingFrameData 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...
 
- Protected Attributes inherited from Nektar::SolverUtils::ALEHelper
Array< OneD, MultiRegions::ExpListSharedPtrm_fieldsALE
 
Array< OneD, Array< OneD, NekDouble > > m_gridVelocity
 
Array< OneD, Array< OneD, NekDouble > > m_gridVelocityTrace
 
std::vector< ALEBaseShPtrm_ALEs
 
bool m_ALESolver = false
 
bool m_ImplicitALESolver = false
 
NekDouble m_prevStageTime = 0.0
 
int m_spaceDim
 

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)
 
void GetCommArray (std::map< int, LibUtilities::CommSharedPtr > &retval)
 Set and retrn a series of communicators for each partition. More...
 
void SetUpDomainInterfaceBCs (SpatialDomains::BoundaryConditions &AllBcs, std::map< int, LibUtilities::CommSharedPtr > &domComm)
 

Private Attributes

Gs::gs_datam_intComm
 Communicator for interfaces. More...
 

Additional Inherited Members

- Static Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
static std::string cmdSetStartTime
 
static std::string cmdSetStartChkNum
 
- 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 []
 
static std::string projectionTypeLookupIds []
 

Detailed Description

Base class for unsteady solvers.

Initialises the arterial subdomains in m_vessels and sets up all domain-linking conditions (bifurcations, intefaces/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 77 of file PulseWaveSystem.h.

Constructor & Destructor Documentation

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

◆ ~PulseWaveSystem()

Nektar::PulseWaveSystem::~PulseWaveSystem ( )
overrideprotecteddefault

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

1594{
1595
1596 NekDouble rho = m_rho;
1597 Array<OneD, NekDouble> W(3);
1598 Array<OneD, NekDouble> P_Au(3);
1599 Array<OneD, NekDouble> W_Au(3);
1600 NekMatrix<NekDouble> invJ(6, 6);
1601 NekVector<NekDouble> f(6);
1602 NekVector<NekDouble> dx(6);
1603
1604 int proceed = 1;
1605 int iter = 0;
1606 int MAX_ITER = 100;
1607
1608 // Forward characteristic
1609 m_pressureArea->GetW1(W[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1610
1611 // Backward characteristics
1612 for (size_t i = 1; i < 3; ++i)
1613 {
1614 m_pressureArea->GetW2(W[i], uu[i], beta[i], Au[i], A_0[i], alpha[i]);
1615 }
1616
1617 // Tolerances for the algorithm
1618 NekDouble Tol = 1.0E-10;
1619
1620 // Newton Iteration
1621 while ((proceed) && (iter < MAX_ITER))
1622 {
1623 LibUtilities::Timer time_BifurcationRiemann;
1624 time_BifurcationRiemann.Start();
1625 iter += 1;
1626
1627 /*
1628 * We solve the six constraint equations via a multivariate Newton
1629 * iteration. Equations are:
1630 * 1. Forward characteristic: W1(A_L, U_L) = W1(Au_L,
1631 * Uu_L)
1632 * 2. Backward characteristic 1: W2(A_R1, U_R1) = W2(Au_R1,
1633 * Uu_R1)
1634 * 3. Backward characteristic 2: W2(A_R2, U_R2) = W2(Au_R2,
1635 * Uu_R2)
1636 * 4. Conservation of mass: Au_L * Uu_L = Au_R1 *
1637 * Uu_R1 + Au_R2 * Uu_R2
1638 * 5. Continuity of total pressure 1: rho * Uu_L * Uu_L / 2 +
1639 * p(Au_L) = rho * Uu_R1 * Uu_R1 / 2 + p(Au_R1)
1640 * 6. Continuity of total pressure 2: rho * Uu_L * Uu_L / 2 +
1641 * p(Au_L) = rho * Uu_R2 * Uu_R2 / 2 + p(Au_R2)
1642 */
1643 for (size_t i = 0; i < 3; ++i)
1644 {
1645 m_pressureArea->GetPressure(P_Au[i], beta[i], Au[i], A_0[i], 0, 0,
1646 alpha[i]);
1647 }
1648
1649 m_pressureArea->GetW1(W_Au[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1650 m_pressureArea->GetW2(W_Au[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
1651 m_pressureArea->GetW2(W_Au[2], uu[2], beta[2], Au[2], A_0[2], alpha[2]);
1652
1653 // Constraint equations set to zero
1654 f[0] = W_Au[0] - W[0];
1655 f[1] = W_Au[1] - W[1];
1656 f[2] = W_Au[2] - W[2];
1657 f[3] = Au[0] * uu[0] - Au[1] * uu[1] - Au[2] * uu[2];
1658 f[4] = uu[0] * uu[0] + 2 * P_Au[0] / rho - uu[1] * uu[1] -
1659 2 * P_Au[1] / rho;
1660 f[5] = uu[0] * uu[0] + 2 * P_Au[0] / rho - uu[2] * uu[2] -
1661 2 * P_Au[2] / rho;
1662
1663 // Inverse Jacobian for x + dx = x - J^(-1) * f(x)
1664 m_pressureArea->GetJacobianInverse(invJ, Au, uu, beta, A_0, alpha,
1665 "Bifurcation");
1666
1667 Multiply(dx, invJ, f);
1668
1669 // Update the solution: x_new = x_old - dx
1670 for (int i = 0; i < 3; ++i)
1671 {
1672 uu[i] -= dx[i];
1673 Au[i] -= dx[i + 3];
1674 }
1675
1676 // Check if the error of the solution is smaller than Tol
1677 if (Dot(dx, dx) < Tol)
1678 {
1679 proceed = 0;
1680 }
1681
1682 // Check if solver converges
1683 if (iter >= MAX_ITER)
1684 {
1685 ASSERTL0(false, "Riemann solver for Bifurcation did not converge.");
1686 }
1687 time_BifurcationRiemann.Stop();
1688 time_BifurcationRiemann.AccumulateRegion(
1689 "PulseWaveSystem::Bifurcation Riemann", 2);
1690 }
1691}
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
PulseWavePressureAreaSharedPtr m_pressureArea
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
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:1193

References Nektar::LibUtilities::Timer::AccumulateRegion(), ASSERTL0, Nektar::LibUtilities::beta, Nektar::Dot(), m_pressureArea, m_rho, Nektar::Multiply(), Nektar::LibUtilities::Timer::Start(), and Nektar::LibUtilities::Timer::Stop().

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

1920{
1921 std::stringstream outname;
1922 outname << m_sessionName << "_" << n << ".chk";
1923
1924 WriteVessels(outname.str());
1925}
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 1445 of file PulseWaveSystem.cpp.

1447{
1448 LibUtilities::Timer time_EnforceInterfaceConditions;
1449 time_EnforceInterfaceConditions.Start();
1450 int dom, bcpos;
1451
1452 int totif =
1453 m_bifurcations.size() + m_mergingJcts.size() + m_vesselIntfcs.size();
1454
1455 Array<OneD, NekDouble> Aut, Au(3 * totif, 0.0);
1456 Array<OneD, NekDouble> uut, uu(3 * totif, 0.0);
1457 Array<OneD, NekDouble> betat, beta(3 * totif, 0.0);
1458 Array<OneD, NekDouble> A_0t, A_0(3 * totif, 0.0);
1459 Array<OneD, NekDouble> alphat, alpha(3 * totif, 0.0);
1460
1461 // Bifurcations Data:
1462 size_t cnt = 0;
1463 for (size_t n = 0; n < m_bifurcations.size(); ++n, ++cnt)
1464 {
1465 for (size_t i = 0; i < m_bifurcations[n].size(); ++i)
1466 {
1467 size_t l = m_bifurcations[n][i]->m_riemannOrd;
1469 m_bifurcations[n][i], fields, Au[l + 3 * cnt], uu[l + 3 * cnt],
1470 beta[l + 3 * cnt], A_0[l + 3 * cnt], alpha[l + 3 * cnt]);
1471 }
1472 }
1473
1474 // Enforce Merging vessles Data:
1475 for (size_t n = 0; n < m_mergingJcts.size(); ++n, ++cnt)
1476 {
1477 // Merged vessel
1478 for (size_t i = 0; i < m_mergingJcts.size(); ++i)
1479 {
1480 int l = m_mergingJcts[n][i]->m_riemannOrd;
1482 m_mergingJcts[n][i], fields, Au[l + 3 * cnt], uu[l + 3 * cnt],
1483 beta[l + 3 * cnt], A_0[l + 3 * cnt], alpha[l + 3 * cnt]);
1484 }
1485 }
1486
1487 // Enforce interface junction between two vessesls
1488 for (size_t n = 0; n < m_vesselIntfcs.size(); ++n, ++cnt)
1489 {
1490 for (size_t i = 0; i < m_vesselIntfcs[n].size(); ++i)
1491 {
1492 int l = m_vesselIntfcs[n][i]->m_riemannOrd;
1494 m_vesselIntfcs[n][i], fields, Au[l + 3 * cnt], uu[l + 3 * cnt],
1495 beta[l + 3 * cnt], A_0[l + 3 * cnt], alpha[l + 3 * cnt]);
1496 }
1497 }
1498
1499 // Gather data if running in parallel
1505
1506 // Enforce Bifurcations:
1507 cnt = 0;
1508 for (size_t n = 0; n < m_bifurcations.size(); ++n, ++cnt)
1509 {
1510 // Solve the Riemann problem for a bifurcation
1511 BifurcationRiemann(Aut = Au + 3 * cnt, uut = uu + 3 * cnt,
1512 betat = beta + 3 * cnt, A_0t = A_0 + 3 * cnt,
1513 alphat = alpha + 3 * cnt);
1514
1515 // Store the values into the right positions:
1516 for (size_t i = 0; i < m_bifurcations[n].size(); ++i)
1517 {
1518 dom = m_bifurcations[n][i]->m_domain;
1519 bcpos = m_bifurcations[n][i]->m_bcPosition;
1520 int l = m_bifurcations[n][i]->m_riemannOrd;
1521 m_vessels[dom * m_nVariables]
1522 ->UpdateBndCondExpansion(bcpos)
1523 ->UpdatePhys()[0] = Au[l + 3 * cnt];
1524 m_vessels[dom * m_nVariables + 1]
1525 ->UpdateBndCondExpansion(bcpos)
1526 ->UpdatePhys()[0] = uu[l + 3 * cnt];
1527 }
1528 }
1529
1530 // Enforce Merging vessles;
1531 for (size_t n = 0; n < m_mergingJcts.size(); ++n, ++cnt)
1532 {
1533 // Solve the Riemann problem for a merging vessel
1534 MergingRiemann(Aut = Au + 3 * cnt, uut = uu + 3 * cnt,
1535 betat = beta + 3 * cnt, A_0t = A_0 + 3 * cnt,
1536 alphat = alpha + 3 * cnt);
1537
1538 // Store the values into the right positions:
1539 for (size_t i = 0; i < m_mergingJcts[n].size(); ++i)
1540 {
1541 int dom = m_mergingJcts[n][i]->m_domain;
1542 int bcpos = m_mergingJcts[n][i]->m_bcPosition;
1543 int l = m_mergingJcts[n][i]->m_riemannOrd;
1544 m_vessels[dom * m_nVariables]
1545 ->UpdateBndCondExpansion(bcpos)
1546 ->UpdatePhys()[0] = Au[l + 3 * cnt];
1547 m_vessels[dom * m_nVariables + 1]
1548 ->UpdateBndCondExpansion(bcpos)
1549 ->UpdatePhys()[0] = uu[l + 3 * cnt];
1550 }
1551 }
1552
1553 // Enforce interface junction between two vessesls
1554 for (size_t n = 0; n < m_vesselIntfcs.size(); ++n, ++cnt)
1555 {
1556 InterfaceRiemann(Aut = Au + 3 * cnt, uut = uu + 3 * cnt,
1557 betat = beta + 3 * cnt, A_0t = A_0 + 3 * cnt,
1558 alphat = alpha + 3 * cnt);
1559
1560 // Store the values into the right positions:
1561 for (size_t i = 0; i < m_vesselIntfcs[n].size(); ++i)
1562 {
1563 int dom = m_vesselIntfcs[n][i]->m_domain;
1564 int bcpos = m_vesselIntfcs[n][i]->m_bcPosition;
1565 int l = m_vesselIntfcs[n][i]->m_riemannOrd;
1566 m_vessels[dom * m_nVariables]
1567 ->UpdateBndCondExpansion(bcpos)
1568 ->UpdatePhys()[0] = Au[l + 3 * cnt];
1569 m_vessels[dom * m_nVariables + 1]
1570 ->UpdateBndCondExpansion(bcpos)
1571 ->UpdatePhys()[0] = uu[l + 3 * cnt];
1572 }
1573 }
1574 time_EnforceInterfaceConditions.Stop();
1575 time_EnforceInterfaceConditions.AccumulateRegion(
1576 "PulseWaveSystem::EnforceInterfaceConditions", 1);
1577}
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.
void InterfaceRiemann(Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0, Array< OneD, NekDouble > &alpha)
Riemann Problem for Interface/Junction.
Gs::gs_data * m_intComm
Communicator for interfaces.
std::vector< std::vector< InterfacePointShPtr > > m_vesselIntfcs
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_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.
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
Definition: GsLib.hpp:278
@ gs_add
Definition: GsLib.hpp:60

References Nektar::LibUtilities::Timer::AccumulateRegion(), Nektar::LibUtilities::beta, BifurcationRiemann(), FillDataFromInterfacePoint(), Gs::Gather(), Gs::gs_add, InterfaceRiemann(), m_bifurcations, m_intComm, m_mergingJcts, m_nVariables, m_vesselIntfcs, m_vessels, MergingRiemann(), Nektar::LibUtilities::Timer::Start(), and Nektar::LibUtilities::Timer::Stop().

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

1417{
1418 LibUtilities::Timer time_FillDataFromInterfacePoint;
1419 time_FillDataFromInterfacePoint.Start();
1420
1421 int omega = I->m_domain;
1422 int traceId = I->m_traceId;
1423 int eid = I->m_elmt;
1424 int vert = I->m_elmtVert;
1425 int vesselID = omega * m_nVariables;
1426 int phys_offset = m_vessels[vesselID]->GetPhys_Offset(eid);
1428 Array<OneD, NekDouble> Tmp(1);
1429
1430 m_vessels[vesselID]->GetExp(eid)->GetTracePhysVals(
1431 vert, dumExp, fields[0] + m_fieldPhysOffset[omega] + phys_offset, Tmp);
1432 A = Tmp[0];
1433 m_vessels[vesselID]->GetExp(eid)->GetTracePhysVals(
1434 vert, dumExp, fields[1] + m_fieldPhysOffset[omega] + phys_offset, Tmp);
1435 u = Tmp[0];
1436
1437 beta = m_beta_trace[omega][traceId];
1438 A_0 = m_A_0_trace[omega][traceId];
1439 alpha = m_alpha_trace[omega][traceId];
1440 time_FillDataFromInterfacePoint.Stop();
1441 time_FillDataFromInterfacePoint.AccumulateRegion(
1442 "PulseWaveSystem::FillDataFromInterfacePoint", 2);
1443}
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:66

References Nektar::LibUtilities::Timer::AccumulateRegion(), Nektar::LibUtilities::beta, m_A_0_trace, m_alpha_trace, m_beta_trace, m_fieldPhysOffset, m_nVariables, m_vessels, Nektar::LibUtilities::Timer::Start(), and Nektar::LibUtilities::Timer::Stop().

Referenced by EnforceInterfaceConditions().

◆ GetCommArray()

void Nektar::PulseWaveSystem::GetCommArray ( std::map< int, LibUtilities::CommSharedPtr > &  retval)
private

Set and retrn a series of communicators for each partition.

Definition at line 347 of file PulseWaveSystem.cpp.

349{
350
351 // set up defualt serial communicator
352 std::string def("default");
353 char *argv = new char[def.length() + 1];
354 std::strcpy(argv, def.c_str());
355 LibUtilities::CommSharedPtr serialComm =
357
358 size_t nprocs = m_comm->GetSize();
359
360 if (nprocs == 1) // serial case
361 {
362 for (auto &d : m_domain)
363 {
364 // serial communicator
365 retval[d.first] = serialComm;
366 // fill out m_domOrder set for serial case
367 m_domOrder.push_back(d.first);
368 }
369 }
370 else // parallel case
371 {
372 size_t rank = m_comm->GetRank();
373
374 // Fill array with domain details
375 size_t dmax = 0;
376 for (auto &d : m_domain)
377 {
378 dmax = ((size_t)d.first > dmax) ? d.first : dmax;
379 }
380 dmax += 1;
381
382 // get domain maximum id + 1;
383 m_comm->AllReduce(dmax, LibUtilities::ReduceMax);
384
385 Array<OneD, int> commtmp(dmax, 0);
386 for (auto &d : m_domain)
387 {
388 commtmp[d.first] = 1;
389 }
390
391 // identfy which domains are split between two procs
392 m_comm->AllReduce(commtmp, LibUtilities::ReduceSum);
393
394 Array<OneD, bool> DoneCreate(nprocs, false);
395
396 // setup communicators for domain partitions greater than 2
397 for (size_t d = 0; d < dmax; ++d)
398 {
399 // set up communicator for cases where more than 2 procs
400 // share domain
401 if (commtmp[d] > 2)
402 {
404 int flag = 0;
405 // set flag to 10 if d is on this processors domain
406 for (auto &d1 : m_domain)
407 {
408 if ((size_t)d1.first == d)
409 {
410 flag = 10;
411 }
412 }
413
414 newcomm = m_comm->CommCreateIf(flag);
415
416 // set up communicator if domain on this processor
417 for (auto &d1 : m_domain)
418 {
419 if ((size_t)d1.first == d)
420 {
421 retval[d] = newcomm;
422 }
423 }
424
425 // turn off commtmp flag since should now be setup
426 commtmp[d] = -1;
427 }
428 }
429
430 // we now should have a list of the number of each
431 // domains/vessel that is split over two processors
432
433 // set up serial communicators for domains on only one
434 // processor and set up map of domains that require parallel
435 // communicator and number of processors in parallel
436
437 // communicator
438 std::map<int, int> SharedProc;
439 Array<OneD, int> commtmp1(dmax, 0);
440
441 for (auto &d : m_domain)
442 {
443 if (commtmp[d.first] == 1) // set to serial communicator
444 {
445 retval[d.first] = serialComm;
446 }
447 else if (commtmp[d.first] == 2) // shared by 2 procs
448 {
449 SharedProc[d.first] = 1;
450 commtmp1[d.first] = rank;
451 }
452 else // do nothing
453 {
454 }
455 }
456
457 // find the maximum number of communicator on all processors
458 size_t nShrProc = SharedProc.size();
459 int commMax = nShrProc;
460 m_comm->AllReduce(commMax, LibUtilities::ReduceMax);
461
462 // find out processor id of each communicator by adding rank
463 // and then subtracting local rank
464
465 m_comm->AllReduce(commtmp1, LibUtilities::ReduceSum);
466
467 // reset SharedProc to hold shared proc id.
468 for (auto &d : SharedProc)
469 {
470 SharedProc[d.first] = commtmp1[d.first] - rank;
471 }
472
473 // now work out a comm ordering;
474 Array<OneD, int> procs(2);
475 int commcall = 0;
476 std::map<int, std::pair<int, int>> CreateComm;
477
478 bool search = true;
479 Array<OneD, int> setdone(3);
480
481 while (search)
482 {
483 size_t dorank = 0; // search index over processors
484 std::set<int> proclist; // has proc been identified for a comm at
485 // this level
486 int flag = 100; // colour for communicators in this search level
487
488 while (dorank < nprocs)
489 {
490 int sharedproc = -1;
491 int ncomms = 0;
492
493 if (DoneCreate[dorank] == false)
494 {
495 // check to see if proc even has a shared domain
496 ncomms = (rank == dorank) ? nShrProc : 0;
497 m_comm->AllReduce(ncomms, LibUtilities::ReduceMax);
498
499 // set DoneCreate to true if n comms required on
500 // proc
501 DoneCreate[dorank] = (ncomms == 0) ? true : false;
502 }
503
504 // requires communicator and needs to create one.
505 if (ncomms)
506 {
507 bool doneCreateOnDoRank = false;
508 bool doneCreateOnSharedProc = false;
509 bool createdComm = false;
510
511 // check not already callng communicator on
512 // this processor for this search sweep
513 if ((proclist.count(dorank) == 0))
514 {
515
516 // give all processors rank id of shared proc
517 if (rank == dorank)
518 {
519 for (auto &d : SharedProc)
520 {
521 if (CreateComm.count(d.first) == 0)
522 {
523 sharedproc = SharedProc[d.first];
524 break;
525 }
526 }
527
528 ASSERTL1(sharedproc != -1,
529 "Failed to fine a new "
530 "processor to set up another "
531 "communicator");
532
533 if (proclist.count(sharedproc) == 0)
534 {
535 // save all communicators on this for
536 // split comm setup
537 for (auto &d : SharedProc)
538 {
539 if (d.second == sharedproc)
540 {
541 CreateComm[d.first] =
542 std::pair<int, int>(commcall,
543 flag + d.first);
544 createdComm = true;
545 }
546 }
547 }
548
549 // set bool to mark as done on dorank
550 if (CreateComm.size() == nShrProc)
551 {
552 doneCreateOnDoRank = true;
553 }
554 }
555 else
556 {
557 sharedproc = -1;
558 }
559
560 m_comm->AllReduce(sharedproc, LibUtilities::ReduceMax);
561 ASSERTL1(
562 sharedproc < nprocs,
563 "Error in obtaining rank "
564 "shared by domain. Sharedproc value is greater "
565 "than number of procs");
566
567 if (proclist.count(sharedproc) == 0)
568 {
569 if ((int)rank == sharedproc)
570 {
571 // save all communicators on this for
572 // split comm setup
573 for (auto &d : SharedProc)
574 {
575 if (d.second == (int)dorank)
576 {
577 CreateComm[d.first] =
578 std::pair<int, int>(commcall,
579 flag + d.first);
580 createdComm = true;
581 }
582 }
583
584 if (CreateComm.size() == nShrProc)
585 {
586 doneCreateOnSharedProc = true;
587 }
588 }
589 }
590 }
591
592 // mark DoneCreate as if complete
593 setdone[0] = (doneCreateOnDoRank) ? 1 : 0;
594 setdone[1] = (doneCreateOnSharedProc) ? 1 : 0;
595 setdone[2] = (createdComm) ? 1 : 0;
596 m_comm->AllReduce(setdone, LibUtilities::ReduceMax);
597 DoneCreate[dorank] = (bool)setdone[0];
598 if (sharedproc != -1)
599 {
600 DoneCreate[sharedproc] = (bool)setdone[1];
601 }
602 if (setdone[2]) // created new communicator
603 {
604 proclist.insert(dorank);
605 proclist.insert(sharedproc);
606 }
607 }
608
609 dorank++;
610 }
611
612 // have we found all comms on all processors.
613 size_t i;
614 for (i = 0; i < nprocs; ++i)
615 {
616 if (DoneCreate[i] == false)
617 {
618 break;
619 }
620 }
621
622 if (i == nprocs)
623 {
624 search = false;
625 }
626 else
627 {
628 commcall++;
629 ASSERTL0(commcall < 4 * commMax,
630 "Failed to find sub communicators "
631 "pattern in 4*commMax searches");
632 }
633 }
634
635 ASSERTL1(CreateComm.size() == nShrProc,
636 "Have not created communicators for all shared procs");
637
638 // determine maxmimum number of CreateComm size
639 size_t maxCreateComm = CreateComm.size();
640 m_comm->AllReduce(maxCreateComm, LibUtilities::ReduceMax);
641
642 // loop over CreateComm list
643 for (size_t i = 0; i < maxCreateComm; ++i)
644 {
646 int flag = 0;
647
648 for (auto &d : CreateComm)
649 {
650 if (d.second.first == (int)i)
651 {
652 flag = d.second.second;
653 }
654 }
655
656 newcomm = m_comm->CommCreateIf(flag);
657
658 for (auto &d : CreateComm)
659 {
660 if (d.second.first == (int)i)
661 {
662 retval[d.first] = newcomm;
663 }
664 }
665 }
666
667 // Finally need to re-order domain list so that we call the
668 // communicators in a non-locking manner. This is done by
669 // uniquely numbering each shared proc/comm interface and
670 // ensuring the domains are ordered so that we call this
671 // numbering in an increasing order
672
673 Array<OneD, NekDouble> numShared(nprocs, 0.0);
674 // determine the number of communicators on this rank where
675 // the share proc has a higher rank id
676 numShared[rank] = nShrProc;
677 m_comm->AllReduce(numShared, LibUtilities::ReduceSum);
678 int totShared = (int)Vmath::Vsum(nprocs, numShared, 1);
679
680 if (totShared)
681 {
682 size_t cnt = 0;
683 Array<OneD, NekDouble> numOffset(nprocs, 0.0);
684 for (auto &s : SharedProc)
685 {
686 if (s.second > (int)rank)
687 {
688 cnt++;
689 }
690 }
691
692 numOffset[rank] = cnt;
693 m_comm->AllReduce(numOffset, LibUtilities::ReduceMax);
694
695 // make numShared into a cumulative list
696 for (size_t i = 1; i < nprocs; ++i)
697 {
698 numShared[i] += numShared[i - 1];
699 numOffset[i] += numOffset[i - 1];
700 }
701 for (size_t i = nprocs - 1; i > 0; --i)
702 {
703 numShared[i] = numShared[i - 1];
704 numOffset[i] = numOffset[i - 1];
705 }
706 numShared[0] = 0;
707 numOffset[0] = 0;
708
709 Array<OneD, NekDouble> shareddom(totShared, -1.0);
710 cnt = 0; // collect a list of domains that each shared commm is
711 // attached to
712 for (auto &s : SharedProc)
713 {
714 shareddom[numShared[rank] + cnt] = s.first;
715 ++cnt;
716 }
717 m_comm->AllReduce(shareddom, LibUtilities::ReduceMax);
718
719 // define a numbering scheme
720 Array<OneD, NekDouble> sharedid(totShared, -1.0);
721 cnt = 0;
722 size_t cnt1 = 0;
723 for (auto &s : SharedProc)
724 {
725 if (s.second > (int)rank)
726 {
727 sharedid[numShared[rank] + cnt] = numOffset[rank] + cnt1;
728
729 // find shared proc offset by matching the domain ids.
730 size_t j;
731 for (j = 0; j < maxCreateComm; ++j)
732 {
733 if ((numShared[s.second] + j < totShared) &&
734 (shareddom[numShared[s.second] + j] == s.first))
735 {
736 break;
737 }
738 }
739
740 sharedid[numShared[s.second] + j] = numOffset[rank] + cnt1;
741 cnt1++;
742 }
743 cnt++;
744 }
745
746 // now communicate ids to all procesors.
747 m_comm->AllReduce(sharedid, LibUtilities::ReduceMax);
748
749 if (rank == 0)
750 {
751 for (size_t i = 0; i < (size_t)totShared; ++i)
752 {
753 ASSERTL1(sharedid[i] != -1.0,
754 "Failed to number shared proc uniquely");
755 }
756 }
757
758 cnt = 0;
759 int maxoffset =
760 (int)Vmath::Vmax(nShrProc, &sharedid[numShared[rank]], 1);
761 m_comm->AllReduce(maxoffset, LibUtilities::ReduceMax);
762 maxoffset++;
763
764 // make a map relating the order of the SharedProc to the domain
765 // id
766 std::map<int, int> ShrToDom;
767 cnt = 0;
768 for (auto &s : SharedProc)
769 {
770 ShrToDom[cnt] = s.first;
771 ++cnt;
772 }
773
774 // Set up a set the domain ids listed in the ordering of
775 // the SharedProc unique numbering
776 std::set<int> doneDom;
777 for (size_t i = 0; i < nShrProc; ++i)
778 {
779 int minId =
780 Vmath::Imin(nShrProc, &sharedid[numShared[rank]], 1);
781 sharedid[numShared[rank] + minId] += maxoffset;
782
783 m_domOrder.push_back(ShrToDom[minId]);
784 doneDom.insert(ShrToDom[minId]);
785 }
786
787 // add all the other
788 for (auto &d : m_domain)
789 {
790 if (doneDom.count(d.first) == 0)
791 {
792 m_domOrder.push_back(d.first);
793 }
794 }
795 }
796 }
797}
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::vector< int > m_domOrder
std::map< int, SpatialDomains::CompositeMap > m_domain
LibUtilities::CommSharedPtr m_comm
Communicator.
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
std::vector< double > d(NPUPPER *NPUPPER)
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.hpp:704
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.hpp:644

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ASSERTL1, Nektar::UnitTests::d(), Vmath::Imin(), Nektar::SolverUtils::EquationSystem::m_comm, m_domain, m_domOrder, Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceSum, Vmath::Vmax(), and Vmath::Vsum().

Referenced by v_InitObject().

◆ GetNdomains()

int Nektar::PulseWaveSystem::GetNdomains ( )
inline

Definition at line 80 of file PulseWaveSystem.h.

81 {
82 return m_nDomains;
83 }

References m_nDomains.

◆ InterfaceRiemann()

void Nektar::PulseWaveSystem::InterfaceRiemann ( 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 Interface/Junction.

Solves the Riemann problem at an interdomain junction/Interface 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 1816 of file PulseWaveSystem.cpp.

1821{
1822
1823 NekDouble rho = m_rho;
1824 Array<OneD, NekDouble> W(2);
1825 Array<OneD, NekDouble> W_Au(2);
1826 Array<OneD, NekDouble> P_Au(2);
1827 NekMatrix<NekDouble> invJ(4, 4);
1828 NekVector<NekDouble> f(4);
1829 NekVector<NekDouble> dx(4);
1830
1831 int proceed = 1;
1832 int iter = 0;
1833 int MAX_ITER = 15;
1834 NekDouble Tol = 1.0E-10;
1835
1836 // Forward and backward characteristics
1837 m_pressureArea->GetW1(W[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1838 m_pressureArea->GetW2(W[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
1839
1840 while ((proceed) && (iter < MAX_ITER))
1841 {
1842 LibUtilities::Timer time_InterfaceRiemann;
1843 time_InterfaceRiemann.Start();
1844 iter += 1;
1845
1846 /*
1847 * We solve the four constraint equations via a multivariate Newton
1848 * iteration. Equations are:
1849 * 1. Forward characteristic: W1(A_L, U_L) = W1(Au_L, Uu_L)
1850 * 2. Backward characteristic: W2(A_R, U_R) = W2(Au_R, Uu_R)
1851 * 3. Conservation of mass: Au_L * Uu_L = Au_R * Uu_R
1852 * 4. Continuity of total pressure: rho * Uu_L * Uu_L / 2 + p(Au_L)
1853 * = rho * Uu_R * Uu_R / 2 + p(Au_R)
1854 */
1855 for (size_t i = 0; i < 2; ++i)
1856 {
1857 m_pressureArea->GetPressure(P_Au[i], beta[i], Au[i], A_0[i], 0, 0,
1858 alpha[i]);
1859 }
1860
1861 m_pressureArea->GetW1(W_Au[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1862 m_pressureArea->GetW2(W_Au[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
1863
1864 // Constraint equations set to zero
1865 f[0] = W_Au[0] - W[0];
1866 f[1] = W_Au[1] - W[1];
1867 f[2] = Au[0] * uu[0] - Au[1] * uu[1];
1868 f[3] = uu[0] * uu[0] + 2 * P_Au[0] / rho - uu[1] * uu[1] -
1869 2 * P_Au[1] / rho;
1870
1871 // Inverse Jacobian for x + dx = x - J^(-1) * f(x)
1872 m_pressureArea->GetJacobianInverse(invJ, Au, uu, beta, A_0, alpha,
1873 "Interface");
1874
1875 Multiply(dx, invJ, f);
1876
1877 // Update solution: x_new = x_old - dx
1878 for (size_t i = 0; i < 2; ++i)
1879 {
1880 uu[i] -= dx[i];
1881 Au[i] -= dx[i + 2];
1882 }
1883
1884 // Check if the error of the solution is smaller than Tol.
1885 if (Dot(dx, dx) < Tol)
1886 {
1887 proceed = 0;
1888 }
1889 time_InterfaceRiemann.Stop();
1890 time_InterfaceRiemann.AccumulateRegion(
1891 "PulseWaveSystem::InterfaceRiemann", 2);
1892 }
1893
1894 if (iter >= MAX_ITER)
1895 {
1896 ASSERTL0(false, "Riemann solver for Interface did not converge");
1897 }
1898}

References Nektar::LibUtilities::Timer::AccumulateRegion(), ASSERTL0, Nektar::LibUtilities::beta, Nektar::Dot(), m_pressureArea, m_rho, Nektar::Multiply(), Nektar::LibUtilities::Timer::Start(), and Nektar::LibUtilities::Timer::Stop().

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

1707{
1708
1709 NekDouble rho = m_rho;
1710 Array<OneD, NekDouble> W(3);
1711 Array<OneD, NekDouble> W_Au(3);
1712 Array<OneD, NekDouble> P_Au(3);
1713 NekMatrix<NekDouble> invJ(6, 6);
1714 NekVector<NekDouble> f(6);
1715 NekVector<NekDouble> dx(6);
1716
1717 int proceed = 1;
1718 int iter = 0;
1719 int MAX_ITER = 15;
1720
1721 // Backward characteristic
1722 m_pressureArea->GetW2(W[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1723
1724 // Forward characteristics
1725 for (size_t i = 1; i < 3; ++i)
1726 {
1727 m_pressureArea->GetW1(W[i], uu[i], beta[i], Au[i], A_0[i], alpha[i]);
1728 }
1729
1730 // Tolerances for the algorithm
1731 NekDouble Tol = 1.0E-10;
1732
1733 // Newton Iteration
1734 while ((proceed) && (iter < MAX_ITER))
1735 {
1736 LibUtilities::Timer time_MergingRiemann;
1737 time_MergingRiemann.Start();
1738 iter += 1;
1739
1740 /*
1741 * We solve the six constraint equations via a multivariate Newton
1742 * iteration. Equations are:
1743 * 1. Backward characteristic: W2(A_R, U_R) = W1(Au_R,
1744 * Uu_R)
1745 * 2. Forward characteristic 1: W1(A_L1, U_L1) = W1(Au_L1,
1746 * Uu_R1)
1747 * 3. Forward characteristic 2: W1(A_L2, U_L2) = W1(Au_L2,
1748 * Uu_L2)
1749 * 4. Conservation of mass: Au_R * Uu_R = Au_L1 *
1750 * Uu_L1 + Au_L2 * Uu_L2
1751 * 5. Continuity of total pressure 1: rho * Uu_R * Uu_R / 2 +
1752 * p(Au_R) = rho * Uu_L1 * Uu_L1 / 2 + p(Au_L1)
1753 * 6. Continuity of total pressure 2: rho * Uu_R * Uu_R / 2 +
1754 * p(Au_R) = rho * Uu_L2 * Uu_L2 / 2 + p(Au_L2)
1755 */
1756 for (size_t i = 0; i < 3; ++i)
1757 {
1758 m_pressureArea->GetPressure(P_Au[i], beta[i], Au[i], A_0[i], 0, 0,
1759 alpha[i]);
1760 }
1761
1762 m_pressureArea->GetW2(W_Au[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1763 m_pressureArea->GetW1(W_Au[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
1764 m_pressureArea->GetW1(W_Au[2], uu[2], beta[2], Au[2], A_0[2], alpha[2]);
1765
1766 // Constraint equations set to zero
1767 f[0] = W_Au[0] - W[0];
1768 f[1] = W_Au[1] - W[1];
1769 f[2] = W_Au[2] - W[2];
1770 f[3] = Au[0] * uu[0] - Au[1] * uu[1] - Au[2] * uu[2];
1771 f[4] = uu[0] * uu[0] + 2 * P_Au[0] / rho - uu[1] * uu[1] -
1772 2 * P_Au[1] / rho;
1773 f[5] = uu[0] * uu[0] + 2 * P_Au[0] / rho - uu[2] * uu[2] -
1774 2 * P_Au[2] / rho;
1775
1776 // Inverse Jacobian for x + dx = x - J^(-1) * f(x)
1777 m_pressureArea->GetJacobianInverse(invJ, Au, uu, beta, A_0, alpha,
1778 "Merge");
1779
1780 Multiply(dx, invJ, f);
1781
1782 // Update the solution: x_new = x_old - dx
1783 for (size_t i = 0; i < 3; ++i)
1784 {
1785 uu[i] -= dx[i];
1786 Au[i] -= dx[i + 3];
1787 }
1788
1789 // Check if the error of the solution is smaller than Tol
1790 if (Dot(dx, dx) < Tol)
1791 {
1792 proceed = 0;
1793 }
1794
1795 // Check if solver converges
1796 if (iter >= MAX_ITER)
1797 {
1798 ASSERTL0(false, "Riemann solver for Merging Flow did not converge");
1799 }
1800 time_MergingRiemann.Stop();
1801 time_MergingRiemann.AccumulateRegion("PulseWaveSystem::MergingRiemann",
1802 2);
1803 }
1804}

References Nektar::LibUtilities::Timer::AccumulateRegion(), ASSERTL0, Nektar::LibUtilities::beta, Nektar::Dot(), m_pressureArea, m_rho, Nektar::Multiply(), Nektar::LibUtilities::Timer::Start(), and Nektar::LibUtilities::Timer::Stop().

Referenced by EnforceInterfaceConditions().

◆ SetUpDomainInterfaceBCs()

void Nektar::PulseWaveSystem::SetUpDomainInterfaceBCs ( SpatialDomains::BoundaryConditions AllBcs,
std::map< int, LibUtilities::CommSharedPtr > &  domComm 
)
private

Definition at line 1074 of file PulseWaveSystem.cpp.

1077{
1078
1080 Allbcs.GetBoundaryRegions();
1081
1082 /* Loop over domain and find out if we have any undefined
1083 * boundary conditions representing interfaces. If so make a
1084 * map based around vid and storing the domains that are
1085 * part of interfaces.
1086 */
1087 int numNewBc = 1;
1088 for (auto &d : m_domOrder)
1089 {
1090
1091 std::map<int, SpatialDomains::GeometrySharedPtr> domvids;
1092
1093 // Loop over each domain and find which vids are only touched
1094 // by one element
1095 for (auto &compIt : m_domain[d])
1096 {
1097 for (size_t j = 0; j < compIt.second->m_geomVec.size(); ++j)
1098 {
1099
1100 // get hold of vids of each segment
1101 for (size_t p = 0; p < 2; ++p)
1102 {
1104 compIt.second->m_geomVec[j]->GetVertex(p);
1105 int vid = vert->GetGlobalID();
1106
1107 // if vid has already been added remove it otherwise add
1108 // into set
1109 if (domvids.count(vid))
1110 {
1111 domvids.erase(vid);
1112 }
1113 else
1114 {
1115 domvids[vid] = vert;
1116 }
1117 }
1118 }
1119 }
1120
1121 ASSERTL1(domvids.size() == 2, "Failed to find two end points "
1122 "of a domain (domvids = " +
1123 std::to_string(domvids.size()) + ")");
1124
1125 size_t nprocs = domComm[d]->GetSize();
1126
1127 if (nprocs > 1) // Remove parallel interfaces
1128 {
1129 size_t rank = domComm[d]->GetRank();
1130 Array<OneD, int> nvids(nprocs, 0);
1131 nvids[rank] = domvids.size();
1132 domComm[d]->AllReduce(nvids, LibUtilities::ReduceSum);
1133
1134 size_t totvids = Vmath::Vsum(nprocs, nvids, 1);
1135
1136 Array<OneD, int> locids(totvids, -1);
1137
1138 size_t cnt = 0;
1139 for (size_t i = 0; i < rank; ++i)
1140 {
1141 cnt += nvids[i];
1142 }
1143
1144 for (auto &i : domvids)
1145 {
1146 locids[cnt++] = i.first;
1147 }
1148
1149 // collect lcoal vids on this domain on all processor
1150 domComm[d]->AllReduce(locids, LibUtilities::ReduceMax);
1151
1152 std::set<int> chkvids;
1153 for (size_t i = 0; i < totvids; ++i)
1154 {
1155 if (chkvids.count(locids[i]))
1156 {
1157 // if this id is on local domain then remove
1158 if (domvids.count(locids[i]))
1159 {
1160 domvids.erase(locids[i]);
1161 }
1162 }
1163 else
1164 {
1165 chkvids.insert(locids[i]);
1166 }
1167 }
1168 }
1169
1170 // Finally check if remaining domvids are already declared, if
1171 // not add them
1172
1173 // Erase from domvids any existing BCs
1174 for (auto &it : bregions)
1175 {
1176 for (auto &bregionIt : *it.second)
1177 {
1178 // can assume that all regions only contain one point in 1D
1179 // Really do not need loop above
1180 int id = bregionIt.second->m_geomVec[0]->GetGlobalID();
1181
1182 if (domvids.count(id))
1183 {
1184 domvids.erase(id);
1185 }
1186 }
1187 }
1188
1189 // Finally add interface conditions to Allbcs
1190 std::vector<std::string> variables = m_session->GetVariables();
1191
1192 // determine the maxmimum bregion index
1193 int maxbregind = 0;
1194 for (auto &b : bregions)
1195 {
1196 maxbregind = (maxbregind > b.first) ? maxbregind : b.first;
1197 }
1198
1199 for (auto &b : domvids)
1200 {
1202 MemoryManager<
1203 SpatialDomains::BoundaryRegion>::AllocateSharedPtr());
1204
1205 // Set up Composite (GemetryVector) to contain vertex and put
1206 // into bRegion
1209 gvec->m_geomVec.push_back(b.second);
1210 (*breg)[b.first] = gvec;
1211
1212 Allbcs.AddBoundaryRegions(maxbregind + numNewBc, breg);
1213
1215 MemoryManager<
1216 SpatialDomains::BoundaryConditionMap>::AllocateSharedPtr();
1217
1218 std::string userDefined = "Interface";
1219 // Set up just boundary condition for this variable.
1220 SpatialDomains::BoundaryConditionShPtr DirichletInterface(
1221 MemoryManager<SpatialDomains::DirichletBoundaryCondition>::
1222 AllocateSharedPtr(m_session, "0", userDefined));
1223
1224 for (size_t i = 0; i < variables.size(); ++i)
1225 {
1226 (*bCondition)[variables[i]] = DirichletInterface;
1227 }
1228
1229 Allbcs.AddBoundaryConditions(maxbregind + numNewBc, bCondition);
1230 ++numNewBc;
1231 }
1232 }
1233}
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
std::map< int, CompositeSharedPtr > BoundaryRegion
Definition: Conditions.h:208
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:211
std::shared_ptr< BoundaryRegion > BoundaryRegionShPtr
Definition: Conditions.h:209
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:213
std::map< std::string, BoundaryConditionShPtr > BoundaryConditionMap
Definition: Conditions.h:218
std::shared_ptr< Composite > CompositeSharedPtr
Definition: MeshGraph.h:135
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:219
std::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:51

References Nektar::SpatialDomains::BoundaryConditions::AddBoundaryConditions(), Nektar::SpatialDomains::BoundaryConditions::AddBoundaryRegions(), Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Nektar::UnitTests::d(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), m_domain, m_domOrder, Nektar::SolverUtils::EquationSystem::m_session, CellMLToNektar.cellml_metadata::p, Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceSum, and Vmath::Vsum().

Referenced by v_InitObject().

◆ SetUpDomainInterfaces()

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

Definition at line 799 of file PulseWaveSystem.cpp.

800{
801 std::map<int, std::vector<InterfacePointShPtr>> VidToDomain;
802
803 // First set up domain to determine how many vertices meet at a
804 // vertex. This allows us to
805
806 /* Loop over domain and find out if we have any undefined
807 * boundary conditions representing interfaces. If so make a
808 * map based around vid and storing the domains that are
809 * part of interfaces.
810 */
811 for (size_t omega = 0; omega < m_nDomains; ++omega)
812 {
813 size_t vesselID = omega * m_nVariables;
814
815 for (size_t i = 0; i < (m_vessels[vesselID]->GetBndConditions()).size();
816 ++i)
817 {
818 if (m_vessels[vesselID]->GetBndConditions()[i]->GetUserDefined() ==
819 "Interface")
820 {
821 // Get Vid of interface
822 int vid = m_vessels[vesselID]
823 ->UpdateBndCondExpansion(i)
824 ->GetExp(0)
825 ->GetGeom()
826 ->GetVid(0);
827
829 m_vessels[vesselID]->GetTrace();
831
832 bool finish = false;
833 // find which elmt, the local vertex and the data
834 // offset of point
835 size_t nExp = m_vessels[vesselID]->GetExpSize();
836 for (size_t n = 0; n < nExp; ++n)
837 {
838 for (size_t p = 0; p < 2; ++p)
839 {
840
841 if (m_vessels[vesselID]
842 ->GetTraceMap()
843 ->GetElmtToTrace()[n][p]
844 ->as<LocalRegions::Expansion>()
845 ->GetGeom()
846 ->GetVid(0) == vid)
847 {
848 int eid = m_vessels[vesselID]
849 ->GetTraceMap()
850 ->GetElmtToTrace()[n][p]
851 ->GetElmtId();
852
853 int tid = m_vessels[vesselID]
854 ->GetTrace()
855 ->GetCoeff_Offset(eid);
856
857 Ipt = MemoryManager<
858 InterfacePoint>::AllocateSharedPtr(vid, omega,
859 n, p, tid,
860 i);
861
862 finish = true;
863 break;
864 }
865 }
866 if (finish == true)
867 {
868 break;
869 }
870 }
871
872 VidToDomain[vid].push_back(Ipt);
873 }
874 }
875 }
876
877 // Set up comms for parallel cases
878 std::map<int, int> domId;
879
880 size_t cnt = 0;
881 for (auto &d : m_domain)
882 {
883 domId[cnt] = d.first;
884 ++cnt;
885 }
886 ASSERTL1(domId.size() == m_nDomains, "Number of domains do not match");
887
888 Gs::gs_data *gscomm;
889 int nvid2dom = VidToDomain.size();
890 Array<OneD, long> tmp(nvid2dom);
891 Array<OneD, NekDouble> nvid(nvid2dom, 0.0);
892
893 cnt = 0;
894 for (auto &v : VidToDomain)
895 {
896 tmp[cnt] = v.first + 1;
897 for (size_t i = 0; i < v.second.size(); ++i)
898 {
899 nvid[cnt] += 1.0;
900 }
901 cnt++;
902 }
903
904 const bool verbose = m_session->DefinesCmdLineArgument("verbose");
905
906 gscomm = Gs::Init(tmp, m_comm->GetRowComm(), verbose);
907 Gs::Gather(nvid, Gs::gs_add, gscomm);
908
909 // Loop over domains and identify how many vessels at a
910 // bifurcation use the beginning of the element at the
911 // junction. This distinguishes a bifurcation and merging junction
912 // since we implicitly assume elements are ordered left ot right
913 // and so at a bifurcation we have two elements meeting this point
914 // at the first vertex
915 Array<OneD, NekDouble> nbeg(nvid2dom, 0.0);
916 cnt = 0;
917 for (auto &v : VidToDomain)
918 {
919 if (nvid[cnt] == 3.0)
920 {
921 for (size_t i = 0; i < v.second.size(); ++i)
922 {
923 // for bifurcations and merging junctions store how
924 // many points at interface are at the beginning or
925 // end
926 if (v.second[i]->m_elmtVert == 0)
927 {
928 nbeg[cnt] += 1.0;
929 }
930 }
931 }
932 ++cnt;
933 }
934 Gs::Gather(nbeg, Gs::gs_add, gscomm);
935
936 // Finally determine the maximum domain id between the two vessels
937 // at an interfaces or the maximum domain id of the daughter
938 // vessels for a bifurcation or the two parent id's for a merging
939 // junction
940 Array<OneD, NekDouble> dom(VidToDomain.size(), 0.0);
941 cnt = 0;
942 for (auto &v : VidToDomain)
943 {
944 if (nvid[cnt] == 2.0)
945 {
946 // store the domain id of the two domains
947 for (size_t i = 0; i < v.second.size(); ++i)
948 {
949 dom[cnt] =
950 std::max(dom[cnt], (NekDouble)domId[v.second[i]->m_domain]);
951 }
952 }
953 else if (nvid[cnt] == 3.0)
954 {
955 // store the domain id of the two daughter vessels if a
956 // bifurcation or the two parent vessels if a merging
957 // junction
958
959 int val = (nbeg[cnt] == 2.0) ? 0 : 1;
960
961 for (size_t i = 0; i < v.second.size(); ++i)
962 {
963 if (v.second[i]->m_elmtVert == val)
964 {
965 dom[cnt] = std::max(
966 dom[cnt], (NekDouble)domId[v.second[i]->m_domain]);
967 }
968 }
969 }
970 ++cnt;
971 }
972 Gs::Gather(dom, Gs::gs_max, gscomm);
973
974 // loop over map and set up Interface information;
975 int v = 0;
976 for (auto &iter : VidToDomain)
977 {
978 ASSERTL1(nvid[v] != 1.0, "Found an interface wth only "
979 "one domain which should not happen");
980
981 if (nvid[v] == 2.0) // Vessel jump interface
982 {
983 for (size_t i = 0; i < iter.second.size(); ++i)
984 {
985 if (domId[iter.second[i]->m_domain] == dom[v])
986 {
987 iter.second[i]->m_riemannOrd = 1;
988 }
989 else
990 {
991 iter.second[i]->m_riemannOrd = 0;
992 }
993 }
994 m_vesselIntfcs.push_back(iter.second);
995 }
996 else if (nvid[v] == 3.0) // Bifurcation or Merging junction.
997 {
998 // Set up Bifurcation information
999 int val = (nbeg[v] == 2.0) ? 1 : 0;
1000
1001 for (size_t i = 0; i < iter.second.size(); ++i)
1002 {
1003 if (iter.second[i]->m_elmtVert == val)
1004 {
1005 iter.second[i]->m_riemannOrd = 0;
1006 }
1007 else
1008 {
1009 if (domId[iter.second[i]->m_domain] == dom[v])
1010 {
1011 iter.second[i]->m_riemannOrd = 2;
1012 }
1013 else
1014 {
1015 iter.second[i]->m_riemannOrd = 1;
1016 }
1017 }
1018 }
1019
1020 if (nbeg[v] == 2.0)
1021 {
1022 m_bifurcations.push_back(iter.second);
1023 }
1024 else
1025 {
1026 m_mergingJcts.push_back(iter.second);
1027 }
1028 }
1029 else
1030 {
1031 ASSERTL0(false, "Unknown junction type");
1032 }
1033 ++v; // incrementing loop over vertices
1034 }
1035
1036 // finally set up a communicator for interface data collection
1037 Array<OneD, long> tmp1(3 * nvid2dom);
1038 if (nvid2dom)
1039 {
1040 Vmath::Zero(3 * nvid2dom, tmp1, 1);
1041 }
1042
1043 cnt = 0;
1044 for (size_t n = 0; n < m_bifurcations.size(); ++n, ++cnt)
1045 {
1046 size_t vid = m_bifurcations[n][0]->m_vid;
1047 for (size_t i = 0; i < 3; ++i)
1048 {
1049 tmp1[3 * cnt + i] = (NekDouble)(3 * vid + i + 1);
1050 }
1051 }
1052
1053 for (size_t n = 0; n < m_mergingJcts.size(); ++n, ++cnt)
1054 {
1055 size_t vid = m_mergingJcts[n][0]->m_vid;
1056 for (size_t i = 0; i < 3; ++i)
1057 {
1058 tmp1[3 * cnt + i] = (NekDouble)(3 * vid + i + 1);
1059 }
1060 }
1061
1062 for (size_t n = 0; n < m_vesselIntfcs.size(); ++n, ++cnt)
1063 {
1064 size_t vid = m_vesselIntfcs[n][0]->m_vid;
1065 for (size_t i = 0; i < 3; ++i)
1066 {
1067 tmp1[3 * cnt + i] = (NekDouble)(3 * vid + i + 1);
1068 }
1069 }
1070
1071 m_intComm = Gs::Init(tmp1, m_comm->GetRowComm(), verbose);
1072}
static gs_data * Init(const Nektar::Array< OneD, long > &pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:190
@ gs_max
Definition: GsLib.hpp:63
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< InterfacePoint > InterfacePointShPtr
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273

References ASSERTL0, ASSERTL1, Nektar::UnitTests::d(), Gs::Gather(), Gs::gs_add, Gs::gs_max, Gs::Init(), m_bifurcations, Nektar::SolverUtils::EquationSystem::m_comm, m_domain, m_intComm, m_mergingJcts, m_nDomains, m_nVariables, Nektar::SolverUtils::EquationSystem::m_session, m_vesselIntfcs, m_vessels, CellMLToNektar.cellml_metadata::p, and Vmath::Zero().

Referenced by v_InitObject().

◆ UpdateVessels()

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

Definition at line 85 of file PulseWaveSystem.h.

86 {
87 return m_vessels;
88 }

References m_vessels.

◆ v_DoInitialise()

void Nektar::PulseWaveSystem::v_DoInitialise ( bool  dumpInitialConditions = false)
overrideprotectedvirtual

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

1244{
1245 if (m_session->GetComm()->GetRank() == 0)
1246 {
1247 std::cout << "Initial Conditions: " << std::endl;
1248 }
1249
1250 /* Loop over all subdomains to initialize all with the Initial
1251 * Conditions read from the inputfile*/
1252 int omega = 0;
1253 for (auto &d : m_domOrder)
1254 {
1255 for (size_t i = 0; i < m_nVariables; ++i)
1256 {
1257 m_fields[i] = m_vessels[m_nVariables * omega + i];
1258 }
1259
1260 if (m_session->GetComm()->GetRank() == 0)
1261 {
1262 std::cout << "Subdomain: " << omega << std::endl;
1263 }
1264
1265 SetInitialConditions(0.0, false, d);
1266 omega++;
1267 }
1268
1269 // Reset variables to first vessels
1270 for (size_t i = 0; i < m_nVariables; ++i)
1271 {
1272 m_fields[i] = m_vessels[i];
1273 }
1274
1275 // Update filters.
1276 for (auto &d : m_domainToFilterIDs)
1277 {
1278 Array<OneD, MultiRegions::ExpListSharedPtr> filterFields(m_nVariables);
1279
1280 int vesselID = m_filterToVesselID[d.second[0]];
1281 for (size_t j = 0; j < m_nVariables; ++j)
1282 {
1283 filterFields[j] = m_vessels[vesselID + j];
1284 }
1285
1286 // loop over filters in domain
1287 for (auto &f : d.second)
1288 {
1289 m_filters[f].second->Update(filterFields, m_time);
1290 }
1291 }
1292}
std::map< int, int > m_filterToVesselID
std::map< int, std::vector< int > > m_domainToFilterIDs
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters

References Nektar::UnitTests::d(), m_domainToFilterIDs, m_domOrder, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::UnsteadySystem::m_filters, m_filterToVesselID, m_nVariables, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_time, m_vessels, and Nektar::SolverUtils::EquationSystem::SetInitialConditions().

◆ v_DoSolve()

void Nektar::PulseWaveSystem::v_DoSolve ( void  )
overrideprotectedvirtual

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

if(m_session->GetComm()->GetRank() == 0) { cout << "Time-integration timing: " << IntegrationTime << " s" << endl << endl; }

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1311 of file PulseWaveSystem.cpp.

1312{
1313 size_t i;
1314 int n;
1315 int nchk = 1;
1316
1317 Array<OneD, Array<OneD, NekDouble>> fields(m_nVariables);
1318
1319 for (i = 0; i < m_nVariables; ++i)
1320 {
1321 fields[i] = m_vessels[i]->UpdatePhys();
1322 m_fields[i]->SetPhysState(false);
1323 }
1324
1325 m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
1326
1327 // Time loop
1328 for (n = 0; n < m_steps; ++n)
1329 {
1330 LibUtilities::Timer time_v_IntStep;
1331 time_v_IntStep.Start();
1332 fields = m_intScheme->TimeIntegrate(n, m_timestep);
1333 m_time += m_timestep;
1334 time_v_IntStep.Stop();
1335 time_v_IntStep.AccumulateRegion("PulseWaveSystem::TimeIntegrationStep",
1336 0);
1337 // Write out status information.
1338 if (m_infosteps && !((n + 1) % m_infosteps) &&
1339 m_session->GetComm()->GetRank() == 0)
1340 {
1341 std::cout << "Steps: " << n + 1 << "\t Time: " << m_time
1342 << "\t Time-step: " << m_timestep << "\t" << std::endl;
1343 }
1344
1345 // Copy Array To Vessel Phys Fields
1346 for (size_t i = 0; i < m_nVariables; ++i)
1347 {
1348 Vmath::Vcopy(fields[i].size(), fields[i], 1,
1349 m_vessels[i]->UpdatePhys(), 1);
1350 }
1351
1352 // Transform data if needed
1353 if (m_checksteps && !((n + 1) % m_checksteps))
1354 {
1355 for (i = 0; i < m_nVariables; ++i)
1356 {
1357 for (size_t omega = 0; omega < m_nDomains; omega++)
1358 {
1359 unsigned int offset = omega * m_nVariables + i;
1360 m_vessels[offset]->FwdTrans(
1361 fields[i] + m_fieldPhysOffset[omega],
1362 m_vessels[offset]->UpdateCoeffs());
1363 }
1364 }
1365 CheckPoint_Output(nchk++);
1366 }
1367
1368 // Update filters.
1369 for (auto &d : m_domainToFilterIDs)
1370 {
1371 Array<OneD, MultiRegions::ExpListSharedPtr> filterFields(
1372 m_nVariables);
1373 int vesselID = m_filterToVesselID[d.second[0]];
1374 for (size_t j = 0; j < m_nVariables; ++j)
1375 {
1376 filterFields[j] = m_vessels[vesselID + j];
1377 }
1378
1379 // loop over filters in domain
1380 for (auto &f : d.second)
1381 {
1382 m_filters[f].second->Update(filterFields, m_time);
1383 }
1384 }
1385
1386 } // end of timeintegration
1387
1388 /** if(m_session->GetComm()->GetRank() == 0)
1389 {
1390 cout << "Time-integration timing: " << IntegrationTime << " s" <<
1391 endl
1392 << endl;
1393 } **/
1394
1395 // Finalise filters.
1396 for (auto &d : m_domainToFilterIDs)
1397 {
1398 Array<OneD, MultiRegions::ExpListSharedPtr> filterFields(m_nVariables);
1399 int vesselID = m_filterToVesselID[d.second[0]];
1400 for (size_t j = 0; j < m_nVariables; ++j)
1401 {
1402 filterFields[j] = m_vessels[vesselID + j];
1403 }
1404
1405 // loop over filters in domain
1406 for (auto &f : d.second)
1407 {
1408 m_filters[f].second->Finalise(filterFields, m_time);
1409 }
1410 }
1411}
void CheckPoint_Output(const int n)
NekDouble m_timestep
Time step size.
int m_infosteps
Number of time steps between outputting status information.
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.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825

References Nektar::LibUtilities::Timer::AccumulateRegion(), CheckPoint_Output(), Nektar::UnitTests::d(), Nektar::SolverUtils::EquationSystem::m_checksteps, m_domainToFilterIDs, m_fieldPhysOffset, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::UnsteadySystem::m_filters, m_filterToVesselID, Nektar::SolverUtils::EquationSystem::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(), and Vmath::Vcopy().

◆ v_InitObject()

void Nektar::PulseWaveSystem::v_InitObject ( bool  DeclareField = false)
overrideprotectedvirtual

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.

Definition at line 80 of file PulseWaveSystem.cpp.

81{
82 // Initialise base class
83 UnsteadySystem::v_InitObject(DeclareField);
84
85 // Read the geometry and the expansion information
86 m_domain = m_graph->GetDomain();
87 m_nDomains = m_domain.size();
88
89 // Determine projectiontype
90 ASSERTL0(m_session->MatchSolverInfo("Projection", "DisContinuous"),
91 "Pulse solver only set up for Discontinuous projections");
94 m_graph->GetMeshDimension() == 1,
95 "Pulse wave solver only set up for expansion dimension equal to 1");
96
97 size_t i;
98 m_nVariables = m_session->GetVariables().size();
99
100 m_fields = Array<OneD, MultiRegions::ExpListSharedPtr>(m_nVariables);
101 m_vessels =
102 Array<OneD, MultiRegions::ExpListSharedPtr>(m_nVariables * m_nDomains);
103
104 /* If the PressureArea property is not specified, default to the Beta
105 * law; it's the most well known and this way previous examples that did
106 * not specify the tube law will still work.
107 */
108 if (m_session->DefinesSolverInfo("PressureArea"))
109 {
111 m_session->GetSolverInfo("PressureArea"), m_vessels, m_session);
112 }
113 else
114 {
116 "Beta", m_vessels, m_session);
117 }
118
119 m_pressure = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
120
121 m_PWV = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
122 m_W1 = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
123 m_W2 = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
124
125 SpatialDomains::BoundaryConditions Allbcs(m_session, m_graph);
126
127 // get parallel communicators as required
128 std::map<int, LibUtilities::CommSharedPtr> domComm;
129 GetCommArray(domComm);
130
131 SetUpDomainInterfaceBCs(Allbcs, domComm);
132
133 // set up a list of the domain to filter ids
134 auto filterMap = m_session->GetFilters();
135 unsigned int fcnt = 0;
136 for (auto &x : filterMap)
137 {
138 ASSERTL0(x.domain != -1, "Filter needs DOMAIN attribute when "
139 "applied to PulseWaveSolver.")
140 m_domainToFilterIDs[x.domain].push_back(fcnt);
141 ++fcnt;
142 }
143
144 // Set up domains and put geometry to be only one space dimension.
145 size_t cnt = 0;
146 // currently I believe we have to set to one space dimension on
147 // constuction since the advection terms are set to differentiate
148 // with respect to x (or the lcoal ordinate) not as an arc length
149 // derivative. So will not get consistent answers wtihout this
150 bool SetToOneSpaceDimension = true;
151
152 for (auto &d : m_domOrder)
153 {
154 // need to check if we have filters on this domain and if so
155 // set them up without SetToOneSpaceDimension if so we can
156 // initialisatio them (i.e. history points)
157 if (m_domainToFilterIDs.count(d))
158 {
159 // make temporary expansion list to not settig to one space
160 // dimension
161 Array<OneD, MultiRegions::ExpListSharedPtr> filterFields(
163 for (size_t j = 0; j < m_nVariables; ++j)
164 {
167 m_session->GetVariable(j), domComm[d],
168 false);
169 }
170
171 // intiialise filters with temporary expansion lsits
172 for (auto &x : m_domainToFilterIDs[d])
173 {
174 m_filters[x].second->SetUpdateOnInitialise(false);
175 m_filters[x].second->Initialise(filterFields, m_time);
176 m_filters[x].second->SetUpdateOnInitialise(true);
177
178 // save vessel start id for update
179 m_filterToVesselID[x] = cnt;
180 }
181 }
182
183 for (size_t j = 0; j < m_nVariables; ++j)
184 {
185 m_vessels[cnt++] =
187 m_session, m_graph, m_domain[d], Allbcs,
188 m_session->GetVariable(j), domComm[d],
189 SetToOneSpaceDimension);
190 }
191 }
192
193 // Reset coeff and phys space to be continuous over all domains
194 int totcoeffs = 0;
195 int totphys = 0;
196 m_fieldPhysOffset = Array<OneD, int>(m_nDomains + 1, 0);
197 for (i = 0; i < m_nDomains; ++i)
198 {
199 totcoeffs += m_vessels[i * m_nVariables]->GetNcoeffs();
200
201 m_fieldPhysOffset[i] = totphys;
202 totphys += m_vessels[i * m_nVariables]->GetTotPoints();
203 }
204 m_fieldPhysOffset[m_nDomains] = totphys;
205
206 for (size_t n = 0; n < m_nVariables; ++n)
207 {
208 Array<OneD, NekDouble> coeffs(totcoeffs, 0.0);
209 Array<OneD, NekDouble> phys(totphys, 0.0);
210 Array<OneD, NekDouble> tmpcoeffs, tmpphys;
211
212 m_vessels[n]->SetCoeffsArray(coeffs);
213 m_vessels[n]->SetPhysArray(phys);
214
215 int cnt = m_vessels[n]->GetNcoeffs();
216 int cnt1 = m_vessels[n]->GetTotPoints();
217
218 for (i = 1; i < m_nDomains; ++i)
219 {
220 m_vessels[i * m_nVariables + n]->SetCoeffsArray(tmpcoeffs =
221 coeffs + cnt);
222 m_vessels[i * m_nVariables + n]->SetPhysArray(tmpphys =
223 phys + cnt1);
224 cnt += m_vessels[i * m_nVariables + n]->GetNcoeffs();
225 cnt1 += m_vessels[i * m_nVariables + n]->GetTotPoints();
226 }
227 }
228
229 // If Discontinuous Galerkin determine upwinding method to use
230 for (size_t i = 0; i < SIZE_UpwindTypePulse; ++i)
231 {
232 bool match;
233 m_session->MatchSolverInfo("UPWINDTYPEPULSE", UpwindTypeMapPulse[i],
234 match, false);
235 if (match)
236 {
238 break;
239 }
240 }
241
242 // Load blood density and external pressure
243 m_session->LoadParameter("rho", m_rho, 0.5);
244 m_session->LoadParameter("pext", m_pext, 0.0);
245
246 int nq = 0;
247 /*
248 * Gets the Material Properties of each arterial segment
249 * specified in the inputfile from section MaterialProperties
250 * Also gets the Area at static equilibrium A_0 specified in the
251 * inputfile.
252 *
253 * Having found these points also extract the values at the
254 * trace points and the normal direction consistent with the
255 * left adjacent definition of Fwd and Bwd
256 */
257 m_beta = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
258 m_beta_trace = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
259 m_gamma = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
260 m_alpha = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
261 m_alpha_trace = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
262 m_A_0 = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
263 m_A_0_trace = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
264 m_trace_fwd_normal = Array<OneD, Array<OneD, NekDouble>>(m_nDomains);
265
266 int omega = 0;
267 for (auto &d : m_domOrder)
268 {
269 nq = m_vessels[2 * omega]->GetNpoints();
270 m_fields[0] = m_vessels[2 * omega];
271
272 m_beta[omega] = Array<OneD, NekDouble>(nq);
273 GetFunction("MaterialProperties")
274 ->Evaluate("beta", m_beta[omega], m_time, d);
275
276 // If the input file doesn't specify viscoelasticity, set it to zero.
277 m_gamma[omega] = Array<OneD, NekDouble>(nq);
278
279 if (m_session->DefinesFunction("Viscoelasticity"))
280 {
281 GetFunction("Viscoelasticity")
282 ->Evaluate("gamma", m_gamma[omega], m_time, d);
283 }
284 else
285 {
286 for (int j = 0; j < nq; ++j)
287 {
288 m_gamma[omega][j] = 0;
289 }
290 }
291
292 /* If the input file doesn't specify strain-stiffening, set it to
293 * 0.5 for elastic behaviour.
294 */
295 m_alpha[omega] = Array<OneD, NekDouble>(nq);
296
297 if (m_session->DefinesFunction("StrainStiffening"))
298 {
299 GetFunction("StrainStiffening")
300 ->Evaluate("alpha", m_alpha[omega], m_time, d);
301 }
302 else
303 {
304 for (int j = 0; j < nq; ++j)
305 {
306 m_alpha[omega][j] = 0.5;
307 }
308 }
309
310 m_A_0[omega] = Array<OneD, NekDouble>(nq);
311 GetFunction("A_0")->Evaluate("A_0", m_A_0[omega], m_time, d);
312
313 int nqTrace = GetTraceTotPoints();
314
315 m_beta_trace[omega] = Array<OneD, NekDouble>(nqTrace);
316 m_fields[0]->ExtractTracePhys(m_beta[omega], m_beta_trace[omega]);
317
318 m_A_0_trace[omega] = Array<OneD, NekDouble>(nqTrace);
319 m_fields[0]->ExtractTracePhys(m_A_0[omega], m_A_0_trace[omega]);
320
321 m_alpha_trace[omega] = Array<OneD, NekDouble>(nqTrace);
322 m_fields[0]->ExtractTracePhys(m_alpha[omega], m_alpha_trace[omega]);
323
324 m_trace_fwd_normal[omega] = Array<OneD, NekDouble>(nqTrace, 0.0);
325
326 MultiRegions::ExpListSharedPtr trace = m_fields[0]->GetTrace();
327 int nelmt_trace = trace->GetExpSize();
328
329 Array<OneD, Array<OneD, NekDouble>> normals(nelmt_trace);
330
331 for (int i = 0; i < nelmt_trace; ++i)
332 {
333 normals[i] = m_trace_fwd_normal[omega] + i;
334 }
335
336 // need to set to 1 for consistency since boundary
337 // conditions may not have coordim = 1
338 trace->GetExp(0)->GetGeom()->SetCoordim(1);
339
340 trace->GetNormals(normals);
341 omega++;
342 }
343
345}
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
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
void GetCommArray(std::map< int, LibUtilities::CommSharedPtr > &retval)
Set and retrn a series of communicators for each partition.
Array< OneD, Array< OneD, NekDouble > > m_alpha
Array< OneD, Array< OneD, NekDouble > > m_PWV
void SetUpDomainInterfaceBCs(SpatialDomains::BoundaryConditions &AllBcs, std::map< int, LibUtilities::CommSharedPtr > &domComm)
Array< OneD, Array< OneD, NekDouble > > m_beta
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
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::UnitTests::d(), Nektar::MultiRegions::eDiscontinuous, GetCommArray(), 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_domain, m_domainToFilterIDs, m_domOrder, m_fieldPhysOffset, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::UnsteadySystem::m_filters, m_filterToVesselID, 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, SetUpDomainInterfaceBCs(), SetUpDomainInterfaces(), Nektar::SIZE_UpwindTypePulse, Nektar::UpwindTypeMapPulse, and Nektar::SolverUtils::UnsteadySystem::v_InitObject().

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

◆ v_L2Error()

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

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 2025 of file PulseWaveSystem.cpp.

2028{
2029 NekDouble L2error = 0.0;
2030 NekDouble L2error_dom;
2031 NekDouble Vol = 0.0;
2032
2033 if (m_NumQuadPointsError == 0)
2034 {
2035 for (size_t omega = 0; omega < m_nDomains; omega++)
2036 {
2037 size_t vesselid = field + omega * m_nVariables;
2038
2039 // since exactsoln is passed for just the first field size we
2040 // need to reset it each domain
2041 Array<OneD, NekDouble> exactsoln(m_vessels[vesselid]->GetNpoints());
2043 EvaluateExactSolution(field, exactsoln, m_time);
2044
2045 if (m_vessels[vesselid]->GetPhysState() == false)
2046 {
2047 m_vessels[vesselid]->BwdTrans(
2048 m_vessels[vesselid]->GetCoeffs(),
2049 m_vessels[vesselid]->UpdatePhys());
2050 }
2051
2052 if (exactsoln.size())
2053 {
2054 L2error_dom = m_vessels[vesselid]->L2(
2055 m_vessels[vesselid]->GetPhys(), exactsoln);
2056 }
2057 else if (m_session->DefinesFunction("ExactSolution"))
2058 {
2059 Array<OneD, NekDouble> exactsoln(
2060 m_vessels[vesselid]->GetNpoints());
2061
2063 m_session->GetFunction("ExactSolution", field, omega);
2064 GetFunction("ExactSolution")
2065 ->Evaluate(m_session->GetVariable(field), exactsoln,
2066 m_time);
2067
2068 L2error_dom = m_vessels[vesselid]->L2(
2069 m_vessels[vesselid]->GetPhys(), exactsoln);
2070 }
2071 else
2072 {
2073 L2error_dom =
2074 m_vessels[vesselid]->L2(m_vessels[vesselid]->GetPhys());
2075 }
2076
2077 if (m_vessels[vesselid]->GetComm()->GetRank())
2078 {
2079 // ensure domains are only summed on local root of
2080 // domain
2081 L2error_dom = 0.0;
2082 }
2083
2084 L2error += L2error_dom * L2error_dom;
2085
2086 if (Normalised == true)
2087 {
2088 Array<OneD, NekDouble> one(m_vessels[vesselid]->GetNpoints(),
2089 1.0);
2090
2091 Vol += m_vessels[vesselid]->Integral(one);
2092 }
2093 }
2094 }
2095 else
2096 {
2097 ASSERTL0(false, "Not set up");
2098 }
2099
2100 m_comm->AllReduce(L2error, LibUtilities::ReduceSum);
2101
2102 if (Normalised == true)
2103 {
2104 m_comm->AllReduce(Vol, LibUtilities::ReduceSum);
2105
2106 L2error = sqrt(L2error / Vol);
2107 }
2108 else
2109 {
2110 L2error = sqrt(L2error);
2111 }
2112
2113 return L2error;
2114}
SOLVER_UTILS_EXPORT int GetNpoints()
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
SOLVER_UTILS_EXPORT void EvaluateExactSolution(int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Evaluates an exact solution.
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:125
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285

References ASSERTL0, Nektar::SolverUtils::EquationSystem::EvaluateExactSolution(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::GetFunction(), Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::SolverUtils::EquationSystem::m_comm, Nektar::SolverUtils::EquationSystem::m_fields, 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 > &  exact = NullNekDouble1DArray 
)
overrideprotectedvirtual

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

2124{
2125 NekDouble LinferrorDom, Linferror = -1.0;
2126
2127 for (size_t omega = 0; omega < m_nDomains; ++omega)
2128 {
2129 size_t vesselid = field + omega * m_nVariables;
2130
2131 // since exactsoln is passed for just the first field size we need
2132 // to reset it each domain
2133 Array<OneD, NekDouble> exactsoln(m_vessels[vesselid]->GetNpoints());
2135 EvaluateExactSolution(field, exactsoln, m_time);
2136
2137 if (m_NumQuadPointsError == 0)
2138 {
2139 if (m_vessels[vesselid]->GetPhysState() == false)
2140 {
2141 m_vessels[vesselid]->BwdTrans(
2142 m_vessels[vesselid]->GetCoeffs(),
2143 m_vessels[vesselid]->UpdatePhys());
2144 }
2145
2146 if (exactsoln.size())
2147 {
2148 LinferrorDom = m_vessels[vesselid]->Linf(
2149 m_vessels[vesselid]->GetPhys(), exactsoln);
2150 }
2151 else if (m_session->DefinesFunction("ExactSolution"))
2152 {
2153 Array<OneD, NekDouble> exactsoln(
2154 m_vessels[vesselid]->GetNpoints());
2155
2156 GetFunction("ExactSolution")
2157 ->Evaluate(m_session->GetVariable(field), exactsoln,
2158 m_time);
2159
2160 LinferrorDom = m_vessels[vesselid]->Linf(
2161 m_vessels[vesselid]->GetPhys(), exactsoln);
2162 }
2163 else
2164 {
2165 LinferrorDom = 0.0;
2166 }
2167
2168 Linferror = (Linferror > LinferrorDom) ? Linferror : LinferrorDom;
2169 }
2170 else
2171 {
2172 ASSERTL0(false, "ErrorExtraPoints not allowed for this solver");
2173 }
2174 }
2175 m_comm->AllReduce(Linferror, LibUtilities::ReduceMax);
2176 return Linferror;
2177}

References ASSERTL0, Nektar::SolverUtils::EquationSystem::EvaluateExactSolution(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::GetFunction(), Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::SolverUtils::EquationSystem::m_comm, Nektar::SolverUtils::EquationSystem::m_fields, m_nDomains, Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError, m_nVariables, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_time, m_vessels, and Nektar::LibUtilities::ReduceMax.

◆ v_Output()

void Nektar::PulseWaveSystem::v_Output ( void  )
overrideprotectedvirtual

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

1905{
1906 /**
1907 * Write the field data to file. The file is named according to the
1908 * session name with the extension .fld appended.
1909 */
1910 std::string outname = m_sessionName + ".fld";
1911
1912 WriteVessels(outname);
1913}

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

1932{
1933 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1934 std::vector<std::string> variables = m_session->GetVariables();
1935
1936 for (size_t n = 0; n < m_nDomains; ++n)
1937 {
1938 m_vessels[n * m_nVariables]->GetFieldDefinitions(FieldDef);
1939 }
1940
1941 std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
1942
1943 size_t nFieldDefPerDomain = FieldDef.size() / m_nDomains;
1944 size_t cnt;
1945
1946 // Copy Data into FieldData and set variable
1947 for (size_t n = 0; n < m_nDomains; ++n)
1948 {
1949 // Outputs area and velocity
1950 for (size_t j = 0; j < m_nVariables; ++j)
1951 {
1952 for (size_t i = 0; i < nFieldDefPerDomain; ++i)
1953 {
1954 cnt = n * nFieldDefPerDomain + i;
1955
1956 FieldDef[cnt]->m_fields.push_back(variables[j]);
1957
1958 m_vessels[n * m_nVariables]->AppendFieldData(
1959 FieldDef[cnt], FieldData[cnt],
1960 m_vessels[n * m_nVariables + j]->UpdateCoeffs());
1961 }
1962 }
1963
1964 // Outputs pressure
1965 Array<OneD, NekDouble> PFwd(m_vessels[n * m_nVariables]->GetNcoeffs());
1966
1967 m_vessels[n * m_nVariables]->FwdTransLocalElmt(m_pressure[n], PFwd);
1968
1969 for (size_t i = 0; i < nFieldDefPerDomain; ++i)
1970 {
1971 cnt = n * nFieldDefPerDomain + i;
1972
1973 FieldDef[cnt]->m_fields.push_back("P");
1974
1975 m_vessels[n * m_nVariables]->AppendFieldData(FieldDef[cnt],
1976 FieldData[cnt], PFwd);
1977 }
1978
1979 if (extraFields)
1980 {
1981 Array<OneD, NekDouble> PWVFwd(
1983 Array<OneD, NekDouble> W1Fwd(
1985 Array<OneD, NekDouble> W2Fwd(
1987
1988 m_vessels[n * m_nVariables]->FwdTransLocalElmt(m_PWV[n], PWVFwd);
1989 m_vessels[n * m_nVariables]->FwdTransLocalElmt(m_W1[n], W1Fwd);
1990 m_vessels[n * m_nVariables]->FwdTransLocalElmt(m_W2[n], W2Fwd);
1991
1992 for (size_t i = 0; i < nFieldDefPerDomain; ++i)
1993 {
1994 cnt = n * nFieldDefPerDomain + i;
1995
1996 FieldDef[cnt]->m_fields.push_back("c");
1997 FieldDef[cnt]->m_fields.push_back("W1");
1998 FieldDef[cnt]->m_fields.push_back("W2");
1999
2000 m_vessels[n * m_nVariables]->AppendFieldData(
2001 FieldDef[cnt], FieldData[cnt], PWVFwd);
2002 m_vessels[n * m_nVariables]->AppendFieldData(
2003 FieldDef[cnt], FieldData[cnt], W1Fwd);
2004 m_vessels[n * m_nVariables]->AppendFieldData(
2005 FieldDef[cnt], FieldData[cnt], W2Fwd);
2006 }
2007 }
2008 }
2009
2010 // Update time in field info if required
2011 if (m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
2012 {
2013 m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
2014 }
2015
2016 m_fld->Write(outname, FieldDef, FieldData, m_fieldMetaDataMap);
2017}
LibUtilities::FieldIOSharedPtr m_fld
Field input/output.
SOLVER_UTILS_EXPORT int GetNcoeffs()
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.

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

Referenced by CheckPoint_Output(), and v_Output().

Member Data Documentation

◆ extraFields

bool Nektar::PulseWaveSystem::extraFields = false
protected

Definition at line 124 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 130 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_C

NekDouble Nektar::PulseWaveSystem::m_C
protected

Definition at line 101 of file PulseWaveSystem.h.

◆ m_currentDomain

size_t Nektar::PulseWaveSystem::m_currentDomain
protected

◆ m_domain

std::map<int, SpatialDomains::CompositeMap> Nektar::PulseWaveSystem::m_domain
protected

◆ m_domainToFilterIDs

std::map<int, std::vector<int> > Nektar::PulseWaveSystem::m_domainToFilterIDs
protected

Definition at line 117 of file PulseWaveSystem.h.

Referenced by v_DoInitialise(), v_DoSolve(), and v_InitObject().

◆ m_domOrder

std::vector<int> Nektar::PulseWaveSystem::m_domOrder
protected

◆ m_fieldPhysOffset

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

Definition at line 97 of file PulseWaveSystem.h.

Referenced by FillDataFromInterfacePoint(), v_DoSolve(), and v_InitObject().

◆ m_filterToVesselID

std::map<int, int> Nektar::PulseWaveSystem::m_filterToVesselID
protected

Definition at line 118 of file PulseWaveSystem.h.

Referenced by v_DoInitialise(), v_DoSolve(), and v_InitObject().

◆ m_gamma

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

Definition at line 109 of file PulseWaveSystem.h.

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

◆ m_intComm

Gs::gs_data* Nektar::PulseWaveSystem::m_intComm
private

Communicator for interfaces.

Definition at line 196 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_mergingJcts

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

Definition at line 131 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_nDomains

size_t Nektar::PulseWaveSystem::m_nDomains
protected

◆ m_nVariables

size_t Nektar::PulseWaveSystem::m_nVariables
protected

◆ m_pext

NekDouble Nektar::PulseWaveSystem::m_pext
protected

Definition at line 99 of file PulseWaveSystem.h.

Referenced by v_InitObject().

◆ m_pout

NekDouble Nektar::PulseWaveSystem::m_pout
protected

Definition at line 103 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 102 of file PulseWaveSystem.h.

◆ m_trace_fwd_normal

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

Definition at line 112 of file PulseWaveSystem.h.

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

◆ m_upwindTypePulse

UpwindTypePulse Nektar::PulseWaveSystem::m_upwindTypePulse
protected

Definition at line 95 of file PulseWaveSystem.h.

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

◆ m_vesselIntfcs

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

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