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

virtual ~PulseWaveSystem ()
 Destructor. More...
 
int GetNdomains ()
 
Array< OneD, MultiRegions::ExpListSharedPtr > UpdateVessels (void)
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
virtual SOLVER_UTILS_EXPORT ~UnsteadySystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Calculate the larger time-step mantaining the problem stable. More...
 
SOLVER_UTILS_EXPORT void SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeSharedPtr & GetTimeIntegrationScheme ()
 Returns the time integration scheme. More...
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeOperators & GetTimeIntegrationSchemeOperators ()
 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 Array< OneD, NekDouble > ErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write checkpoint file of custom data fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write input fields to the given filename. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
 Input field data from the given file to multiple domains. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, std::vector< std::string > &fieldStr, Array< OneD, Array< OneD, NekDouble > > &coeffs)
 Output a field. Input field data into array from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, MultiRegions::ExpListSharedPtr &pField, std::string &pFieldName)
 Output a field. Input field data into ExpList from the given file. More...
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary. More...
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMap & UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetFinalTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, const Array< OneD, const NekDouble > &input)
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT 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...
 
SOLVER_UTILS_EXPORT bool ParallelInTime ()
 Check if solver use Parallel-in-Time. More...
 

Protected Member Functions

 PulseWaveSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises PulseWaveSystem class members. More...
 
virtual void v_InitObject (bool DeclareField=false) override
 
virtual void v_DoInitialise (bool dumpInitialConditions=false) override
 Sets up initial conditions. More...
 
virtual 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...
 
virtual void v_Output (void) override
 
void CheckPoint_Output (const int n)
 
virtual 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...
 
virtual 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...
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. More...
 
virtual 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...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true) override
 Sets up initial conditions. More...
 
virtual 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 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::ExpListSharedPtr > m_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::CompositeMap > m_domain
 
std::vector< int > m_domOrder
 
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::SessionFunctionSharedPtr > m_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
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, NekDouble > m_movingFrameVelsxyz
 Moving frame of reference velocities. More...
 
Array< OneD, NekDouble > m_movingFrameTheta
 Moving frame of reference angles with respect to the. More...
 
boost::numeric::ublas::matrix< NekDouble > m_movingFrameProjMat
 Projection matrix for transformation between inertial and moving. More...
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error. More...
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous) More...
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous) More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
int m_HomoDirec
 number of homogenous directions More...
 

Private Member Functions

void SetUpDomainInterfaces (void)
 
void FillDataFromInterfacePoint (InterfacePointShPtr &I, const Array< OneD, const Array< OneD, NekDouble > > &field, NekDouble &A, NekDouble &u, NekDouble &beta, NekDouble &A_0, NekDouble &alpha)
 
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_data * m_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 ( )
virtual

Destructor.

Destructor

Definition at line 75 of file PulseWaveSystem.cpp.

76{
77}

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

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

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

