Nektar++
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
Nektar::PulseWaveSystem Class Reference

Base class for unsteady solvers. More...

#include <PulseWaveSystem.h>

Inheritance diagram for Nektar::PulseWaveSystem:
[legend]

Public Member Functions

int GetNdomains ()
 
Array< OneD, MultiRegions::ExpListSharedPtr > UpdateVessels (void)
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT ~UnsteadySystem () override=default
 Destructor.
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Calculate the larger time-step mantaining the problem stable.
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeSharedPtr & GetTimeIntegrationScheme ()
 Returns the time integration scheme.
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeOperators & GetTimeIntegrationSchemeOperators ()
 Returns the time integration scheme operators.
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor.
 
SOLVER_UTILS_EXPORT void InitObject (bool DeclareField=true)
 Initialises the members of this object.
 
SOLVER_UTILS_EXPORT void DoInitialise (bool dumpInitialConditions=true)
 Perform any initialisation necessary before solving the problem.
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem.
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space.
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space.
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve.
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name.
 
template<class T >
std::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name.
 
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession ()
 Get Session name.
 
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure ()
 Get pressure field if available.
 
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.
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda.
 
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction (std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
 Get a SessionFunction by name.
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields.
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution.
 
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.
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields.
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation.
 
SOLVER_UTILS_EXPORT NekDouble H1Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the H1 error between fields and a given exact solution.
 
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].
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields.
 
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.
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields.
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename.
 
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.
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file.
 
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.
 
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.
 
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.
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary.
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMap & UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated.
 
SOLVER_UTILS_EXPORT NekDouble GetTime ()
 Return final time.
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, const Array< OneD, const NekDouble > &input)
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > & UpdatePhysField (const int i)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT int GetInfoSteps ()
 
SOLVER_UTILS_EXPORT void SetInfoSteps (int num)
 
SOLVER_UTILS_EXPORT void SetIterationNumberPIT (int num)
 
SOLVER_UTILS_EXPORT void SetWindowNumberPIT (int num)
 
SOLVER_UTILS_EXPORT Array< OneD, const Array< OneD, NekDouble > > GetTraceNormals ()
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time.
 
SOLVER_UTILS_EXPORT bool NegatedOp ()
 Identify if operator is negated in DoSolve.
 
- Public Member Functions inherited from Nektar::SolverUtils::ALEHelper
virtual ~ALEHelper ()=default
 
