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)
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void SetUpTraceNormals (void)
 
SOLVER_UTILS_EXPORT void InitObject (bool DeclareField=true)
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise ()
 Perform any initialisation necessary before solving the problem. More...
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem. More...
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space. More...
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space. More...
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve. More...
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name. More...
 
template<class T >
std::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name. More...
 
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession ()
 Get Session name. More...
 
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure ()
 Get pressure field if available. More...
 
SOLVER_UTILS_EXPORT void ExtraFldOutput (std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
 
SOLVER_UTILS_EXPORT void PrintSummary (std::ostream &out)
 Print a summary of parameters and solver characteristics. More...
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda. More...
 
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction (std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
 Get a SessionFunction by name. More...
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields. More...
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the L2 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > ErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
 Write checkpoint file of custom data fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
 Write input fields to the given filename. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
 Input field data from the given file to multiple domains. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, std::vector< std::string > &fieldStr, Array< OneD, Array< OneD, NekDouble >> &coeffs)
 Output a field. Input field data into array from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, MultiRegions::ExpListSharedPtr &pField, std::string &pFieldName)
 Output a field. Input field data into ExpList from the given file. More...
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary. More...
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMap & UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetFinalTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, 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 int GetPararealIterationNumber ()
 
SOLVER_UTILS_EXPORT void SetPararealIterationNumber (int num)
 
SOLVER_UTILS_EXPORT bool GetUseInitialCondition ()
 
SOLVER_UTILS_EXPORT void SetUseInitialCondition (bool 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...
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp ()
 Virtual function to 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 () 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...
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. 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)
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time, int &nchk)
 
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff (const Array< OneD, Array< OneD, NekDouble >> vel, StdRegions::VarCoeffMap &varCoeffMap)
 Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h time velocity. More...
 
virtual SOLVER_UTILS_EXPORT bool v_UpdateTimeStepCheck ()
 
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_TransCoeffToPhys ()
 Virtual function for transformation to physical space. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff ()
 Virtual function for transformation to coefficient space. More...
 
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution (unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure (void)
 
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
 

Protected Attributes

Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
 
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
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
int m_nanSteps
 
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
NekDouble m_epsilon
 
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used. More...
 
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used. More...
 
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used. More...
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
NekDouble m_steadyStateRes = 1.0
 
NekDouble m_steadyStateRes0 = 1.0
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check. More...
 
std::ofstream m_errFile
 
std::vector< int > m_intVariables
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
NekDouble m_filterTimeWarning
 Number of time steps between outputting status information. More...
 
NekDouble m_TimeIntegLambda = 0.0
 coefff of spacial derivatives(rhs or m_F in GLM) in calculating the residual of the whole equation(used in unsteady time integrations) More...
 
bool m_flagImplicitItsStatistics
 
bool m_flagImplicitSolver = false
 
Array< OneD, NekDouble > m_magnitdEstimat
 estimate the magnitude of each conserved varibles More...
 
Array< OneD, NekDouble > m_locTimeStep
 local time step(notice only for jfnk other see m_cflSafetyFactor) More...
 
NekDouble m_inArrayNorm = -1.0
 
int m_TotLinItePerStep = 0
 
int m_StagesPerStep = 1
 
bool m_flagUpdatePreconMat
 
int m_maxLinItePerNewton
 
int m_TotNewtonIts = 0
 
int m_TotLinIts = 0
 
int m_TotImpStages = 0
 
bool m_CalcPhysicalAV = true
 flag to update artificial viscosity More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
bool m_verbose
 
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_timestepMax = -1.0
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
NekDouble m_lastCheckTime
 
NekDouble m_TimeIncrementFactor
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_infosteps
 Number of time steps between outputting status information. More...
 
int m_pararealIter
 Number of parareal 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_useInitialCondition
 Flag to determine if IC are used. 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

- Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 
NekDouble m_cflNonAcoustic
 
NekDouble m_CFLGrowth
 CFL growth rate. More...
 
NekDouble m_CFLEnd
 maximun cfl in cfl growth More...
 
- 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 []
 

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:68

Member Function Documentation

◆ BifurcationRiemann()

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

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

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

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

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

1345 {
1346  LibUtilities::Timer time_FillDataFromInterfacePoint;
1347  time_FillDataFromInterfacePoint.Start();
1348 
1349  int omega = I->m_domain;
1350  int traceId = I->m_traceId;
1351  int eid = I->m_elmt;
1352  int vert = I->m_elmtVert;
1353  int vesselID = omega * m_nVariables;
1354  int phys_offset = m_vessels[vesselID]->GetPhys_Offset(eid);
1355  LocalRegions::ExpansionSharedPtr dumExp;
1356  Array<OneD, NekDouble> Tmp(1);
1357 
1358  m_vessels[vesselID]->GetExp(eid)->GetTracePhysVals(
1359  vert, dumExp, fields[0] + m_fieldPhysOffset[omega] + phys_offset, Tmp);
1360  A = Tmp[0];
1361  m_vessels[vesselID]->GetExp(eid)->GetTracePhysVals(
1362  vert, dumExp, fields[1] + m_fieldPhysOffset[omega] + phys_offset, Tmp);
1363  u = Tmp[0];
1364 
1365  beta = m_beta_trace[omega][traceId];
1366  A_0 = m_A_0_trace[omega][traceId];
1367  alpha = m_alpha_trace[omega][traceId];
1368  time_FillDataFromInterfacePoint.Stop();
1369  time_FillDataFromInterfacePoint.AccumulateRegion(
1370  "PulseWaveSystem::FillDataFromInterfacePoint", 2);
1371 }
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:365
Nektar::LibUtilities::CommSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:54
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:895
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:1023
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:945

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ASSERTL1, 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  }

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

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

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

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

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:368
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:53

