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, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT Array< OneD, const Array< OneD, NekDouble > > GetTraceNormals ()
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp ()
 Virtual function to identify if operator is negated in DoSolve. More...
 

Protected Member Functions

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

Protected Attributes

Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
 
int m_nDomains
 
int m_currentDomain
 
int m_nVariables
 
UpwindTypePulse m_upwindTypePulse
 
Array< OneD, int > m_fieldPhysOffset
 
NekDouble m_rho
 
NekDouble m_pext
 
NekDouble m_C
 
NekDouble m_RT
 
NekDouble m_pout
 
Array< OneD, Array< OneD, NekDouble > > m_A_0
 
Array< OneD, Array< OneD, NekDouble > > m_A_0_trace
 
Array< OneD, Array< OneD, NekDouble > > m_beta
 
Array< OneD, Array< OneD, NekDouble > > m_beta_trace
 
Array< OneD, Array< OneD, NekDouble > > m_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_infosteps
 Number of time steps between outputting status information. More...
 
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
int m_nanSteps
 
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
NekDouble m_epsilon
 
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used. More...
 
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used. More...
 
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used. More...
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
NekDouble m_steadyStateRes = 1.0
 
NekDouble m_steadyStateRes0 = 1.0
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check. More...
 
std::ofstream m_errFile
 
std::vector< int > m_intVariables
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
NekDouble m_filterTimeWarning
 Number of time steps between outputting status information. More...
 
NekDouble m_TimeIntegLambda = 0.0
 coefff of spacial derivatives(rhs or m_F in GLM) in calculating the residual of the whole equation(used in unsteady time integrations) More...
 
bool m_flagImplicitItsStatistics
 
bool m_flagImplicitSolver = false
 
Array< OneD, 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
 
bool m_root
 
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_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
Array< OneD, NekDouble > m_movingFrameVelsxyz
 Moving frame of reference velocities. More...
 
Array< OneD, NekDouble > m_movingFrameTheta
 Moving frame of reference angles with respect to the. More...
 
boost::numeric::ublas::matrix< NekDouble > m_movingFrameProjMat
 Projection matrix for transformation between inertial and moving. More...
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error. More...
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous) More...
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous) More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
int m_HomoDirec
 number of homogenous directions More...
 

Private Member Functions

void SetUpDomainInterfaces (void)
 
void FillDataFromInterfacePoint (InterfacePointShPtr &I, const Array< OneD, const Array< OneD, NekDouble >> &field, NekDouble &A, NekDouble &u, NekDouble &beta, NekDouble &A_0, NekDouble &alpha)
 
void GetCommArray (std::map< int, LibUtilities::CommSharedPtr > &retval)
 Set and retrn a series of communicators for each partition. More...
 
void SetUpDomainInterfaceBCs (SpatialDomains::BoundaryConditions &AllBcs, std::map< int, LibUtilities::CommSharedPtr > &domComm)
 

Private Attributes

Gs::gs_data * m_intComm
 Communicator for interfaces. More...
 

Additional Inherited Members