virtual SOLVER_UTILS_EXPORT void v_ALEInitObject (int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 
SOLVER_UTILS_EXPORT void InitObject (int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 
virtual SOLVER_UTILS_EXPORT void v_UpdateGridVelocity (const NekDouble &time)
 
virtual SOLVER_UTILS_EXPORT void v_ALEPreMultiplyMass (Array< OneD, Array< OneD, NekDouble > > &fields)
 
SOLVER_UTILS_EXPORT void ALEDoElmtInvMass (Array< OneD, Array< OneD, NekDouble > > &traceNormals, Array< OneD, Array< OneD, NekDouble > > &fields, NekDouble time)
 Update m_fields with u^n by multiplying by inverse mass matrix. That's then used in e.g. checkpoint output and L^2 error calculation.
 
SOLVER_UTILS_EXPORT void ALEDoElmtInvMassBwdTrans (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
SOLVER_UTILS_EXPORT void MoveMesh (const NekDouble &time, Array< OneD, Array< OneD, NekDouble > > &traceNormals)
 
SOLVER_UTILS_EXPORT void ResetMatricesNormal (Array< OneD, Array< OneD, NekDouble > > &traceNormals)
 
SOLVER_UTILS_EXPORT void UpdateNormalsFlag ()
 
const Array< OneD, const Array< OneD, NekDouble > > & GetGridVelocity ()
 
bool & GetUpdateNormalsFlag ()
 
SOLVER_UTILS_EXPORT const Array< OneD, const Array< OneD, NekDouble > > & GetGridVelocityTrace ()
 
SOLVER_UTILS_EXPORT void ExtraFldOutputGridVelocity (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 
SOLVER_UTILS_EXPORT void ExtraFldOutputGrid (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 

Protected Member Functions

 PulseWaveSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises PulseWaveSystem class members.
 
 ~PulseWaveSystem () override=default
 
void v_InitObject (bool DeclareField=false) override
 
void v_DoInitialise (bool dumpInitialConditions=false) override
 Sets up initial conditions.
 
void v_DoSolve () override
 Solves an unsteady problem.
 
void LinkSubdomains (Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &fields)
 Links the subdomains.
 
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.
 
void MergingRiemann (Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0, Array< OneD, NekDouble > &alpha)
 Riemann Problem for Merging Flow.
 
void InterfaceRiemann (Array< OneD, NekDouble > &Au, Array< OneD, NekDouble > &uu, Array< OneD, NekDouble > &beta, Array< OneD, NekDouble > &A_0, Array< OneD, NekDouble > &alpha)
 Riemann Problem for Interface/Junction.
 
void v_Output (void) override
 
void CheckPoint_Output (const int n)
 
NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false) override
 Compute the L2 error between fields and a given exact solution.
 
NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray) override
 Compute the L_inf error between fields and a given exact solution.
 
void WriteVessels (const std::string &outname)
 Write input fields to the given filename.
 
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.
 
SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class.
 
SOLVER_UTILS_EXPORT void v_DoSolve () override
 Solves an unsteady problem.
 
virtual SOLVER_UTILS_EXPORT void v_PrintStatusInformation (const int step, const NekDouble cpuTime)
 Print Status Information.
 
virtual SOLVER_UTILS_EXPORT void v_PrintSummaryStatistics (const NekDouble intTime)
 Print Summary Statistics.
 
SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true) override
 Sets up initial conditions.
 
SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &s) override
 Print a summary of time stepping parameters.
 
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.
 
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_RequireFwdTrans ()
 
virtual SOLVER_UTILS_EXPORT void v_SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
virtual SOLVER_UTILS_EXPORT bool v_UpdateTimeStepCheck ()
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control.
 
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.
 
SOLVER_UTILS_EXPORT void DoDummyProjection (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
 Perform dummy projection.
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members.
 
virtual SOLVER_UTILS_EXPORT NekDouble v_H1Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the H_1 error computation between fields and a given exact solution.
 
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys ()
 Virtual function for transformation to physical space.
 
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff ()
 Virtual function for transformation to coefficient space.
 
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 bool v_NegatedOp (void)
 Virtual function to identify if operator is negated in DoSolve.
 
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 

Protected Attributes

Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
 
size_t m_nDomains
 
size_t m_currentDomain
 
size_t m_nVariables
 
UpwindTypePulse m_upwindTypePulse
 
Array< OneD, int > m_fieldPhysOffset
 
NekDouble m_rho
 
NekDouble m_pext
 
NekDouble m_C
 
NekDouble m_RT
 
NekDouble m_pout
 
Array< OneD, Array< OneD, NekDouble > > m_A_0
 
Array< OneD, Array< OneD, NekDouble > > m_A_0_trace
 
Array< OneD, Array< OneD, NekDouble > > m_beta
 
Array< OneD, Array< OneD, NekDouble > > m_beta_trace
 
Array< OneD, Array< OneD, NekDouble > > m_gamma
 
Array< OneD, Array< OneD, NekDouble > > m_alpha
 
Array< OneD, Array< OneD, NekDouble > > m_alpha_trace
 
Array< OneD, Array< OneD, NekDouble > > m_trace_fwd_normal
 
std::map< int, SpatialDomains::CompositeMap > m_domain
 
std::vector< int > m_domOrder
 
std::map< int, std::vector< int > > m_domainToFilterIDs
 
std::map< int, int > m_filterToVesselID
 
Array< OneD, Array< OneD, NekDouble > > m_pressure
 
PulseWavePressureAreaSharedPtr m_pressureArea
 
bool extraFields = false
 
Array< OneD, Array< OneD, NekDouble > > m_PWV
 
Array< OneD, Array< OneD, NekDouble > > m_W1
 
Array< OneD, Array< OneD, NekDouble > > m_W2
 
std::vector< std::vector< InterfacePointShPtr > > m_vesselIntfcs
 
std::vector< std::vector< InterfacePointShPtr > > m_bifurcations
 
std::vector< std::vector< InterfacePointShPtr > > m_mergingJcts
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
 Wrapper to the time integration scheme.
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use.
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check.
 
std::vector< int > m_intVariables
 
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1).
 
NekDouble m_CFLGrowth
 CFL growth rate.
 
NekDouble m_CFLEnd
 Maximun cfl in cfl growth.
 
int m_abortSteps
 Number of steps between checks for abort conditions.
 
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used.
 
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used.
 
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used.
 
int m_steadyStateSteps
 Check for steady state at step interval.
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at.
 
int m_filtersInfosteps
 Number of time steps between outputting filters information.
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state.
 
std::ofstream m_errFile
 
NekDouble m_epsilon
 Diffusion coefficient.
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator.
 
bool m_verbose
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader.
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtr > m_sessionFunctions
 Map of known SessionFunctions.
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output.
 
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
 Array holding all dependent variables.
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object.
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh.
 
std::string m_sessionName
 Name of the session.
 
NekDouble m_time
 Current time of simulation.
 
int m_initialStep
 Number of the step where the simulation should begin.
 
NekDouble m_fintime
 Finish time of the simulation.
 
NekDouble m_timestep
 Time step size.
 
NekDouble m_lambda
 Lambda constant in real system if one required.
 
NekDouble m_checktime
 Time between checkpoints.
 
NekDouble m_lastCheckTime
 
NekDouble m_TimeIncrementFactor
 
int m_nchk
 Number of checkpoints written so far.
 
int m_steps
 Number of steps to take.
 
int m_checksteps
 Number of steps between checkpoints.
 
int m_infosteps
 Number of time steps between outputting status information.
 
int m_iterPIT = 0
 Number of parallel-in-time time iteration.
 
int m_windowPIT = 0
 Index of windows for parallel-in-time time iteration.
 
int m_spacedim
 Spatial dimension (>= expansion dim).
 
int m_expdim
 Expansion dimension.
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used.
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used.
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used.
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform.
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations.
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous.
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction.
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity.
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields.
 
Array< OneD, NekDouble > m_movingFrameData
 Moving reference frame status in the inertial frame X, Y, Z, Theta_x, Theta_y, Theta_z, U, V, W, Omega_x, Omega_y, Omega_z, A_x, A_y, A_z, DOmega_x, DOmega_y, DOmega_z, pivot_x, pivot_y, pivot_z.
 
std::vector< std::string > m_strFrameData
 variable name in m_movingFrameData
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error.
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous)
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous)
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous)
 