References Nektar::SpatialDomains::BoundaryConditions::AddBoundaryConditions(), Nektar::SpatialDomains::BoundaryConditions::AddBoundaryRegions(), Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, 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::gs_max
@ gs_max
Definition: GsLib.hpp:65
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
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
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:492
Gs::gs_data
Definition: GsLib.hpp:151

References ASSERTL0, ASSERTL1, 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  }

◆ v_DoInitialise()

void Nektar::PulseWaveSystem::v_DoInitialise ( void  )
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 
1222  if (m_session->GetComm()->GetRank() == 0)
1223  {
1224  cout << "Initial Conditions: " << endl;
1225  }
1226 
1227  /* Loop over all subdomains to initialize all with the Initial
1228  * Conditions read from the inputfile*/
1229  int omega = 0;
1230  for (auto &d : m_domOrder)
1231  {
1232  for (size_t i = 0; i < 2; ++i)
1233  {
1234  m_fields[i] = m_vessels[m_nVariables * omega + i];
1235  }
1236 
1237  if (m_session->GetComm()->GetRank() == 0)
1238  {
1239  cout << "Subdomain: " << omega << endl;
1240  }
1241 
1242  SetInitialConditions(0.0, 0, d);
1243  omega++;
1244  }
1245 
1246  // Reset 2 variables to first vessels
1247  for (size_t i = 0; i < 2; ++i)
1248  {
1249  m_fields[i] = m_vessels[i];
1250  }
1251 }
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:375
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:684

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

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

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::PulseWaveSystemOutput, and Nektar::PulseWavePropagation.

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:379
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:743
Nektar::SolverUtils::EquationSystem::m_projectionType
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:436
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:717
Nektar::SolverUtils::UnsteadySystem::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:79
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::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 1949 of file PulseWaveSystem.cpp.

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

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

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

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

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

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

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(), 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(), 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 Sat Mar 25 2023 12:35:56 for Nektar++ by doxygen 1.9.1