- 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 (int 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 (int 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:341

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  int cnt = 0;
1391  for (int n = 0; n < m_bifurcations.size(); ++n, ++cnt)
1392  {
1393  for (int i = 0; i < m_bifurcations[n].size(); ++i)
1394  {
1395  int 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 (int n = 0; n < m_mergingJcts.size(); ++n, ++cnt)
1404  {
1405  // Merged vessel
1406  for (int 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 (int n = 0; n < m_vesselIntfcs.size(); ++n, ++cnt)
1417  {
1418  for (int 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 (int 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 (int 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 (int 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 (int 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 (int 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 (int 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  int 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  int rank = m_comm->GetRank();
350 
351  // Fill array with domain details
352  int dmax = 0;
353  for (auto &d : m_domain)
354  {
355  dmax = (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 (int 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 (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 (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  int 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  int 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 (rank == sharedproc)
546  {
547  // save all communicators on this for
548  // split comm setup
549  for (auto &d : SharedProc)
550  {
551  if (d.second == 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  int foundallcomms = 0;
589  int i;
590  for (i = 0; i < nprocs; ++i)
591  {
592  if (DoneCreate[i] == false)
593  {
594  break;
595  }
596  }
597 
598  if (i == nprocs)
599  {
600  search = false;
601  }
602  else
603  {
604  commcall++;
605  ASSERTL0(commcall < 4 * commMax,
606  "Failed to find sub communicators "
607  "pattern in 4*commMax searches");
608  }
609  }
610 
611  ASSERTL1(CreateComm.size() == nShrProc,
612  "Have not created communicators for all shared procs");
613 
614  // determine maxmimum number of CreateComm size
615  int maxCreateComm = CreateComm.size();
616  m_comm->AllReduce(maxCreateComm, LibUtilities::ReduceMax);
617 
618  // loop over CreateComm list
619  for (int i = 0; i < maxCreateComm; ++i)
620  {
621  LibUtilities::CommSharedPtr newcomm;
622  int flag = 0;
623 
624  for (auto &d : CreateComm)
625  {
626  if (d.second.first == i)
627  {
628  flag = d.second.second;
629  }
630  }
631 
632  newcomm = m_comm->CommCreateIf(flag);
633 
634  for (auto &d : CreateComm)
635  {
636  if (d.second.first == i)
637  {
638  retval[d.first] = newcomm;
639  }
640  }
641  }
642 
643  // Finally need to re-order domain list so that we call the
644  // communicators in a non-locking manner. This is done by
645  // uniquely numbering each shared proc/comm interface and
646  // ensuring the domains are ordered so that we call this
647  // numbering in an increasing order
648 
649  Array<OneD, NekDouble> numShared(nprocs, 0.0);
650  // determine the number of communicators on this rank where
651  // the share proc has a higher rank id
652  numShared[rank] = nShrProc;
653  m_comm->AllReduce(numShared, LibUtilities::ReduceSum);
654  int totShared = (int)Vmath::Vsum(nprocs, numShared, 1);
655 
656  if (totShared)
657  {
658  int cnt = 0;
659  Array<OneD, NekDouble> numOffset(nprocs, 0.0);
660  for (auto &s : SharedProc)
661  {
662  if (s.second > rank)
663  {
664  cnt++;
665  }
666  }
667 
668  numOffset[rank] = cnt;
669  m_comm->AllReduce(numOffset, LibUtilities::ReduceMax);
670 
671  // make numShared into a cumulative list
672  for (int i = 1; i < nprocs; ++i)
673  {
674  numShared[i] += numShared[i - 1];
675  numOffset[i] += numOffset[i - 1];
676  }
677  for (int i = nprocs - 1; i > 0; --i)
678  {
679  numShared[i] = numShared[i - 1];
680  numOffset[i] = numOffset[i - 1];
681  }
682  numShared[0] = 0;
683  numOffset[0] = 0;
684 
685  Array<OneD, NekDouble> shareddom(totShared, -1.0);
686  cnt = 0; // collect a list of domains that each shared commm is
687  // attached to
688  for (auto &s : SharedProc)
689  {
690  shareddom[numShared[rank] + cnt] = s.first;
691  ++cnt;
692  }
693  m_comm->AllReduce(shareddom, LibUtilities::ReduceMax);
694 
695  // define a numbering scheme
696  Array<OneD, NekDouble> sharedid(totShared, -1.0);
697  cnt = 0;
698  int cnt1 = 0;
699  for (auto &s : SharedProc)
700  {
701  if (s.second > rank)
702  {
703  sharedid[numShared[rank] + cnt] = numOffset[rank] + cnt1;
704 
705  // find shared proc offset by matching the domain ids.
706  int j;
707  for (j = 0; j < maxCreateComm; ++j)
708  {
709  if ((numShared[s.second] + j < totShared) &&
710  (shareddom[numShared[s.second] + j] == s.first))
711  {
712  break;
713  }
714  }
715 
716  sharedid[numShared[s.second] + j] = numOffset[rank] + cnt1;
717  cnt1++;
718  }
719  cnt++;
720  }
721 
722  // now communicate ids to all procesors.
723  m_comm->AllReduce(sharedid, LibUtilities::ReduceMax);
724 
725  if (rank == 0)
726  {
727  for (int i = 0; i < totShared; ++i)
728  {
729  ASSERTL1(sharedid[i] != -1.0,
730  "Failed to number shared proc uniquely");
731  }
732  }
733 
734  cnt = 0;
735  int maxoffset =
736  (int)Vmath::Vmax(nShrProc, &sharedid[numShared[rank]], 1);
737  m_comm->AllReduce(maxoffset, LibUtilities::ReduceMax);
738  maxoffset++;
739 
740  // make a map relating the order of the SharedProc to the domain id
741  map<int, int> ShrToDom;
742  cnt = 0;
743  for (auto &s : SharedProc)
744  {
745  ShrToDom[cnt] = s.first;
746  ++cnt;
747  }
748 
749  // Set up a set the domain ids listed in the ordering of
750  // the SharedProc unique numbering
751  set<int> doneDom;
752  NekDouble maxdom = Vmath::Vmax(totShared, shareddom, 1);
753  maxdom++;
754  for (int i = 0; i < nShrProc; ++i)
755  {
756  int minId =
757  Vmath::Imin(nShrProc, &sharedid[numShared[rank]], 1);
758  sharedid[numShared[rank] + minId] += maxoffset;
759 
760  m_domOrder.push_back(ShrToDom[minId]);
761  doneDom.insert(ShrToDom[minId]);
762  }
763 
764  // add all the other
765  for (auto &d : m_domain)
766  {
767  if (doneDom.count(d.first) == 0)
768  {
769  m_domOrder.push_back(d.first);
770  }
771  }
772  }
773  }
774 }
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:324
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 (int 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 (int 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 (int 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 (int 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 (int 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 1051 of file PulseWaveSystem.cpp.

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

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

777 {
778  map<int, std::vector<InterfacePointShPtr>> VidToDomain;
779 
780  // First set up domain to determine how many vertices meet at a
781  // vertex. This allows us to
782 
783  /* Loop over domain and find out if we have any undefined
784  * boundary conditions representing interfaces. If so make a
785  * map based around vid and storing the domains that are
786  * part of interfaces.
787  */
788  for (int omega = 0; omega < m_nDomains; ++omega)
789  {
790  int vesselID = omega * m_nVariables;
791 
792  for (int i = 0; i < (m_vessels[vesselID]->GetBndConditions()).size();
793  ++i)
794  {
795  if (m_vessels[vesselID]->GetBndConditions()[i]->GetUserDefined() ==
796  "Interface")
797  {
798  // Get Vid of interface
799  int vid = m_vessels[vesselID]
800  ->UpdateBndCondExpansion(i)
801  ->GetExp(0)
802  ->GetGeom()
803  ->GetVid(0);
804 
805  MultiRegions::ExpListSharedPtr trace =
806  m_vessels[vesselID]->GetTrace();
807  InterfacePointShPtr Ipt;
808 
809  bool finish = false;
810  // find which elmt, the local vertex and the data
811  // offset of point
812  for (int n = 0; n < m_vessels[vesselID]->GetExpSize(); ++n)
813  {
814  for (int 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  int 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 (int 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 (int 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 (int 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 (int 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 (int 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 (int 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 (int n = 0; n < m_bifurcations.size(); ++n, ++cnt)
1021  {
1022  int vid = m_bifurcations[n][0]->m_vid;
1023  for (int i = 0; i < 3; ++i)
1024  {
1025  tmp1[3 * cnt + i] = (NekDouble)(3 * vid + i + 1);
1026  }
1027  }
1028 
1029  for (int n = 0; n < m_mergingJcts.size(); ++n, ++cnt)
1030  {
1031  int vid = m_mergingJcts[n][0]->m_vid;
1032  for (int i = 0; i < 3; ++i)
1033  {
1034  int offset = m_mergingJcts[n][i]->m_riemannOrd;
1035  tmp1[3 * cnt + i] = (NekDouble)(3 * vid + i + 1);
1036  }
1037  }
1038 
1039  for (int n = 0; n < m_vesselIntfcs.size(); ++n, ++cnt)
1040  {
1041  int vid = m_vesselIntfcs[n][0]->m_vid;
1042  for (int i = 0; i < 3; ++i)
1043  {
1044  tmp1[3 * cnt + i] = (NekDouble)(3 * vid + i + 1);
1045  }
1046  }
1047 
1048  m_intComm = Gs::Init(tmp1, m_comm->GetRowComm(), verbose);
1049 }
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  )
protectedvirtual

Sets up initial conditions.

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

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1220 of file PulseWaveSystem.cpp.

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

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  )
protectedvirtual

Solves an unsteady problem.

NEEDS Updating:

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

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

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1271 of file PulseWaveSystem.cpp.

1272 {
1273  NekDouble IntegrationTime = 0.0;
1274  int i;
1275  int n;
1276  int nchk = 1;
1277 
1278  Array<OneD, Array<OneD, NekDouble>> fields(m_nVariables);
1279 
1280  for (int i = 0; i < m_nVariables; ++i)
1281  {
1282  fields[i] = m_vessels[i]->UpdatePhys();
1283  m_fields[i]->SetPhysState(false);
1284  }
1285 
1286  m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
1287 
1288  // Time loop
1289  for (n = 0; n < m_steps; ++n)
1290  {
1291  LibUtilities::Timer time_v_IntStep;
1292  time_v_IntStep.Start();
1293  fields = m_intScheme->TimeIntegrate(n, m_timestep, m_ode);
1294  m_time += m_timestep;
1295  time_v_IntStep.Stop();
1296  time_v_IntStep.AccumulateRegion("PulseWaveSystem::TimeIntegrationStep",
1297  0);
1298  // IntegrationTime += timer.TimePerTest(1);
1299 
1300  // Write out status information.
1301  if (m_session->GetComm()->GetRank() == 0 && !((n + 1) % m_infosteps))
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 (!((n + 1) % m_checksteps))
1309  {
1310  for (i = 0; i < m_nVariables; ++i)
1311  {
1312  int cnt = 0;
1313  for (int 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 (int 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:349
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:343
Nektar::SolverUtils::EquationSystem::m_steps
int m_steps
Number of steps to take.
Definition: EquationSystem.h:364
Nektar::SolverUtils::EquationSystem::m_checksteps
int m_checksteps
Number of steps between checkpoints.
Definition: EquationSystem.h:366
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:85
Nektar::SolverUtils::UnsteadySystem::m_intScheme
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
Definition: UnsteadySystem.h:83
Nektar::SolverUtils::UnsteadySystem::m_infosteps
int m_infosteps
Number of time steps between outputting status information.
Definition: UnsteadySystem.h:76
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::UnsteadySystem::m_infosteps, Nektar::SolverUtils::UnsteadySystem::m_intScheme, m_nDomains, m_nVariables, Nektar::SolverUtils::UnsteadySystem::m_ode, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_steps, Nektar::SolverUtils::EquationSystem::m_time, Nektar::SolverUtils::EquationSystem::m_timestep, m_vessels, Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), and Vmath::Vcopy().

◆ v_InitObject()

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

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

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Reimplemented in Nektar::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  int 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  int 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 (int 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 (int 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 (int i = 0; i < (int)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:1051
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:339
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:695
Nektar::SolverUtils::EquationSystem::m_projectionType
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:390
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:681
Nektar::SolverUtils::UnsteadySystem::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true)
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::PulseWaveSystemOutput::v_InitObject(), and Nektar::PulseWavePropagation::v_InitObject().

◆ v_L2Error()

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

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1949 of file PulseWaveSystem.cpp.

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

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 
)
protectedvirtual

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

Compute the error in the L_inf-norm

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 2046 of file PulseWaveSystem.cpp.

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

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  )
protectedvirtual

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

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 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 (int 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  int nFieldDefPerDomain = FieldDef.size() / m_nDomains;
1868  int cnt;
1869 
1870  // Copy Data into FieldData and set variable
1871  for (int n = 0; n < m_nDomains; ++n)
1872  {
1873  // Outputs area and velocity
1874  for (int j = 0; j < m_nVariables; ++j)
1875  {
1876  for (int 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 (int 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 (int 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:333
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:662
Nektar::SolverUtils::EquationSystem::m_fieldMetaDataMap
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
Definition: EquationSystem.h:397

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

int 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

int 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

int 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 Tue Aug 23 2022 07:27:18 for Nektar++ by doxygen 1.9.1