int m_npointsX
 number of points in X direction (if homogeneous)
 
int m_npointsY
 number of points in Y direction (if homogeneous)
 
int m_npointsZ
 number of points in Z direction (if homogeneous)
 
int m_HomoDirec
 number of homogenous directions
 
- Protected Attributes inherited from Nektar::SolverUtils::ALEHelper
Array< OneD, MultiRegions::ExpListSharedPtr > m_fieldsALE
 
Array< OneD, Array< OneD, NekDouble > > m_gridVelocity
 
Array< OneD, Array< OneD, NekDouble > > m_gridVelocityTrace
 
std::vector< ALEBaseShPtr > m_ALEs
 
bool m_ALESolver = false
 
bool m_meshDistorted = false
 
bool m_implicitALESolver = false
 
bool m_updateNormals = false
 
NekDouble m_prevStageTime = 0.0
 
int m_spaceDim
 

Private Member Functions

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

Private Attributes

Gs::gs_data * m_intComm
 Communicator for interfaces.
 

Additional Inherited Members

- Static Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
static std::string cmdSetStartTime
 
static std::string cmdSetStartChkNum
 
- Protected Types inherited from Nektar::SolverUtils::EquationSystem
enum  HomogeneousType { eHomogeneous1D , eHomogeneous2D , eHomogeneous3D , eNotHomogeneous }
 Parameter for homogeneous expansions. More...
 
- Static Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
static std::string equationSystemTypeLookupIds []
 
static std::string projectionTypeLookupIds []
 

Detailed Description

Base class for unsteady solvers.

Initialises the arterial subdomains in m_vessels and sets up all domain-linking conditions (bifurcations, intefaces/junctions, merging flows). Detects the network structure and assigns boundary conditons. Also provides the underlying timestepping framework for pulse wave solvers including the general timestepping routines.

Definition at line 77 of file PulseWaveSystem.h.

Constructor & Destructor Documentation

◆ PulseWaveSystem()

Nektar::PulseWaveSystem::PulseWaveSystem ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const SpatialDomains::MeshGraphSharedPtr &  pGraph 
)
protected