1523{
1524
1525 NekDouble rho = m_rho;
1526 Array<OneD, NekDouble> W(3);
1527 Array<OneD, NekDouble> P_Au(3);
1528 Array<OneD, NekDouble> W_Au(3);
1529 NekMatrix<NekDouble> invJ(6, 6);
1530 NekVector<NekDouble> f(6);
1531 NekVector<NekDouble> dx(6);
1532
1533 int proceed = 1;
1534 int iter = 0;
1535 int MAX_ITER = 100;
1536
1537 // Forward characteristic
1538 m_pressureArea->GetW1(W[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1539
1540 // Backward characteristics
1541 for (size_t i = 1; i < 3; ++i)
1542 {
1543 m_pressureArea->GetW2(W[i], uu[i], beta[i], Au[i], A_0[i], alpha[i]);
1544 }
1545
1546 // Tolerances for the algorithm
1547 NekDouble Tol = 1.0E-10;
1548
1549 // Newton Iteration
1550 while ((proceed) && (iter < MAX_ITER))
1551 {
1552 LibUtilities::Timer time_BifurcationRiemann;
1553 time_BifurcationRiemann.Start();
1554 iter += 1;
1555
1556 /*
1557 * We solve the six constraint equations via a multivariate Newton
1558 * iteration. Equations are:
1559 * 1. Forward characteristic: W1(A_L, U_L) = W1(Au_L, Uu_L)
1560 * 2. Backward characteristic 1: W2(A_R1, U_R1) = W2(Au_R1, Uu_R1)
1561 * 3. Backward characteristic 2: W2(A_R2, U_R2) = W2(Au_R2, Uu_R2)
1562 * 4. Conservation of mass: Au_L * Uu_L = Au_R1 * Uu_R1 +
1563 * Au_R2 * Uu_R2
1564 * 5. Continuity of total pressure 1: rho * Uu_L * Uu_L / 2 + p(Au_L)
1565 * = rho * Uu_R1 * Uu_R1 / 2 + p(Au_R1)
1566 * 6. Continuity of total pressure 2: rho * Uu_L * Uu_L / 2 + p(Au_L)
1567 * = rho * Uu_R2 * Uu_R2 / 2 + p(Au_R2)
1568 */
1569 for (size_t i = 0; i < 3; ++i)
1570 {
1571 m_pressureArea->GetPressure(P_Au[i], beta[i], Au[i], A_0[i], 0, 0,
1572 alpha[i]);
1573 }
1574
1575 m_pressureArea->GetW1(W_Au[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1576 m_pressureArea->GetW2(W_Au[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
1577 m_pressureArea->GetW2(W_Au[2], uu[2], beta[2], Au[2], A_0[2], alpha[2]);
1578
1579 // Constraint equations set to zero
1580 f[0] = W_Au[0] - W[0];
1581 f[1] = W_Au[1] - W[1];
1582 f[2] = W_Au[2] - W[2];
1583 f[3] = Au[0] * uu[0] - Au[1] * uu[1] - Au[2] * uu[2];
1584 f[4] = uu[0] * uu[0] + 2 * P_Au[0] / rho - uu[1] * uu[1] -
1585 2 * P_Au[1] / rho;
1586 f[5] = uu[0] * uu[0] + 2 * P_Au[0] / rho - uu[2] * uu[2] -
1587 2 * P_Au[2] / rho;
1588
1589 // Inverse Jacobian for x + dx = x - J^(-1) * f(x)
1590 m_pressureArea->GetJacobianInverse(invJ, Au, uu, beta, A_0, alpha,
1591 "Bifurcation");
1592
1593 Multiply(dx, invJ, f);
1594
1595 // Update the solution: x_new = x_old - dx
1596 for (int i = 0; i < 3; ++i)
1597 {
1598 uu[i] -= dx[i];
1599 Au[i] -= dx[i + 3];
1600 }
1601
1602 // Check if the error of the solution is smaller than Tol
1603 if (Dot(dx, dx) < Tol)
1604 {
1605 proceed = 0;
1606 }
1607
1608 // Check if solver converges
1609 if (iter >= MAX_ITER)
1610 {
1611 ASSERTL0(false, "Riemann solver for Bifurcation did not converge.");
1612 }
1613 time_BifurcationRiemann.Stop();
1614 time_BifurcationRiemann.AccumulateRegion(
1615 "PulseWaveSystem::Bifurcation Riemann", 2);
1616 }
1617}
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Nektar::PulseWaveSystem::m_pressureArea
PulseWavePressureAreaSharedPtr m_pressureArea
Definition: PulseWaveSystem.h:122
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
Nektar::Multiply
void Multiply(NekMatrix< ResultDataType, StandardMatrixTag > &result, const NekMatrix< LhsDataType, LhsMatrixType > &lhs, const ResultDataType &rhs)
Definition: MatrixOperations.cpp:41
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::Dot
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 1844 of file PulseWaveSystem.cpp.

1845{
1846 std::stringstream outname;
1847 outname << m_sessionName << "_" << n << ".chk";
1848
1849 WriteVessels(outname.str());
1850}
Nektar::PulseWaveSystem::WriteVessels
void WriteVessels(const std::string &outname)
Write input fields to the given filename.
Definition: PulseWaveSystem.cpp:1856
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:369

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

1376{
1377 LibUtilities::Timer time_EnforceInterfaceConditions;
1378 time_EnforceInterfaceConditions.Start();
1379 int dom, bcpos;
1380
1381 int totif =
1382 m_bifurcations.size() + m_mergingJcts.size() + m_vesselIntfcs.size();
1383
1384 Array<OneD, NekDouble> Aut, Au(3 * totif, 0.0);
1385 Array<OneD, NekDouble> uut, uu(3 * totif, 0.0);
1386 Array<OneD, NekDouble> betat, beta(3 * totif, 0.0);
1387 Array<OneD, NekDouble> A_0t, A_0(3 * totif, 0.0);
1388 Array<OneD, NekDouble> alphat, alpha(3 * totif, 0.0);
1389
1390 // Bifurcations Data:
1391 size_t cnt = 0;
1392 for (size_t n = 0; n < m_bifurcations.size(); ++n, ++cnt)
1393 {
1394 for (size_t i = 0; i < m_bifurcations[n].size(); ++i)
1395 {
1396 size_t l = m_bifurcations[n][i]->m_riemannOrd;
1397 FillDataFromInterfacePoint(
1398 m_bifurcations[n][i], fields, Au[l + 3 * cnt], uu[l + 3 * cnt],
1399 beta[l + 3 * cnt], A_0[l + 3 * cnt], alpha[l + 3 * cnt]);
1400 }
1401 }
1402
1403 // Enforce Merging vessles Data:
1404 for (size_t n = 0; n < m_mergingJcts.size(); ++n, ++cnt)
1405 {
1406 // Merged vessel
1407 for (size_t i = 0; i < m_mergingJcts.size(); ++i)
1408 {
1409 int l = m_mergingJcts[n][i]->m_riemannOrd;
1410 FillDataFromInterfacePoint(
1411 m_mergingJcts[n][i], fields, Au[l + 3 * cnt], uu[l + 3 * cnt],
1412 beta[l + 3 * cnt], A_0[l + 3 * cnt], alpha[l + 3 * cnt]);
1413 }
1414 }
1415
1416 // Enforce interface junction between two vessesls
1417 for (size_t n = 0; n < m_vesselIntfcs.size(); ++n, ++cnt)
1418 {
1419 for (size_t i = 0; i < m_vesselIntfcs[n].size(); ++i)
1420 {
1421 int l = m_vesselIntfcs[n][i]->m_riemannOrd;
1422 FillDataFromInterfacePoint(
1423 m_vesselIntfcs[n][i], fields, Au[l + 3 * cnt], uu[l + 3 * cnt],
1424 beta[l + 3 * cnt], A_0[l + 3 * cnt], alpha[l + 3 * cnt]);
1425 }
1426 }
1427
1428 // Gather data if running in parallel
1429 Gs::Gather(Au, Gs::gs_add, m_intComm);
1430 Gs::Gather(uu, Gs::gs_add, m_intComm);
1431 Gs::Gather(beta, Gs::gs_add, m_intComm);
1432 Gs::Gather(A_0, Gs::gs_add, m_intComm);
1433 Gs::Gather(alpha, Gs::gs_add, m_intComm);
1434
1435 // Enforce Bifurcations:
1436 cnt = 0;
1437 for (size_t n = 0; n < m_bifurcations.size(); ++n, ++cnt)
1438 {
1439 // Solve the Riemann problem for a bifurcation
1440 BifurcationRiemann(Aut = Au + 3 * cnt, uut = uu + 3 * cnt,
1441 betat = beta + 3 * cnt, A_0t = A_0 + 3 * cnt,
1442 alphat = alpha + 3 * cnt);
1443
1444 // Store the values into the right positions:
1445 for (size_t i = 0; i < m_bifurcations[n].size(); ++i)
1446 {
1447 dom = m_bifurcations[n][i]->m_domain;
1448 bcpos = m_bifurcations[n][i]->m_bcPosition;
1449 int l = m_bifurcations[n][i]->m_riemannOrd;
1450 m_vessels[dom * m_nVariables]
1451 ->UpdateBndCondExpansion(bcpos)
1452 ->UpdatePhys()[0] = Au[l + 3 * cnt];
1453 m_vessels[dom * m_nVariables + 1]
1454 ->UpdateBndCondExpansion(bcpos)
1455 ->UpdatePhys()[0] = uu[l + 3 * cnt];
1456 }
1457 }
1458
1459 // Enforce Merging vessles;
1460 for (size_t n = 0; n < m_mergingJcts.size(); ++n, ++cnt)
1461 {
1462 // Solve the Riemann problem for a merging vessel
1463 MergingRiemann(Aut = Au + 3 * cnt, uut = uu + 3 * cnt,
1464 betat = beta + 3 * cnt, A_0t = A_0 + 3 * cnt,
1465 alphat = alpha + 3 * cnt);
1466
1467 // Store the values into the right positions:
1468 for (size_t i = 0; i < m_mergingJcts[n].size(); ++i)
1469 {
1470 int dom = m_mergingJcts[n][i]->m_domain;
1471 int bcpos = m_mergingJcts[n][i]->m_bcPosition;
1472 int l = m_mergingJcts[n][i]->m_riemannOrd;
1473 m_vessels[dom * m_nVariables]
1474 ->UpdateBndCondExpansion(bcpos)
1475 ->UpdatePhys()[0] = Au[l + 3 * cnt];
1476 m_vessels[dom * m_nVariables + 1]
1477 ->UpdateBndCondExpansion(bcpos)
1478 ->UpdatePhys()[0] = uu[l + 3 * cnt];
1479 }
1480 }
1481
1482 // Enforce interface junction between two vessesls
1483 for (size_t n = 0; n < m_vesselIntfcs.size(); ++n, ++cnt)
1484 {
1485 InterfaceRiemann(Aut = Au + 3 * cnt, uut = uu + 3 * cnt,
1486 betat = beta + 3 * cnt, A_0t = A_0 + 3 * cnt,
1487 alphat = alpha + 3 * cnt);
1488
1489 // Store the values into the right positions:
1490 for (size_t i = 0; i < m_vesselIntfcs[n].size(); ++i)
1491 {
1492 int dom = m_vesselIntfcs[n][i]->m_domain;
1493 int bcpos = m_vesselIntfcs[n][i]->m_bcPosition;
1494 int l = m_vesselIntfcs[n][i]->m_riemannOrd;
1495 m_vessels[dom * m_nVariables]
1496 ->UpdateBndCondExpansion(bcpos)
1497 ->UpdatePhys()[0] = Au[l + 3 * cnt];
1498 m_vessels[dom * m_nVariables + 1]
1499 ->UpdateBndCondExpansion(bcpos)
1500 ->UpdatePhys()[0] = uu[l + 3 * cnt];
1501 }
1502 }
1503 time_EnforceInterfaceConditions.Stop();
1504 time_EnforceInterfaceConditions.AccumulateRegion(
1505 "PulseWaveSystem::EnforceInterfaceConditions", 1);
1506}
Nektar::PulseWaveSystem::MergingRiemann
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.
Definition: PulseWaveSystem.cpp:1628
Nektar::PulseWaveSystem::InterfaceRiemann
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.
Definition: PulseWaveSystem.cpp:1741
Nektar::PulseWaveSystem::m_intComm
Gs::gs_data * m_intComm
Communicator for interfaces.
Definition: PulseWaveSystem.h:195
Nektar::PulseWaveSystem::m_vesselIntfcs
std::vector< std::vector< InterfacePointShPtr > > m_vesselIntfcs
Definition: PulseWaveSystem.h:130
Nektar::PulseWaveSystem::m_vessels
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
Definition: PulseWaveSystem.h:94
Nektar::PulseWaveSystem::FillDataFromInterfacePoint
void FillDataFromInterfacePoint(InterfacePointShPtr &I, const Array< OneD, const Array< OneD, NekDouble > > &field, NekDouble &A, NekDouble &u, NekDouble &beta, NekDouble &A_0, NekDouble &alpha)
Definition: PulseWaveSystem.cpp:1342
Nektar::PulseWaveSystem::m_bifurcations
std::vector< std::vector< InterfacePointShPtr > > m_bifurcations
Definition: PulseWaveSystem.h:131
Nektar::PulseWaveSystem::m_mergingJcts
std::vector< std::vector< InterfacePointShPtr > > m_mergingJcts
Definition: PulseWaveSystem.h:132
Nektar::PulseWaveSystem::BifurcationRiemann
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.
Definition: PulseWaveSystem.cpp:1518
Gs::Gather
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:270
Gs::gs_add
@ gs_add
Definition: GsLib.hpp:62

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

1346{
1347 LibUtilities::Timer time_FillDataFromInterfacePoint;
1348 time_FillDataFromInterfacePoint.Start();
1349
1350 int omega = I->m_domain;
1351 int traceId = I->m_traceId;
1352 int eid = I->m_elmt;
1353 int vert = I->m_elmtVert;
1354 int vesselID = omega * m_nVariables;
1355 int phys_offset = m_vessels[vesselID]->GetPhys_Offset(eid);
1356 LocalRegions::ExpansionSharedPtr dumExp;
1357 Array<OneD, NekDouble> Tmp(1);
1358
1359 m_vessels[vesselID]->GetExp(eid)->GetTracePhysVals(
1360 vert, dumExp, fields[0] + m_fieldPhysOffset[omega] + phys_offset, Tmp);
1361 A = Tmp[0];
1362 m_vessels[vesselID]->GetExp(eid)->GetTracePhysVals(
1363 vert, dumExp, fields[1] + m_fieldPhysOffset[omega] + phys_offset, Tmp);
1364 u = Tmp[0];
1365
1366 beta = m_beta_trace[omega][traceId];
1367 A_0 = m_A_0_trace[omega][traceId];
1368 alpha = m_alpha_trace[omega][traceId];
1369 time_FillDataFromInterfacePoint.Stop();
1370 time_FillDataFromInterfacePoint.AccumulateRegion(
1371 "PulseWaveSystem::FillDataFromInterfacePoint", 2);
1372}
Nektar::PulseWaveSystem::m_fieldPhysOffset
Array< OneD, int > m_fieldPhysOffset
Definition: PulseWaveSystem.h:100
Nektar::PulseWaveSystem::m_beta_trace
Array< OneD, Array< OneD, NekDouble > > m_beta_trace
Definition: PulseWaveSystem.h:111
Nektar::PulseWaveSystem::m_alpha_trace
Array< OneD, Array< OneD, NekDouble > > m_alpha_trace
Definition: PulseWaveSystem.h:114
Nektar::PulseWaveSystem::m_A_0_trace
Array< OneD, Array< OneD, NekDouble > > m_A_0_trace
Definition: PulseWaveSystem.h:109
Nektar::LocalRegions::ExpansionSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68

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

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

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 83 of file PulseWaveSystem.h.

84 {
85 return m_nDomains;
86 }

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

1746{
1747
1748 NekDouble rho = m_rho;
1749 Array<OneD, NekDouble> W(2);
1750 Array<OneD, NekDouble> W_Au(2);
1751 Array<OneD, NekDouble> P_Au(2);
1752 NekMatrix<NekDouble> invJ(4, 4);
1753 NekVector<NekDouble> f(4);
1754 NekVector<NekDouble> dx(4);
1755
1756 int proceed = 1;
1757 int iter = 0;
1758 int MAX_ITER = 15;
1759 NekDouble Tol = 1.0E-10;
1760
1761 // Forward and backward characteristics
1762 m_pressureArea->GetW1(W[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1763 m_pressureArea->GetW2(W[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
1764
1765 while ((proceed) && (iter < MAX_ITER))
1766 {
1767 LibUtilities::Timer time_InterfaceRiemann;
1768 time_InterfaceRiemann.Start();
1769 iter += 1;
1770
1771 /*
1772 * We solve the four constraint equations via a multivariate Newton
1773 * iteration. Equations are:
1774 * 1. Forward characteristic: W1(A_L, U_L) = W1(Au_L, Uu_L)
1775 * 2. Backward characteristic: W2(A_R, U_R) = W2(Au_R, Uu_R)
1776 * 3. Conservation of mass: Au_L * Uu_L = Au_R * Uu_R
1777 * 4. Continuity of total pressure: rho * Uu_L * Uu_L / 2 + p(Au_L) =
1778 * rho * Uu_R * Uu_R / 2 + p(Au_R)
1779 */
1780 for (size_t i = 0; i < 2; ++i)
1781 {
1782 m_pressureArea->GetPressure(P_Au[i], beta[i], Au[i], A_0[i], 0, 0,
1783 alpha[i]);
1784 }
1785
1786 m_pressureArea->GetW1(W_Au[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1787 m_pressureArea->GetW2(W_Au[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
1788
1789 // Constraint equations set to zero
1790 f[0] = W_Au[0] - W[0];
1791 f[1] = W_Au[1] - W[1];
1792 f[2] = Au[0] * uu[0] - Au[1] * uu[1];
1793 f[3] = uu[0] * uu[0] + 2 * P_Au[0] / rho - uu[1] * uu[1] -
1794 2 * P_Au[1] / rho;
1795
1796 // Inverse Jacobian for x + dx = x - J^(-1) * f(x)
1797 m_pressureArea->GetJacobianInverse(invJ, Au, uu, beta, A_0, alpha,
1798 "Interface");
1799
1800 Multiply(dx, invJ, f);
1801
1802 // Update solution: x_new = x_old - dx
1803 for (size_t i = 0; i < 2; ++i)
1804 {
1805 uu[i] -= dx[i];
1806 Au[i] -= dx[i + 2];
1807 }
1808
1809 // Check if the error of the solution is smaller than Tol.
1810 if (Dot(dx, dx) < Tol)
1811 {
1812 proceed = 0;
1813 }
1814 time_InterfaceRiemann.Stop();
1815 time_InterfaceRiemann.AccumulateRegion(
1816 "PulseWaveSystem::InterfaceRiemann", 2);
1817 }
1818
1819 if (iter >= MAX_ITER)
1820 {
1821 ASSERTL0(false, "Riemann solver for Interface did not converge");
1822 }
1823}

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

1633{
1634
1635 NekDouble rho = m_rho;
1636 Array<OneD, NekDouble> W(3);
1637 Array<OneD, NekDouble> W_Au(3);
1638 Array<OneD, NekDouble> P_Au(3);
1639 NekMatrix<NekDouble> invJ(6, 6);
1640 NekVector<NekDouble> f(6);
1641 NekVector<NekDouble> dx(6);
1642
1643 int proceed = 1;
1644 int iter = 0;
1645 int MAX_ITER = 15;
1646
1647 // Backward characteristic
1648 m_pressureArea->GetW2(W[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1649
1650 // Forward characteristics
1651 for (size_t i = 1; i < 3; ++i)
1652 {
1653 m_pressureArea->GetW1(W[i], uu[i], beta[i], Au[i], A_0[i], alpha[i]);
1654 }
1655
1656 // Tolerances for the algorithm
1657 NekDouble Tol = 1.0E-10;
1658
1659 // Newton Iteration
1660 while ((proceed) && (iter < MAX_ITER))
1661 {
1662 LibUtilities::Timer time_MergingRiemann;
1663 time_MergingRiemann.Start();
1664 iter += 1;
1665
1666 /*
1667 * We solve the six constraint equations via a multivariate Newton
1668 * iteration. Equations are:
1669 * 1. Backward characteristic: W2(A_R, U_R) = W1(Au_R, Uu_R)
1670 * 2. Forward characteristic 1: W1(A_L1, U_L1) = W1(Au_L1,
1671 * Uu_R1)
1672 * 3. Forward characteristic 2: W1(A_L2, U_L2) = W1(Au_L2,
1673 * Uu_L2)
1674 * 4. Conservation of mass: Au_R * Uu_R = Au_L1 * Uu_L1 +
1675 * Au_L2 * Uu_L2
1676 * 5. Continuity of total pressure 1: rho * Uu_R * Uu_R / 2 + p(Au_R)
1677 * = rho * Uu_L1 * Uu_L1 / 2 + p(Au_L1)
1678 * 6. Continuity of total pressure 2: rho * Uu_R * Uu_R / 2 + p(Au_R)
1679 * = rho * Uu_L2 * Uu_L2 / 2 + p(Au_L2)
1680 */
1681 for (size_t i = 0; i < 3; ++i)
1682 {
1683 m_pressureArea->GetPressure(P_Au[i], beta[i], Au[i], A_0[i], 0, 0,
1684 alpha[i]);
1685 }
1686
1687 m_pressureArea->GetW2(W_Au[0], uu[0], beta[0], Au[0], A_0[0], alpha[0]);
1688 m_pressureArea->GetW1(W_Au[1], uu[1], beta[1], Au[1], A_0[1], alpha[1]);
1689 m_pressureArea->GetW1(W_Au[2], uu[2], beta[2], Au[2], A_0[2], alpha[2]);
1690
1691 // Constraint equations set to zero
1692 f[0] = W_Au[0] - W[0];
1693 f[1] = W_Au[1] - W[1];
1694 f[2] = W_Au[2] - W[2];
1695 f[3] = Au[0] * uu[0] - Au[1] * uu[1] - Au[2] * uu[2];
1696 f[4] = uu[0] * uu[0] + 2 * P_Au[0] / rho - uu[1] * uu[1] -
1697 2 * P_Au[1] / rho;
1698 f[5] = uu[0] * uu[0] + 2 * P_Au[0] / rho - uu[2] * uu[2] -
1699 2 * P_Au[2] / rho;
1700
1701 // Inverse Jacobian for x + dx = x - J^(-1) * f(x)
1702 m_pressureArea->GetJacobianInverse(invJ, Au, uu, beta, A_0, alpha,
1703 "Merge");
1704
1705 Multiply(dx, invJ, f);
1706
1707 // Update the solution: x_new = x_old - dx
1708 for (size_t i = 0; i < 3; ++i)
1709 {
1710 uu[i] -= dx[i];
1711 Au[i] -= dx[i + 3];
1712 }
1713
1714 // Check if the error of the solution is smaller than Tol
1715 if (Dot(dx, dx) < Tol)
1716 {
1717 proceed = 0;
1718 }
1719
1720 // Check if solver converges
1721 if (iter >= MAX_ITER)
1722 {
1723 ASSERTL0(false, "Riemann solver for Merging Flow did not converge");
1724 }
1725 time_MergingRiemann.Stop();
1726 time_MergingRiemann.AccumulateRegion("PulseWaveSystem::MergingRiemann",
1727 2);
1728 }
1729}

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

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

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

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

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 88 of file PulseWaveSystem.h.

89 {
90 return m_vessels;
91 }

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

1220{
1221 boost::ignore_unused(dumpInitialConditions);
1222
1223 if (m_session->GetComm()->GetRank() == 0)
1224 {
1225 cout << "Initial Conditions: " << endl;
1226 }
1227
1228 /* Loop over all subdomains to initialize all with the Initial
1229 * Conditions read from the inputfile*/
1230 int omega = 0;
1231 for (auto &d : m_domOrder)
1232 {
1233 for (size_t i = 0; i < 2; ++i)
1234 {
1235 m_fields[i] = m_vessels[m_nVariables * omega + i];
1236 }
1237
1238 if (m_session->GetComm()->GetRank() == 0)
1239 {
1240 cout << "Subdomain: " << omega << endl;
1241 }
1242
1243 SetInitialConditions(0.0, false, d);
1244 omega++;
1245 }
1246
1247 // Reset 2 variables to first vessels
1248 for (size_t i = 0; i < 2; ++i)
1249 {
1250 m_fields[i] = m_vessels[i];
1251 }
1252}
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:363
Nektar::SolverUtils::EquationSystem::SetInitialConditions
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
Definition: EquationSystem.h:675

References Nektar::UnitTests::d(), m_domOrder, Nektar::SolverUtils::EquationSystem::m_fields, m_nVariables, Nektar::SolverUtils::EquationSystem::m_session, 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 1271 of file PulseWaveSystem.cpp.

1272{
1273 // NekDouble IntegrationTime = 0.0;
1274 size_t i;
1275 int n;
1276 int nchk = 1;
1277
1278 Array<OneD, Array<OneD, NekDouble>> fields(m_nVariables);
1279
1280 for (i = 0; i < m_nVariables; ++i)
1281 {
1282 fields[i] = m_vessels[i]->UpdatePhys();
1283 m_fields[i]->SetPhysState(false);
1284 }
1285
1286 m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
1287
1288 // Time loop
1289 for (n = 0; n < m_steps; ++n)
1290 {
1291 LibUtilities::Timer time_v_IntStep;
1292 time_v_IntStep.Start();
1293 fields = m_intScheme->TimeIntegrate(n, m_timestep);
1294 m_time += m_timestep;
1295 time_v_IntStep.Stop();
1296 time_v_IntStep.AccumulateRegion("PulseWaveSystem::TimeIntegrationStep",
1297 0);
1298 // IntegrationTime += timer.TimePerTest(1);
1299
1300 // Write out status information.
1301 if (m_infosteps && !((n + 1) % m_infosteps) &&
1302 m_session->GetComm()->GetRank() == 0)
1303 {
1304 cout << "Steps: " << n + 1 << "\t Time: " << m_time
1305 << "\t Time-step: " << m_timestep << "\t" << endl;
1306 }
1307
1308 // Transform data if needed
1309 if (m_checksteps && !((n + 1) % m_checksteps))
1310 {
1311 for (i = 0; i < m_nVariables; ++i)
1312 {
1313 size_t cnt = 0;
1314 for (size_t omega = 0; omega < m_nDomains; omega++)
1315 {
1316 m_vessels[omega * m_nVariables + i]->FwdTrans(
1317 fields[i] + cnt,
1318 m_vessels[omega * m_nVariables + i]->UpdateCoeffs());
1319 cnt += m_vessels[omega * m_nVariables + i]->GetTotPoints();
1320 }
1321 }
1322 CheckPoint_Output(nchk++);
1323 }
1324
1325 } // end of timeintegration
1326
1327 // Copy Array To Vessel Phys Fields
1328 for (size_t i = 0; i < m_nVariables; ++i)
1329 {
1330 Vmath::Vcopy(fields[i].size(), fields[i], 1, m_vessels[i]->UpdatePhys(),
1331 1);
1332 }
1333
1334 /** if(m_session->GetComm()->GetRank() == 0)
1335 {
1336 cout << "Time-integration timing: " << IntegrationTime << " s" <<
1337 endl
1338 << endl;
1339 } **/
1340}
Nektar::PulseWaveSystem::CheckPoint_Output
void CheckPoint_Output(const int n)
Definition: PulseWaveSystem.cpp:1844
Nektar::SolverUtils::EquationSystem::m_timestep
NekDouble m_timestep
Time step size.
Definition: EquationSystem.h:377
Nektar::SolverUtils::EquationSystem::m_infosteps
int m_infosteps
Number of time steps between outputting status information.
Definition: EquationSystem.h:394
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:371
Nektar::SolverUtils::EquationSystem::m_steps
int m_steps
Number of steps to take.
Definition: EquationSystem.h:390
Nektar::SolverUtils::EquationSystem::m_checksteps
int m_checksteps
Number of steps between checkpoints.
Definition: EquationSystem.h:392
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:79
Nektar::SolverUtils::UnsteadySystem::m_intScheme
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
Definition: UnsteadySystem.h:77
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1191

References Nektar::LibUtilities::Timer::AccumulateRegion(), CheckPoint_Output(), Nektar::SolverUtils::EquationSystem::m_checksteps, Nektar::SolverUtils::EquationSystem::m_fields, 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.

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

Definition at line 88 of file PulseWaveSystem.cpp.

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

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_domOrder, m_fieldPhysOffset, Nektar::SolverUtils::EquationSystem::m_fields, m_gamma, Nektar::SolverUtils::EquationSystem::m_graph, m_nDomains, m_nVariables, m_pext, m_pressure, m_pressureArea, Nektar::SolverUtils::EquationSystem::m_projectionType, m_PWV, m_rho, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_time, m_trace_fwd_normal, m_upwindTypePulse, m_vessels, m_W1, m_W2, SetUpDomainInterfaceBCs(), SetUpDomainInterfaces(), Nektar::SIZE_UpwindTypePulse, Nektar::UpwindTypeMapPulse, and Nektar::SolverUtils::UnsteadySystem::v_InitObject().

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

◆ v_L2Error()

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

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1950 of file PulseWaveSystem.cpp.

1953{
1954 boost::ignore_unused(exact);
1955
1956 NekDouble L2error = 0.0;
1957 NekDouble L2error_dom;
1958 NekDouble Vol = 0.0;
1959
1960 if (m_NumQuadPointsError == 0)
1961 {
1962 for (size_t omega = 0; omega < m_nDomains; omega++)
1963 {
1964 size_t vesselid = field + omega * m_nVariables;
1965
1966 // since exactsoln is passed for just the first field size we need
1967 // to reset it each domain
1968 Array<OneD, NekDouble> exactsoln(m_vessels[vesselid]->GetNpoints());
1969 m_fields[field] = m_vessels[omega * m_nVariables];
1970 EvaluateExactSolution(field, exactsoln, m_time);
1971
1972 if (m_vessels[vesselid]->GetPhysState() == false)
1973 {
1974 m_vessels[vesselid]->BwdTrans(
1975 m_vessels[vesselid]->GetCoeffs(),
1976 m_vessels[vesselid]->UpdatePhys());
1977 }
1978
1979 if (exactsoln.size())
1980 {
1981 L2error_dom = m_vessels[vesselid]->L2(
1982 m_vessels[vesselid]->GetPhys(), exactsoln);
1983 }
1984 else if (m_session->DefinesFunction("ExactSolution"))
1985 {
1986 Array<OneD, NekDouble> exactsoln(
1987 m_vessels[vesselid]->GetNpoints());
1988
1989 LibUtilities::EquationSharedPtr vEqu =
1990 m_session->GetFunction("ExactSolution", field, omega);
1991 GetFunction("ExactSolution")
1992 ->Evaluate(m_session->GetVariable(field), exactsoln,
1993 m_time);
1994
1995 L2error_dom = m_vessels[vesselid]->L2(
1996 m_vessels[vesselid]->GetPhys(), exactsoln);
1997 }
1998 else
1999 {
2000 L2error_dom =
2001 m_vessels[vesselid]->L2(m_vessels[vesselid]->GetPhys());
2002 }
2003
2004 if (m_vessels[vesselid]->GetComm()->GetRank())
2005 {
2006 // ensure domains are only summed on local root of
2007 // domain
2008 L2error_dom = 0.0;
2009 }
2010
2011 L2error += L2error_dom * L2error_dom;
2012
2013 if (Normalised == true)
2014 {
2015 Array<OneD, NekDouble> one(m_vessels[vesselid]->GetNpoints(),
2016 1.0);
2017
2018 Vol += m_vessels[vesselid]->Integral(one);
2019 }
2020 }
2021 }
2022 else
2023 {
2024 ASSERTL0(false, "Not set up");
2025 }
2026
2027 m_comm->AllReduce(L2error, LibUtilities::ReduceSum);
2028
2029 if (Normalised == true)
2030 {
2031 m_comm->AllReduce(Vol, LibUtilities::ReduceSum);
2032
2033 L2error = sqrt(L2error / Vol);
2034 }
2035 else
2036 {
2037 L2error = sqrt(L2error);
2038 }
2039
2040 return L2error;
2041}
Nektar::SolverUtils::EquationSystem::EvaluateExactSolution
SOLVER_UTILS_EXPORT void EvaluateExactSolution(int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Evaluates an exact solution.
Definition: EquationSystem.h:683
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:775
Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Definition: EquationSystem.h:443
Nektar::LibUtilities::EquationSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:129
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References ASSERTL0, Nektar::SolverUtils::EquationSystem::EvaluateExactSolution(), 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 2049 of file PulseWaveSystem.cpp.

2051{
2052 boost::ignore_unused(exact);
2053
2054 NekDouble LinferrorDom, Linferror = -1.0;
2055
2056 for (size_t omega = 0; omega < m_nDomains; ++omega)
2057 {
2058 size_t vesselid = field + omega * m_nVariables;
2059
2060 // since exactsoln is passed for just the first field size we need
2061 // to reset it each domain
2062 Array<OneD, NekDouble> exactsoln(m_vessels[vesselid]->GetNpoints());
2063 m_fields[field] = m_vessels[omega * m_nVariables];
2064 EvaluateExactSolution(field, exactsoln, m_time);
2065
2066 if (m_NumQuadPointsError == 0)
2067 {
2068 if (m_vessels[vesselid]->GetPhysState() == false)
2069 {
2070 m_vessels[vesselid]->BwdTrans(
2071 m_vessels[vesselid]->GetCoeffs(),
2072 m_vessels[vesselid]->UpdatePhys());
2073 }
2074
2075 if (exactsoln.size())
2076 {
2077 LinferrorDom = m_vessels[vesselid]->Linf(
2078 m_vessels[vesselid]->GetPhys(), exactsoln);
2079 }
2080 else if (m_session->DefinesFunction("ExactSolution"))
2081 {
2082 Array<OneD, NekDouble> exactsoln(
2083 m_vessels[vesselid]->GetNpoints());
2084
2085 GetFunction("ExactSolution")
2086 ->Evaluate(m_session->GetVariable(field), exactsoln,
2087 m_time);
2088
2089 LinferrorDom = m_vessels[vesselid]->Linf(
2090 m_vessels[vesselid]->GetPhys(), exactsoln);
2091 }
2092 else
2093 {
2094 LinferrorDom = 0.0;
2095 }
2096
2097 Linferror = (Linferror > LinferrorDom) ? Linferror : LinferrorDom;
2098 }
2099 else
2100 {
2101 ASSERTL0(false, "ErrorExtraPoints not allowed for this solver");
2102 }
2103 }
2104 m_comm->AllReduce(Linferror, LibUtilities::ReduceMax);
2105 return Linferror;
2106}

References ASSERTL0, Nektar::SolverUtils::EquationSystem::EvaluateExactSolution(), 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 1829 of file PulseWaveSystem.cpp.

1830{
1831 /**
1832 * Write the field data to file. The file is named according to the session
1833 * name with the extension .fld appended.
1834 */
1835 std::string outname = m_sessionName + ".fld";
1836
1837 WriteVessels(outname);
1838}

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

1857{
1858 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1859 std::vector<std::string> variables = m_session->GetVariables();
1860
1861 for (size_t n = 0; n < m_nDomains; ++n)
1862 {
1863 m_vessels[n * m_nVariables]->GetFieldDefinitions(FieldDef);
1864 }
1865
1866 std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
1867
1868 size_t nFieldDefPerDomain = FieldDef.size() / m_nDomains;
1869 size_t cnt;
1870
1871 // Copy Data into FieldData and set variable
1872 for (size_t n = 0; n < m_nDomains; ++n)
1873 {
1874 // Outputs area and velocity
1875 for (size_t j = 0; j < m_nVariables; ++j)
1876 {
1877 for (size_t i = 0; i < nFieldDefPerDomain; ++i)
1878 {
1879 cnt = n * nFieldDefPerDomain + i;
1880
1881 FieldDef[cnt]->m_fields.push_back(variables[j]);
1882
1883 m_vessels[n * m_nVariables]->AppendFieldData(
1884 FieldDef[cnt], FieldData[cnt],
1885 m_vessels[n * m_nVariables + j]->UpdateCoeffs());
1886 }
1887 }
1888
1889 // Outputs pressure
1890 Array<OneD, NekDouble> PFwd(m_vessels[n * m_nVariables]->GetNcoeffs());
1891
1892 m_vessels[n * m_nVariables]->FwdTransLocalElmt(m_pressure[n], PFwd);
1893
1894 for (size_t i = 0; i < nFieldDefPerDomain; ++i)
1895 {
1896 cnt = n * nFieldDefPerDomain + i;
1897
1898 FieldDef[cnt]->m_fields.push_back("P");
1899
1900 m_vessels[n * m_nVariables]->AppendFieldData(FieldDef[cnt],
1901 FieldData[cnt], PFwd);
1902 }
1903
1904 if (extraFields)
1905 {
1906 Array<OneD, NekDouble> PWVFwd(
1907 m_vessels[n * m_nVariables]->GetNcoeffs());
1908 Array<OneD, NekDouble> W1Fwd(
1909 m_vessels[n * m_nVariables]->GetNcoeffs());
1910 Array<OneD, NekDouble> W2Fwd(
1911 m_vessels[n * m_nVariables]->GetNcoeffs());
1912
1913 m_vessels[n * m_nVariables]->FwdTransLocalElmt(m_PWV[n], PWVFwd);
1914 m_vessels[n * m_nVariables]->FwdTransLocalElmt(m_W1[n], W1Fwd);
1915 m_vessels[n * m_nVariables]->FwdTransLocalElmt(m_W2[n], W2Fwd);
1916
1917 for (size_t i = 0; i < nFieldDefPerDomain; ++i)
1918 {
1919 cnt = n * nFieldDefPerDomain + i;
1920
1921 FieldDef[cnt]->m_fields.push_back("c");
1922 FieldDef[cnt]->m_fields.push_back("W1");
1923 FieldDef[cnt]->m_fields.push_back("W2");
1924
1925 m_vessels[n * m_nVariables]->AppendFieldData(
1926 FieldDef[cnt], FieldData[cnt], PWVFwd);
1927 m_vessels[n * m_nVariables]->AppendFieldData(
1928 FieldDef[cnt], FieldData[cnt], W1Fwd);
1929 m_vessels[n * m_nVariables]->AppendFieldData(
1930 FieldDef[cnt], FieldData[cnt], W2Fwd);
1931 }
1932 }
1933 }
1934
1935 // Update time in field info if required
1936 if (m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
1937 {
1938 m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
1939 }
1940
1941 m_fld->Write(outname, FieldDef, FieldData, m_fieldMetaDataMap);
1942}
Nektar::SolverUtils::EquationSystem::m_fld
LibUtilities::FieldIOSharedPtr m_fld
Field input/output.
Definition: EquationSystem.h:361
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:707
Nektar::SolverUtils::EquationSystem::m_fieldMetaDataMap
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
Definition: EquationSystem.h:429

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 125 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

Definition at line 108 of file PulseWaveSystem.h.

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

◆ m_A_0_trace

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

Definition at line 109 of file PulseWaveSystem.h.

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

◆ m_alpha

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

Definition at line 113 of file PulseWaveSystem.h.

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

◆ m_alpha_trace

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

Definition at line 114 of file PulseWaveSystem.h.

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

◆ m_beta

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

Definition at line 110 of file PulseWaveSystem.h.

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

◆ m_beta_trace

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

Definition at line 111 of file PulseWaveSystem.h.

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

◆ m_bifurcations

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

Definition at line 131 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_C

NekDouble Nektar::PulseWaveSystem::m_C
protected

Definition at line 104 of file PulseWaveSystem.h.

◆ m_currentDomain

size_t Nektar::PulseWaveSystem::m_currentDomain
protected

Definition at line 96 of file PulseWaveSystem.h.

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

◆ m_domain

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

Definition at line 118 of file PulseWaveSystem.h.

Referenced by GetCommArray(), SetUpDomainInterfaceBCs(), SetUpDomainInterfaces(), and v_InitObject().

◆ m_domOrder

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

Definition at line 119 of file PulseWaveSystem.h.

Referenced by GetCommArray(), SetUpDomainInterfaceBCs(), v_DoInitialise(), and v_InitObject().

◆ m_fieldPhysOffset

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

Definition at line 100 of file PulseWaveSystem.h.

Referenced by FillDataFromInterfacePoint(), and v_InitObject().

◆ m_gamma

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

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

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_mergingJcts

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

Definition at line 132 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_nDomains

size_t Nektar::PulseWaveSystem::m_nDomains
protected

Definition at line 95 of file PulseWaveSystem.h.

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

◆ m_nVariables

size_t Nektar::PulseWaveSystem::m_nVariables
protected

Definition at line 97 of file PulseWaveSystem.h.

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

◆ m_pext

NekDouble Nektar::PulseWaveSystem::m_pext
protected

Definition at line 102 of file PulseWaveSystem.h.

Referenced by v_InitObject().

◆ m_pout

NekDouble Nektar::PulseWaveSystem::m_pout
protected

Definition at line 106 of file PulseWaveSystem.h.

◆ m_pressure

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

Definition at line 121 of file PulseWaveSystem.h.

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

◆ m_pressureArea

PulseWavePressureAreaSharedPtr Nektar::PulseWaveSystem::m_pressureArea
protected

Definition at line 122 of file PulseWaveSystem.h.

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

◆ m_PWV

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

Definition at line 126 of file PulseWaveSystem.h.

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

◆ m_rho

NekDouble Nektar::PulseWaveSystem::m_rho
protected

Definition at line 101 of file PulseWaveSystem.h.

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

◆ m_RT

NekDouble Nektar::PulseWaveSystem::m_RT
protected

Definition at line 105 of file PulseWaveSystem.h.

◆ m_trace_fwd_normal

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

Definition at line 115 of file PulseWaveSystem.h.

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

◆ m_upwindTypePulse

UpwindTypePulse Nektar::PulseWaveSystem::m_upwindTypePulse
protected

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

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_vessels

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

Definition at line 94 of file PulseWaveSystem.h.

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

◆ m_W1

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

Definition at line 127 of file PulseWaveSystem.h.

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

◆ m_W2

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

Definition at line 128 of file PulseWaveSystem.h.

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

Generated on Wed Jan 31 2024 08:26:53 for Nektar++ by doxygen 1.9.4