Initialises PulseWaveSystem class members.

Processes SolverInfo parameters from the session file and sets up timestepping-specific code.

Parameters
m_SessionSession object to read parameters from.

Definition at line 64 of file PulseWaveSystem.cpp.

67 : UnsteadySystem(pSession, pGraph)
68{
69}
Nektar::SolverUtils::UnsteadySystem
Base class for unsteady solvers.
Definition UnsteadySystem.h:47

◆ ~PulseWaveSystem()

Nektar::PulseWaveSystem::~PulseWaveSystem ( )
overrideprotecteddefault

Member Function Documentation

◆ BifurcationRiemann()

void Nektar::PulseWaveSystem::BifurcationRiemann ( Array< OneD, NekDouble > &  Au,
Array< OneD, NekDouble > &  uu,
Array< OneD, NekDouble > &  beta,
Array< OneD, NekDouble > &  A_0,
Array< OneD, NekDouble > &  alpha 
)
protected

Riemann Problem for Bifurcation.

Solves the Riemann problem at a bifurcation by assuming subsonic flow at both sides of the boundary and by applying conservation of mass and continuity of the total pressure \( \frac{p}{rho} + \frac{u^{2}}{2}. \) The other 3 missing equations come from the characteristic variables. For further information see "Pulse WavePropagation in the human vascular system" Section 3.4.4

Definition at line 1588 of file PulseWaveSystem.cpp.

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

1919{
1920 std::stringstream outname;
1921 outname << m_sessionName << "_" << n << ".chk";
1922
1923 WriteVessels(outname.str());
1924}
Nektar::PulseWaveSystem::WriteVessels
void WriteVessels(const std::string &outname)
Write input fields to the given filename.
Definition PulseWaveSystem.cpp:1930
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition EquationSystem.h:473

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

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

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

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

◆ FillDataFromInterfacePoint()

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

Definition at line 1412 of file PulseWaveSystem.cpp.

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

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

Referenced by EnforceInterfaceConditions().

◆ GetCommArray()

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

Set and retrn a series of communicators for each partition.

Definition at line 347 of file PulseWaveSystem.cpp.

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

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

81 {
82 return m_nDomains;
83 }

References m_nDomains.

◆ InterfaceRiemann()

void Nektar::PulseWaveSystem::InterfaceRiemann ( Array< OneD, NekDouble > &  Au,
Array< OneD, NekDouble > &  uu,
Array< OneD, NekDouble > &  beta,
Array< OneD, NekDouble > &  A_0,
Array< OneD, NekDouble > &  alpha 
)
protected

Riemann Problem for Interface/Junction.

Solves the Riemann problem at an interdomain junction/Interface by assuming subsonic flow at both sides of the boundary and by applying conservation of mass and continuity of the total pressure \( \frac{p}{rho} + \frac{u^{2}}{2}. \) The other 2 missing equations come from the characteristic variables. For further information see "Pulse WavePropagation in the human vascular system" Section 3.4.

Definition at line 1815 of file PulseWaveSystem.cpp.

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

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

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

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

Referenced by EnforceInterfaceConditions().

◆ SetUpDomainInterfaceBCs()

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

Definition at line 1074 of file PulseWaveSystem.cpp.

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

References Nektar::SpatialDomains::BoundaryConditions::AddBoundaryConditions(), Nektar::SpatialDomains::BoundaryConditions::AddBoundaryRegions(), Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), Nektar::SpatialDomains::Geometry::GetGlobalID(), Nektar::SpatialDomains::Geometry::GetVertex(), m_domain, m_domOrder, Nektar::SolverUtils::EquationSystem::m_session, Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceSum, and Vmath::Vsum().

Referenced by v_InitObject().

◆ SetUpDomainInterfaces()

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

Definition at line 799 of file PulseWaveSystem.cpp.

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

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, and Vmath::Zero().

Referenced by v_InitObject().

◆ UpdateVessels()

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

Definition at line 85 of file PulseWaveSystem.h.

86 {
87 return m_vessels;
88 }

References m_vessels.

◆ v_DoInitialise()

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

Sets up initial conditions.

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1241 of file PulseWaveSystem.cpp.

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

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

◆ v_DoSolve()

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

Solves an unsteady problem.

NEEDS Updating:

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

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1310 of file PulseWaveSystem.cpp.

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

References Nektar::LibUtilities::Timer::AccumulateRegion(), CheckPoint_Output(), Nektar::SolverUtils::EquationSystem::m_checksteps, m_domainToFilterIDs, m_fieldPhysOffset, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::UnsteadySystem::m_filters, m_filterToVesselID, Nektar::SolverUtils::EquationSystem::m_infosteps, Nektar::SolverUtils::UnsteadySystem::m_intScheme, m_nDomains, m_nVariables, Nektar::SolverUtils::UnsteadySystem::m_ode, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_steps, Nektar::SolverUtils::EquationSystem::m_time, Nektar::SolverUtils::EquationSystem::m_timestep, m_vessels, Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), and Vmath::Vcopy().

◆ v_InitObject()

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

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 80 of file PulseWaveSystem.cpp.

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

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

◆ v_L2Error()

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

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 2024 of file PulseWaveSystem.cpp.

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

References ASSERTL0, Nektar::SolverUtils::EquationSystem::EvaluateExactSolution(), Nektar::SolverUtils::EquationSystem::GetFunction(), Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::SolverUtils::EquationSystem::m_comm, Nektar::SolverUtils::EquationSystem::m_fields, m_nDomains, Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError, m_nVariables, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_time, m_vessels, Nektar::LibUtilities::ReduceSum, and tinysimd::sqrt().

◆ v_LinfError()

NekDouble Nektar::PulseWaveSystem::v_LinfError ( unsigned int  field,
const Array< OneD, NekDouble > &  exact = NullNekDouble1DArray 
)
overrideprotectedvirtual

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

Compute the error in the L_inf-norm

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 2121 of file PulseWaveSystem.cpp.

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

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

◆ v_Output()

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

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

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1903 of file PulseWaveSystem.cpp.

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

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

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

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

Referenced by CheckPoint_Output(), and v_Output().

Member Data Documentation

◆ extraFields

bool Nektar::PulseWaveSystem::extraFields = false
protected

Definition at line 124 of file PulseWaveSystem.h.

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

◆ m_A_0

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

Definition at line 105 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 106 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 110 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 111 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 107 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 108 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 130 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_C

NekDouble Nektar::PulseWaveSystem::m_C
protected

Definition at line 101 of file PulseWaveSystem.h.

◆ m_currentDomain

size_t Nektar::PulseWaveSystem::m_currentDomain
protected

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

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

◆ m_domainToFilterIDs

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

Definition at line 117 of file PulseWaveSystem.h.

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

◆ m_domOrder

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

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

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

◆ m_filterToVesselID

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

Definition at line 118 of file PulseWaveSystem.h.

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

◆ m_gamma

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

Definition at line 109 of file PulseWaveSystem.h.

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

◆ m_intComm

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

Communicator for interfaces.

Definition at line 196 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_mergingJcts

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

Definition at line 131 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_nDomains

size_t Nektar::PulseWaveSystem::m_nDomains
protected

Definition at line 92 of file PulseWaveSystem.h.

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

◆ m_nVariables

size_t Nektar::PulseWaveSystem::m_nVariables
protected

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

Referenced by v_InitObject().

◆ m_pout

NekDouble Nektar::PulseWaveSystem::m_pout
protected

Definition at line 103 of file PulseWaveSystem.h.

◆ m_pressure

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

Definition at line 120 of file PulseWaveSystem.h.

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

◆ m_pressureArea

PulseWavePressureAreaSharedPtr Nektar::PulseWaveSystem::m_pressureArea
protected

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

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

◆ m_rho

NekDouble Nektar::PulseWaveSystem::m_rho
protected

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

◆ m_trace_fwd_normal

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

Definition at line 112 of file PulseWaveSystem.h.

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

◆ m_upwindTypePulse

UpwindTypePulse Nektar::PulseWaveSystem::m_upwindTypePulse
protected

Definition at line 95 of file PulseWaveSystem.h.

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

◆ m_vesselIntfcs

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

Definition at line 129 of file PulseWaveSystem.h.

Referenced by EnforceInterfaceConditions(), and SetUpDomainInterfaces().

◆ m_vessels

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

Definition at line 91 of file PulseWaveSystem.h.

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

◆ m_W1

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

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

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

Generated on Sun Nov 9 2025 20:17:03 for Nektar++ by doxygen 1.9.8