Nektar++
Public Member Functions | Static Public Member Functions | Public Attributes | Static Public Attributes | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | Friends | List of all members
Nektar::MMFSWE Class Reference

#include <MMFSWE.h>

Inheritance diagram for Nektar::MMFSWE:
[legend]

Public Member Functions

virtual ~MMFSWE ()
 Destructor. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::MMFSystem
SOLVER_UTILS_EXPORT MMFSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 
virtual SOLVER_UTILS_EXPORT ~MMFSystem ()
 
SOLVER_UTILS_EXPORT void MMFInitObject (const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem=-1)
 
SOLVER_UTILS_EXPORT void CopyBoundaryTrace (const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const BoundaryCopyType BDCopyType, const int var=0, const std::string btype="NoUserDefined")
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
virtual SOLVER_UTILS_EXPORT ~UnsteadySystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble >> &inarray)
 Calculate the larger time-step mantaining the problem stable. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void SetUpTraceNormals (void)
 
SOLVER_UTILS_EXPORT void InitObject ()
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise ()
 Perform any initialisation necessary before solving the problem. More...
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem. More...
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space. More...
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space. More...
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve. More...
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name. More...
 
template<class T >
std::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name. More...
 
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession ()
 Get Session name. More...
 
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure ()
 Get pressure field if available. More...
 
SOLVER_UTILS_EXPORT void ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 
SOLVER_UTILS_EXPORT void PrintSummary (std::ostream &out)
 Print a summary of parameters and solver characteristics. More...
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda. More...
 
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction (std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
 Get a SessionFunction by name. More...
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields. More...
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the L2 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > ErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write checkpoint file of custom data fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write input fields to the given filename. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
 Input field data from the given file to multiple domains. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, std::vector< std::string > &fieldStr, Array< OneD, Array< OneD, NekDouble > > &coeffs)
 Output a field. Input field data into array from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, MultiRegions::ExpListSharedPtr &pField, std::string &pFieldName)
 Output a field. Input field data into ExpList from the given file. More...
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary. More...
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMap & UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetFinalTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp ()
 Virtual function to identify if operator is negated in DoSolve. More...
 

Static Public Member Functions

static SolverUtils::EquationSystemSharedPtr create (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Creates an instance of this class. More...
 

Public Attributes

TestType m_TestType
 
- Public Attributes inherited from Nektar::SolverUtils::MMFSystem
NekDouble m_pi
 
int m_shapedim
 
SurfaceType m_surfaceType
 
UpwindType m_upwindType
 
TestMaxwellType m_TestMaxwellType
 
PolType m_PolType
 
IncType m_IncType
 
Array< OneD, NekDouble > m_MMFfactors
 
- Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 

Static Public Attributes

static std::string className
 Name of class. More...
 

Protected Member Functions

 MMFSWE (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Session reader. More...
 
void DoOdeRhs (const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
 Compute the RHS. More...
 
void DoOdeProjection (const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
 Compute the projection. More...
 
void SteadyZonalFlow (unsigned int field, Array< OneD, NekDouble > &outfield)
 
void UnsteadyZonalFlow (unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
 
void IsolatedMountainFlow (unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
 
void UnstableJetFlow (unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
 
void RossbyWave (unsigned int field, Array< OneD, NekDouble > &outfield)
 
NekDouble ComputeUnstableJetEta (const NekDouble theta)
 
NekDouble ComputeUnstableJetuphi (const NekDouble theta)
 
NekDouble ComputeMass (const Array< OneD, const NekDouble > &eta)
 
NekDouble ComputeEnergy (const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
 
NekDouble ComputeEnstrophy (const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
 
void ComputeVorticity (const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
 
void GetNormalVelocity (Array< OneD, NekDouble > &traceVn)
 Get the normal velocity. More...
 
void ComputeNablaCdotVelocity (Array< OneD, NekDouble > &vellc)
 
void PrimitiveToConservative (const Array< OneD, const Array< OneD, NekDouble >> &physin, Array< OneD, Array< OneD, NekDouble >> &physout)
 
void ConservativeToPrimitive (const Array< OneD, const Array< OneD, NekDouble >> &physin, Array< OneD, Array< OneD, NekDouble >> &physout)
 
void PrimitiveToConservative ()
 
void ConservativeToPrimitive ()
 
void WeakDGSWEDirDeriv (const Array< OneD, Array< OneD, NekDouble >> &InField, Array< OneD, Array< OneD, NekDouble >> &OutField)
 
void AddDivForGradient (Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
 
void NumericalSWEFlux (Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
 
void GetSWEFluxVector (const int i, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
 
void RiemannSolverHLLC (const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
 
void Computehhuhvflux (NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, NekDouble hstar, NekDouble &hflux, NekDouble &huflux, NekDouble &hvflux)
 
void AverageFlux (const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
 
void LaxFriedrichFlux (const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
 
void RusanovFlux (const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
 
void ComputeMagAndDot (const int index, NekDouble &MageF1, NekDouble &MageF2, NekDouble &MageB1, NekDouble &MageB2, NekDouble &eF1_cdot_eB1, NekDouble &eF1_cdot_eB2, NekDouble &eF2_cdot_eB1, NekDouble &eF2_cdot_eB2)
 
void AddCoriolis (Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
 
void AddElevationEffect (Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
 
void AddRotation (Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
 
void Compute_demdt_cdot_ek (const int indm, const int indk, const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, NekDouble > &outarray)
 
void TestVorticityComputation ()
 
void EvaluateWaterDepth (void)
 
void EvaluateCoriolis (void)
 
void EvaluateCoriolisForZonalFlow (Array< OneD, NekDouble > &outarray)
 
void EvaluateStandardCoriolis (Array< OneD, NekDouble > &outarray)
 
void SetBoundaryConditions (Array< OneD, Array< OneD, NekDouble >> &inarray, NekDouble time)
 
void WallBoundary2D (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble >> &physarray)
 
virtual void v_InitObject ()
 Initialise the object. More...
 
virtual void v_DoSolve ()
 Solves an unsteady problem. More...
 
virtual void v_DoInitialise ()
 Sets up initial conditions. More...
 
virtual void v_GenerateSummary (SolverUtils::SummaryList &s)
 Print Summary. More...
 
virtual void v_SetInitialConditions (const NekDouble initialtime, bool dumpInitialConditions, const int domain)
 
virtual void v_EvaluateExactSolution (unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 
virtual NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised)
 Virtual function for the L_2 error computation between fields and a given exact solution. More...
 
virtual NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln)
 Virtual function for the L_inf error computation between fields and a given exact solution. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::MMFSystem
void SetUpMovingFrames (const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem)
 
void CheckMovingFrames (const Array< OneD, const Array< OneD, NekDouble >> &movingframes)
 
SOLVER_UTILS_EXPORT void ComputencdotMF ()
 
SOLVER_UTILS_EXPORT void ComputeDivCurlMF ()
 
SOLVER_UTILS_EXPORT void ComputeMFtrace ()
 
SOLVER_UTILS_EXPORT void VectorDotProd (const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, NekDouble > &v3)
 
SOLVER_UTILS_EXPORT void VectorCrossProd (const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, Array< OneD, NekDouble >> &v3)
 
SOLVER_UTILS_EXPORT void VectorCrossProd (const Array< OneD, NekDouble > &v1, const Array< OneD, NekDouble > &v2, Array< OneD, NekDouble > &v3)
 
SOLVER_UTILS_EXPORT void ComputeCurl (const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes (const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
 
SOLVER_UTILS_EXPORT void DeriveCrossProductMF (Array< OneD, Array< OneD, NekDouble >> &CrossProductMF)
 
SOLVER_UTILS_EXPORT void ComputeNtimesMF ()
 
SOLVER_UTILS_EXPORT void ComputeNtimesFz (const int dir, const Array< OneD, Array< OneD, NekDouble >> &Fwd, const Array< OneD, Array< OneD, NekDouble >> &Bwd, const Array< OneD, const NekDouble > &imFwd, const Array< OneD, const NekDouble > &imBwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
 
SOLVER_UTILS_EXPORT void ComputeNtimesF12 (const Array< OneD, Array< OneD, NekDouble >> &Fwd, const Array< OneD, Array< OneD, NekDouble >> &Bwd, const Array< OneD, const NekDouble > &im1Fwd, const Array< OneD, const NekDouble > &im1Bwd, const Array< OneD, const NekDouble > &im2Fwd, const Array< OneD, const NekDouble > &im2Bwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
 
SOLVER_UTILS_EXPORT void ComputeNtimestimesdFz (const int dir, const Array< OneD, Array< OneD, NekDouble >> &Fwd, const Array< OneD, Array< OneD, NekDouble >> &Bwd, const Array< OneD, const NekDouble > &imFwd, const Array< OneD, const NekDouble > &imBwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
 
SOLVER_UTILS_EXPORT void ComputeNtimestimesdF12 (const Array< OneD, Array< OneD, NekDouble >> &Fwd, const Array< OneD, Array< OneD, NekDouble >> &Bwd, const Array< OneD, const NekDouble > &im1Fwd, const Array< OneD, const NekDouble > &im1Bwd, const Array< OneD, const NekDouble > &im2Fwd, const Array< OneD, const NekDouble > &im2Bwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
 
SOLVER_UTILS_EXPORT void CartesianToSpherical (const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
 
SOLVER_UTILS_EXPORT void ComputeZimYim (Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
 
SOLVER_UTILS_EXPORT void AdddedtMaxwell (const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
 
SOLVER_UTILS_EXPORT void GetMaxwellFluxVector (const int var, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
 
SOLVER_UTILS_EXPORT void GetMaxwellFlux1D (const int var, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
 
SOLVER_UTILS_EXPORT void GetMaxwellFlux2D (const int var, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
 
SOLVER_UTILS_EXPORT void LaxFriedrichMaxwellFlux1D (Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
 
SOLVER_UTILS_EXPORT void UpwindMaxwellFlux1D (Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
 
SOLVER_UTILS_EXPORT void AverageMaxwellFlux1D (Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
 
SOLVER_UTILS_EXPORT void NumericalMaxwellFlux (Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd, const NekDouble time=0.0)
 
SOLVER_UTILS_EXPORT void NumericalMaxwellFluxTM (Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd, const NekDouble time)
 
SOLVER_UTILS_EXPORT void NumericalMaxwellFluxTE (Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd, const NekDouble time)
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetIncidentField (const int var, const NekDouble time)
 
SOLVER_UTILS_EXPORT void Computedemdxicdote ()
 
SOLVER_UTILS_EXPORT NekDouble AvgInt (const Array< OneD, const NekDouble > &inarray)
 
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt (const Array< OneD, const NekDouble > &inarray)
 
SOLVER_UTILS_EXPORT NekDouble AbsIntegral (const Array< OneD, const NekDouble > &inarray)
 
SOLVER_UTILS_EXPORT NekDouble RootMeanSquare (const Array< OneD, const NekDouble > &inarray)
 
SOLVER_UTILS_EXPORT NekDouble VectorAvgMagnitude (const Array< OneD, const Array< OneD, NekDouble >> &inarray)
 
SOLVER_UTILS_EXPORT void GramSchumitz (const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, Array< OneD, NekDouble >> &outarray, bool KeepTheMagnitude=true)
 
SOLVER_UTILS_EXPORT void BubbleSort (Array< OneD, NekDouble > &refarray, Array< OneD, NekDouble > &sortarray)
 
- Protected Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises UnsteadySystem class members. More...
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. More...
 
virtual SOLVER_UTILS_EXPORT void v_AppendOutput1D (Array< OneD, Array< OneD, NekDouble >> &solution1D)
 Print the solution at each solution point in a txt file. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_GetTimeStep (const Array< OneD, const Array< OneD, NekDouble >> &inarray)
 Return the timestep to be used for the next step in the time-marching loop. More...
 
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_RequireFwdTrans ()
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time, int &nchk)
 
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff (const Array< OneD, Array< OneD, NekDouble >> vel, StdRegions::VarCoeffMap &varCoeffMap)
 Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h time velocity. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys ()
 Virtual function for transformation to physical space. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff ()
 Virtual function for transformation to coefficient space. More...
 
virtual SOLVER_UTILS_EXPORT void v_Output (void)
 
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure (void)
 
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 

Protected Attributes

Array< OneD, NekDouble > m_depth
 Still water depth. More...
 
Array< OneD, Array< OneD, NekDouble > > m_Derivdepth
 
Array< OneD, NekDouble > m_coriolis
 Coriolis force. More...
 
int m_AddCoriolis
 
int m_AddRotation
 
NekDouble m_g
 
NekDouble m_alpha
 
NekDouble m_u0
 
NekDouble m_Omega
 
NekDouble m_H0
 
NekDouble m_Hvar
 
NekDouble m_k2
 
NekDouble m_hs0
 
NekDouble m_uthetamax
 
NekDouble m_theta0
 
NekDouble m_theta1
 
NekDouble m_en
 
NekDouble m_hbar
 
NekDouble m_angfreq
 
NekDouble m_K
 
NekDouble m_Vorticity0
 
NekDouble m_Mass0
 
NekDouble m_Energy0
 
NekDouble m_Enstrophy0
 
bool m_primitive
 Indicates if variables are primitive or conservative. More...
 
Array< OneD, Array< OneD, NekDouble > > m_velocity
 Advection velocity. More...
 
Array< OneD, NekDouble > m_traceVn
 
Array< OneD, Array< OneD, NekDouble > > m_veldotMF
 
Array< OneD, NekDouble > m_vellc
 
int m_planeNumber
 
- Protected Attributes inherited from Nektar::SolverUtils::MMFSystem
NekDouble m_alpha
 
NekDouble m_Incfreq
 
int m_SmoothFactor
 
NekDouble m_SFinit
 
Array< OneD, Array< OneD, NekDouble > > m_movingframes
 
Array< OneD, Array< OneD, NekDouble > > m_surfaceNormal
 
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
 
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
 
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
 
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFBwd
 
Array< OneD, Array< OneD, NekDouble > > m_DivMF
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFFwd
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFBwd
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFFwd
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFBwd
 
Array< OneD, Array< OneD, NekDouble > > m_ZimFwd
 
Array< OneD, Array< OneD, NekDouble > > m_ZimBwd
 
Array< OneD, Array< OneD, NekDouble > > m_YimFwd
 
Array< OneD, Array< OneD, NekDouble > > m_YimBwd
 
Array< OneD, Array< OneD, NekDouble > > m_epsvec
 
Array< OneD, Array< OneD, NekDouble > > m_muvec
 
Array< OneD, Array< OneD, NekDouble > > m_negepsvecminus1
 
Array< OneD, Array< OneD, NekDouble > > m_negmuvecminus1
 
Array< OneD, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > > m_dedxi_cdot_e
 
SpatialDomains::GeomMMF m_MMFdir
 
Array< OneD, NekDouble > m_MFlength
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
int m_infosteps
 Number of time steps between outputting status information. More...
 
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
int m_nanSteps
 
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
 
NekDouble m_epsilon
 
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used. More...
 
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used. More...
 
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used. More...
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check. More...
 
std::ofstream m_errFile
 
std::vector< int > m_intVariables
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
NekDouble m_filterTimeWarning
 Number of time steps between outputting status information. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtr > m_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error. More...
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous) More...
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous) More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
int m_HomoDirec
 number of homogenous directions More...
 

Private Member Functions

void TestSWE2Dproblem (const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
 
void Checkpoint_Output_Cartesian (std::string outname)
 

Private Attributes

int m_RossbyDisturbance
 
int m_PurturbedJet
 

Friends

class MemoryManager< MMFSWE >
 

Additional Inherited Members

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

Detailed Description

Definition at line 61 of file MMFSWE.h.

Constructor & Destructor Documentation

◆ ~MMFSWE()

Nektar::MMFSWE::~MMFSWE ( )
virtual

Destructor.

Unsteady linear advection equation destructor.

Definition at line 244 of file MMFSWE.cpp.

245 {
246 }

◆ MMFSWE()

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

Session reader.

Definition at line 51 of file MMFSWE.cpp.

References m_planeNumber.

53  : UnsteadySystem(pSession, pGraph), MMFSystem(pSession, pGraph)
54 {
55  m_planeNumber = 0;
56 }
Nektar::MMFSWE::m_planeNumber
int m_planeNumber
Definition: MMFSWE.h:115
Nektar::SolverUtils::UnsteadySystem::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.
Definition: UnsteadySystem.cpp:68
Nektar::SolverUtils::MMFSystem::MMFSystem
SOLVER_UTILS_EXPORT MMFSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition: MMFSystem.cpp:43

Member Function Documentation

◆ AddCoriolis()

void Nektar::MMFSWE::AddCoriolis ( Array< OneD, Array< OneD, NekDouble >> &  physarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray 
)
protected

Definition at line 1199 of file MMFSWE.cpp.

References m_coriolis, m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Neg(), Vmath::Vadd(), and Vmath::Vmul().

Referenced by DoOdeRhs().

1201 {
1202  int ncoeffs = outarray[0].num_elements();
1203  int nq = physarray[0].num_elements();
1204 
1205  Array<OneD, NekDouble> h(nq);
1206  Array<OneD, NekDouble> tmp(nq);
1207  Array<OneD, NekDouble> tmpc(ncoeffs);
1208 
1209  // physarray is primitive
1210  // conservative formulation compute h
1211  // h = \eta + d
1212  Vmath::Vadd(nq, physarray[0], 1, m_depth, 1, h, 1);
1213 
1214  int indx = 0;
1215  for (int j = 0; j < m_shapedim; ++j)
1216  {
1217  if (j == 0)
1218  {
1219  indx = 2;
1220  }
1221 
1222  else if (j == 1)
1223  {
1224  indx = 1;
1225  }
1226 
1227  // add to hu equation
1228  Vmath::Vmul(nq, m_coriolis, 1, physarray[indx], 1, tmp, 1);
1229  Vmath::Vmul(nq, h, 1, tmp, 1, tmp, 1);
1230 
1231  if (j == 1)
1232  {
1233  Vmath::Neg(nq, tmp, 1);
1234  }
1235 
1236  // N \cdot (e^1 \times e^2 )
1237  // Vmath::Vmul(nq, &m_MF1crossMF2dotSN[0], 1, &tmp[0], 1, &tmp[0], 1);
1238  m_fields[0]->IProductWRTBase(tmp, tmpc);
1239  Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1240  }
1241 }
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::MMFSWE::m_coriolis
Array< OneD, NekDouble > m_coriolis
Coriolis force.
Definition: MMFSWE.h:90
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ AddDivForGradient()

void Nektar::MMFSWE::AddDivForGradient ( Array< OneD, Array< OneD, NekDouble >> &  physarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray 
)
protected

Definition at line 550 of file MMFSWE.cpp.

References GetSWEFluxVector(), Nektar::SolverUtils::MMFSystem::m_DivMF, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Neg(), Vmath::Vadd(), and Vmath::Vmul().

Referenced by DoOdeRhs().

552 {
553  // routine works for both primitive and conservative formulations
554  int ncoeffs = outarray[0].num_elements();
555  int nq = physarray[0].num_elements();
556 
557  Array<OneD, NekDouble> h(nq);
558  Array<OneD, NekDouble> tmp(nq);
559  Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
560  for (int i = 0; i < m_shapedim; ++i)
561  {
562  fluxvector[i] = Array<OneD, NekDouble>(nq);
563  }
564 
565  // Get 0.5 g H*H / || e^m ||^2
566  GetSWEFluxVector(3, physarray, fluxvector);
567 
568  Array<OneD, NekDouble> tmpc(ncoeffs);
569  for (int j = 0; j < m_shapedim; ++j)
570  {
571  Vmath::Vmul(nq, &fluxvector[0][0], 1, &m_DivMF[j][0], 1, &tmp[0], 1);
572 
573  Vmath::Neg(nq, &tmp[0], 1);
574  m_fields[0]->IProductWRTBase(tmp, tmpc);
575 
576  Vmath::Vadd(ncoeffs, outarray[j + 1], 1, tmpc, 1, outarray[j + 1], 1);
577  }
578 }
Nektar::SolverUtils::MMFSystem::m_DivMF
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Definition: MMFSystem.h:189
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::MMFSWE::GetSWEFluxVector
void GetSWEFluxVector(const int i, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
Definition: MMFSWE.cpp:580
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ AddElevationEffect()

void Nektar::MMFSWE::AddElevationEffect ( Array< OneD, Array< OneD, NekDouble >> &  physarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray 
)
protected

Definition at line 1244 of file MMFSWE.cpp.

References m_depth, m_Derivdepth, Nektar::SolverUtils::EquationSystem::m_fields, m_g, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vmul().

Referenced by DoOdeRhs().

1246 {
1247  int ncoeffs = outarray[0].num_elements();
1248  int nq = physarray[0].num_elements();
1249 
1250  Array<OneD, NekDouble> h(nq);
1251  Array<OneD, NekDouble> tmp(nq);
1252  Array<OneD, NekDouble> tmpc(ncoeffs);
1253 
1254  // physarray is primitive
1255  // conservative formulation compute h
1256  // h = \eta + d
1257  Vmath::Vadd(nq, physarray[0], 1, m_depth, 1, h, 1);
1258 
1259  for (int j = 0; j < m_shapedim; ++j)
1260  {
1261  Vmath::Vmul(nq, &h[0], 1, &m_Derivdepth[j][0], 1, &tmp[0], 1);
1262  Vmath::Smul(nq, m_g, tmp, 1, tmp, 1);
1263 
1264  m_fields[0]->IProductWRTBase(tmp, tmpc);
1265 
1266  Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1267  }
1268 }
Nektar::MMFSWE::m_Derivdepth
Array< OneD, Array< OneD, NekDouble > > m_Derivdepth
Definition: MMFSWE.h:87
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ AddRotation()

void Nektar::MMFSWE::AddRotation ( Array< OneD, Array< OneD, NekDouble >> &  physarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray 
)
protected

Definition at line 1273 of file MMFSWE.cpp.

References Compute_demdt_cdot_ek(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Neg(), Vmath::Vadd(), Vmath::Vmul(), and Vmath::Vvtvp().

Referenced by DoOdeRhs().

1275 {
1276  // routine works for both primitive and conservative formulations
1277  int ncoeffs = outarray[0].num_elements();
1278  int nq = physarray[0].num_elements();
1279 
1280  // Compute h
1281  Array<OneD, NekDouble> h(nq);
1282  Vmath::Vadd(nq, &physarray[0][0], 1, &m_depth[0], 1, &h[0], 1);
1283 
1284  Array<OneD, NekDouble> de0dt_cdot_e0;
1285  Array<OneD, NekDouble> de0dt_cdot_e1;
1286  Array<OneD, NekDouble> de1dt_cdot_e0;
1287  Array<OneD, NekDouble> de1dt_cdot_e1;
1288  Compute_demdt_cdot_ek(0, 0, physarray, de0dt_cdot_e0);
1289  Compute_demdt_cdot_ek(1, 0, physarray, de1dt_cdot_e0);
1290  Compute_demdt_cdot_ek(0, 1, physarray, de0dt_cdot_e1);
1291  Compute_demdt_cdot_ek(1, 1, physarray, de1dt_cdot_e1);
1292 
1293  Array<OneD, NekDouble> Rott1(nq);
1294  Array<OneD, NekDouble> Rott2(nq);
1295  Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e0, 1, Rott1, 1);
1296  Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e1, 1, Rott2, 1);
1297  Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e0, 1, Rott1, 1, Rott1, 1);
1298  Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e1, 1, Rott2, 1, Rott2, 1);
1299 
1300  // Multiply H and \partial \phi / \partial t which is assumed to be u_{\phi}
1301  Vmath::Vmul(nq, &h[0], 1, &Rott1[0], 1, &Rott1[0], 1);
1302  Vmath::Vmul(nq, &h[0], 1, &Rott2[0], 1, &Rott2[0], 1);
1303 
1304  Vmath::Neg(nq, Rott1, 1);
1305  Vmath::Neg(nq, Rott2, 1);
1306 
1307  Array<OneD, NekDouble> tmpc1(ncoeffs);
1308  Array<OneD, NekDouble> tmpc2(ncoeffs);
1309  m_fields[0]->IProductWRTBase(Rott1, tmpc1);
1310  m_fields[0]->IProductWRTBase(Rott2, tmpc2);
1311 
1312  Vmath::Vadd(ncoeffs, tmpc1, 1, outarray[1], 1, outarray[1], 1);
1313  Vmath::Vadd(ncoeffs, tmpc2, 1, outarray[2], 1, outarray[2], 1);
1314 }
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::MMFSWE::Compute_demdt_cdot_ek
void Compute_demdt_cdot_ek(const int indm, const int indk, const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1320
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ AverageFlux()

void Nektar::MMFSWE::AverageFlux ( const int  index,
NekDouble  hL,
NekDouble  uL,
NekDouble  vL,
NekDouble  hR,
NekDouble  uR,
NekDouble  vR,
Array< OneD, NekDouble > &  numfluxF,
Array< OneD, NekDouble > &  numfluxB 
)
protected

Definition at line 958 of file MMFSWE.cpp.

References ComputeMagAndDot(), and m_g.

Referenced by NumericalSWEFlux().

962 {
963  NekDouble MageF1, MageF2, MageB1, MageB2;
964  NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
965  NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
966 
967  NekDouble g = m_g;
968  NekDouble uRF, vRF, uLB, vLB;
969 
970  ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
971  eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
972 
973  // uRF = uR component in moving frames e^{Fwd}
974  // vRF = vR component in moving frames e^{Fwd}
975  uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
976  vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
977 
978  numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
979  numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
980  numfluxF[2] = 0.0;
981 
982  numfluxF[3] =
983  0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
984  numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
985  numfluxF[5] = 0.0;
986 
987  numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
988  numfluxF[7] =
989  0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
990  numfluxF[8] = 0.0;
991 
992  // uLB = uL component in moving frames e^{Bwd}
993  // vLB = vL component in moving frames e^{Bwd}
994  uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
995  vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
996 
997  numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
998  numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
999  numfluxB[2] = 0.0;
1000 
1001  numfluxB[3] =
1002  0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1003  numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1004  numfluxB[5] = 0.0;
1005 
1006  numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1007  numfluxB[7] =
1008  0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1009  numfluxB[8] = 0.0;
1010 }
Nektar::MMFSWE::ComputeMagAndDot
void ComputeMagAndDot(const int index, NekDouble &MageF1, NekDouble &MageF2, NekDouble &MageB1, NekDouble &MageB2, NekDouble &eF1_cdot_eB1, NekDouble &eF1_cdot_eB2, NekDouble &eF2_cdot_eB1, NekDouble &eF2_cdot_eB2)
Definition: MMFSWE.cpp:1161
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94

◆ Checkpoint_Output_Cartesian()

void Nektar::MMFSWE::Checkpoint_Output_Cartesian ( std::string  outname)
private

Definition at line 2767 of file MMFSWE.cpp.

References ComputeVorticity(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, m_H0, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Neg(), Vmath::Sadd(), Vmath::Smul(), Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vvtvp(), and Nektar::SolverUtils::EquationSystem::WriteFld().

Referenced by v_SetInitialConditions().

2768 {
2769  int nq = m_fields[0]->GetTotPoints();
2770  int ncoeffs = m_fields[0]->GetNcoeffs();
2771 
2772  NekDouble rad_earth = 6.37122 * 1000000;
2773 
2774  int nvariables = 7;
2775 
2776  // vector u in Cartesian coordinates
2777  std::vector<std::string> variables(nvariables);
2778 
2779  variables[0] = "eta";
2780  variables[1] = "hstar";
2781  variables[2] = "vorticity";
2782  variables[3] = "ux";
2783  variables[4] = "uy";
2784  variables[5] = "uz";
2785  variables[6] = "null";
2786 
2787  // Obtain \vec{u} in cartesian coordinate
2788  Array<OneD, Array<OneD, NekDouble>> fieldphys(nvariables);
2789  std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvariables);
2790  for (int i = 0; i < nvariables; ++i)
2791  {
2792  fieldphys[i] = Array<OneD, NekDouble>(nq, 0.0);
2793  fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2794  }
2795 
2796  Vmath::Smul(nq, rad_earth, &(m_fields[0]->GetPhys())[0], 1,
2797  &fieldphys[0][0], 1);
2798  // Vmath::Vadd(nq, &(m_fields[0]->GetPhys())[0], 1, &m_depth[0], 1,
2799  // &fieldphys[1][0], 1);
2800 
2801  Vmath::Vcopy(nq, &m_depth[0], 1, &fieldphys[1][0], 1);
2802  Vmath::Neg(nq, &fieldphys[1][0], 1);
2803  Vmath::Sadd(nq, m_H0, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2804  Vmath::Smul(nq, rad_earth, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2805 
2806  Array<OneD, NekDouble> utmp(nq);
2807  Array<OneD, NekDouble> vtmp(nq);
2808 
2809  Vmath::Vcopy(nq, &(m_fields[1]->GetPhys())[0], 1, &utmp[0], 1);
2810  Vmath::Vcopy(nq, &(m_fields[2]->GetPhys())[0], 1, &vtmp[0], 1);
2811 
2812  ComputeVorticity(utmp, vtmp, fieldphys[2]);
2813  // u_x = u e^1_x + v e^2_x
2814  int indx = 3;
2815  for (int k = 0; k < m_spacedim; ++k)
2816  {
2817  Vmath::Vmul(nq, &utmp[0], 1, &m_movingframes[0][k * nq], 1,
2818  &fieldphys[k + indx][0], 1);
2819  Vmath::Vvtvp(nq, &vtmp[0], 1, &m_movingframes[1][k * nq], 1,
2820  &fieldphys[k + indx][0], 1, &fieldphys[k + indx][0], 1);
2821  }
2822 
2823  for (int i = 0; i < nvariables; ++i)
2824  {
2825  m_fields[0]->FwdTrans(fieldphys[i], fieldcoeffs[i]);
2826  }
2827 
2828  WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2829 }
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::MMFSWE::m_H0
NekDouble m_H0
Definition: MMFSWE.h:95
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:180
Nektar::MMFSWE::ComputeVorticity
void ComputeVorticity(const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
Definition: MMFSWE.cpp:2210
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:318
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1124
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ Compute_demdt_cdot_ek()

void Nektar::MMFSWE::Compute_demdt_cdot_ek ( const int  indm,
const int  indk,
const Array< OneD, const Array< OneD, NekDouble >> &  physarray,
Array< OneD, NekDouble > &  outarray 
)
protected

Definition at line 1320 of file MMFSWE.cpp.

References Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_shapedim, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vcopy(), Vmath::Vmul(), and Vmath::Vvtvp().

Referenced by AddRotation().

1324 {
1325  int j, k;
1326  int nq = m_fields[0]->GetNpoints();
1327 
1328  Array<OneD, NekDouble> tmp(nq, 0.0);
1329  Array<OneD, NekDouble> tmpderiv(nq);
1330 
1331  outarray = Array<OneD, NekDouble>(nq, 0.0);
1332  for (j = 0; j < m_shapedim; ++j)
1333  {
1334  for (k = 0; k < m_spacedim; ++k)
1335  {
1336  // Compute d e^m / d \xi_1 and d e^m / d \xi_2
1337  Vmath::Vcopy(nq, &m_movingframes[indm][k * nq], 1, &tmp[0], 1);
1338  m_fields[0]->PhysDirectionalDeriv(m_movingframes[j], tmp, tmpderiv);
1339 
1340  Vmath::Vmul(nq, &physarray[j + 1][0], 1, &tmpderiv[0], 1,
1341  &tmpderiv[0], 1);
1342 
1343  Vmath::Vvtvp(nq, &tmpderiv[0], 1, &m_movingframes[indk][k * nq], 1,
1344  &outarray[0], 1, &outarray[0], 1);
1345  }
1346  }
1347 }
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:180
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ ComputeEnergy()

NekDouble Nektar::MMFSWE::ComputeEnergy ( const Array< OneD, const NekDouble > &  eta,
const Array< OneD, const NekDouble > &  u,
const Array< OneD, const NekDouble > &  v 
)
protected

Definition at line 2150 of file MMFSWE.cpp.

References m_depth, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, Vmath::Sadd(), Vmath::Smul(), Vmath::Vadd(), Vmath::Vmul(), Vmath::Vsub(), and Vmath::Vvtvp().

Referenced by v_DoSolve(), and v_SetInitialConditions().

2153 {
2154  int nq = m_fields[0]->GetTotPoints();
2155 
2156  Array<OneD, NekDouble> tmp(nq);
2157  Array<OneD, NekDouble> htmp(nq);
2158  Array<OneD, NekDouble> hstmp(nq);
2159 
2160  Vmath::Vmul(nq, u, 1, u, 1, tmp, 1);
2161  Vmath::Vvtvp(nq, v, 1, v, 1, tmp, 1, tmp, 1);
2162  Vmath::Vmul(nq, eta, 1, tmp, 1, tmp, 1);
2163 
2164  Vmath::Sadd(nq, m_H0, eta, 1, htmp, 1);
2165  Vmath::Vmul(nq, htmp, 1, htmp, 1, htmp, 1);
2166 
2167  Vmath::Sadd(nq, -1.0 * m_H0, m_depth, 1, hstmp, 1);
2168  Vmath::Vmul(nq, hstmp, 1, hstmp, 1, hstmp, 1);
2169 
2170  Vmath::Vsub(nq, htmp, 1, hstmp, 1, htmp, 1);
2171  Vmath::Smul(nq, m_g, htmp, 1, htmp, 1);
2172 
2173  Vmath::Vadd(nq, htmp, 1, tmp, 1, tmp, 1);
2174  Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2175 
2176  return m_fields[0]->PhysIntegral(tmp);
2177 }
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::MMFSWE::m_H0
NekDouble m_H0
Definition: MMFSWE.h:95
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:318
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ ComputeEnstrophy()

NekDouble Nektar::MMFSWE::ComputeEnstrophy ( const Array< OneD, const NekDouble > &  eta,
const Array< OneD, const NekDouble > &  u,
const Array< OneD, const NekDouble > &  v 
)
protected

Definition at line 2179 of file MMFSWE.cpp.

References ComputeVorticity(), m_coriolis, m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Smul(), Vmath::Vadd(), Vmath::Vdiv(), and Vmath::Vmul().

Referenced by v_DoSolve(), and v_SetInitialConditions().

2182 {
2183  int nq = m_fields[0]->GetTotPoints();
2184 
2185  Array<OneD, NekDouble> hstartmp(nq);
2186  Array<OneD, NekDouble> tmp(nq);
2187 
2188  Vmath::Vadd(nq, eta, 1, m_depth, 1, hstartmp, 1);
2189 
2190  Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2191  Array<OneD, Array<OneD, NekDouble>> Curlu(m_spacedim);
2192  for (int i = 0; i < m_spacedim; ++i)
2193  {
2194  uvec[i] = Array<OneD, NekDouble>(nq);
2195  Curlu[i] = Array<OneD, NekDouble>(nq);
2196  }
2197 
2198  ComputeVorticity(u, v, tmp);
2199 
2200  Vmath::Vadd(nq, m_coriolis, 1, tmp, 1, tmp, 1);
2201  Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
2202  Vmath::Vdiv(nq, tmp, 1, hstartmp, 1, tmp, 1);
2203  Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2204 
2205  return m_fields[0]->PhysIntegral(tmp);
2206 }
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:244
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::MMFSWE::ComputeVorticity
void ComputeVorticity(const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
Definition: MMFSWE.cpp:2210
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Nektar::MMFSWE::m_coriolis
Array< OneD, NekDouble > m_coriolis
Coriolis force.
Definition: MMFSWE.h:90
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ Computehhuhvflux()

void Nektar::MMFSWE::Computehhuhvflux ( NekDouble  hL,
NekDouble  uL,
NekDouble  vL,
NekDouble  hR,
NekDouble  uR,
NekDouble  vR,
NekDouble  hstar,
NekDouble &  hflux,
NekDouble &  huflux,
NekDouble &  hvflux 
)
protected

Definition at line 892 of file MMFSWE.cpp.

References m_g.

Referenced by RiemannSolverHLLC().

896 {
897  NekDouble g = m_g;
898 
899  NekDouble hC, huC, hvC, SL, SR, Sstar;
900  NekDouble cL = sqrt(g * hL);
901  NekDouble cR = sqrt(g * hR);
902 
903  // Compute SL
904  if (hstar > hL)
905  SL = uL - cL * sqrt(0.5 * ((hstar * hstar + hstar * hL) / (hL * hL)));
906  else
907  SL = uL - cL;
908 
909  // Compute SR
910  if (hstar > hR)
911  SR = uR + cR * sqrt(0.5 * ((hstar * hstar + hstar * hR) / (hR * hR)));
912  else
913  SR = uR + cR;
914 
915  if (fabs(hR * (uR - SR) - hL * (uL - SL)) <= 1.0e-15)
916  Sstar = 0.0;
917  else
918  Sstar = (SL * hR * (uR - SR) - SR * hL * (uL - SL)) /
919  (hR * (uR - SR) - hL * (uL - SL));
920 
921  if (SL >= 0)
922  {
923  hflux = hL * uL;
924  huflux = uL * uL * hL + 0.5 * g * hL * hL;
925  hvflux = hL * uL * vL;
926  }
927  else if (SR <= 0)
928  {
929  hflux = hR * uR;
930  huflux = uR * uR * hR + 0.5 * g * hR * hR;
931  hvflux = hR * uR * vR;
932  }
933  else
934  {
935  if ((SL < 0) && (Sstar >= 0))
936  {
937  hC = hL * ((SL - uL) / (SL - Sstar));
938  huC = hC * Sstar;
939  hvC = hC * vL;
940 
941  hflux = hL * uL + SL * (hC - hL);
942  huflux = (uL * uL * hL + 0.5 * g * hL * hL) + SL * (huC - hL * uL);
943  hvflux = (uL * vL * hL) + SL * (hvC - hL * vL);
944  }
945  else
946  {
947  hC = hR * ((SR - uR) / (SR - Sstar));
948  huC = hC * Sstar;
949  hvC = hC * vR;
950 
951  hflux = hR * uR + SR * (hC - hR);
952  huflux = (uR * uR * hR + 0.5 * g * hR * hR) + SR * (huC - hR * uR);
953  hvflux = (uR * vR * hR) + SR * (hvC - hR * vR);
954  }
955  }
956 }
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94

◆ ComputeMagAndDot()

void Nektar::MMFSWE::ComputeMagAndDot ( const int  index,
NekDouble &  MageF1,
NekDouble &  MageF2,
NekDouble &  MageB1,
NekDouble &  MageB2,
NekDouble &  eF1_cdot_eB1,
NekDouble &  eF1_cdot_eB2,
NekDouble &  eF2_cdot_eB1,
NekDouble &  eF2_cdot_eB2 
)
protected

Definition at line 1161 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::m_MFtraceBwd, and Nektar::SolverUtils::MMFSystem::m_MFtraceFwd.

Referenced by AverageFlux(), LaxFriedrichFlux(), and RusanovFlux().

1166 {
1167  NekDouble MF1x, MF1y, MF1z, MF2x, MF2y, MF2z;
1168  NekDouble MB1x, MB1y, MB1z, MB2x, MB2y, MB2z;
1169 
1170  MF1x = m_MFtraceFwd[0][0][index];
1171  MF1y = m_MFtraceFwd[0][1][index];
1172  MF1z = m_MFtraceFwd[0][2][index];
1173 
1174  MF2x = m_MFtraceFwd[1][0][index];
1175  MF2y = m_MFtraceFwd[1][1][index];
1176  MF2z = m_MFtraceFwd[1][2][index];
1177 
1178  MB1x = m_MFtraceBwd[0][0][index];
1179  MB1y = m_MFtraceBwd[0][1][index];
1180  MB1z = m_MFtraceBwd[0][2][index];
1181 
1182  MB2x = m_MFtraceBwd[1][0][index];
1183  MB2y = m_MFtraceBwd[1][1][index];
1184  MB2z = m_MFtraceBwd[1][2][index];
1185 
1186  // MFtrace = MFtrace [ j*spacedim + k ], j = shape, k = sapce
1187  MageF1 = MF1x * MF1x + MF1y * MF1y + MF1z * MF1z;
1188  MageF2 = MF2x * MF2x + MF2y * MF2y + MF2z * MF2z;
1189  MageB1 = MB1x * MB1x + MB1y * MB1y + MB1z * MB1z;
1190  MageB2 = MB2x * MB2x + MB2y * MB2y + MB2z * MB2z;
1191 
1192  eF1_cdot_eB1 = MF1x * MB1x + MF1y * MB1y + MF1z * MB1z;
1193  eF1_cdot_eB2 = MF1x * MB2x + MF1y * MB2y + MF1z * MB2z;
1194  eF2_cdot_eB1 = MF2x * MB1x + MF2y * MB1y + MF2z * MB1z;
1195  eF2_cdot_eB2 = MF2x * MB2x + MF2y * MB2y + MF2z * MB2z;
1196 }
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::MMFSystem::m_MFtraceBwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
Definition: MMFSystem.h:194
Nektar::SolverUtils::MMFSystem::m_MFtraceFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
Definition: MMFSystem.h:193

◆ ComputeMass()

NekDouble Nektar::MMFSWE::ComputeMass ( const Array< OneD, const NekDouble > &  eta)
protected

Definition at line 2140 of file MMFSWE.cpp.

References m_depth, Nektar::SolverUtils::EquationSystem::m_fields, and Vmath::Vadd().

Referenced by v_DoSolve(), and v_SetInitialConditions().

2141 {
2142  int nq = m_fields[0]->GetTotPoints();
2143 
2144  Array<OneD, NekDouble> tmp(nq);
2145  Vmath::Vadd(nq, eta, 1, m_depth, 1, tmp, 1);
2146 
2147  return m_fields[0]->PhysIntegral(tmp);
2148 }
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302

◆ ComputeNablaCdotVelocity()

void Nektar::MMFSWE::ComputeNablaCdotVelocity ( Array< OneD, NekDouble > &  vellc)
protected

Definition at line 2231 of file MMFSWE.cpp.

References Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_shapedim, Nektar::SolverUtils::EquationSystem::m_spacedim, m_velocity, Vmath::Vvtvp(), and Vmath::Zero().

2232 {
2233  int nq = m_fields[0]->GetNpoints();
2234 
2235  Array<OneD, NekDouble> velcoeff(nq, 0.0);
2236 
2237  Array<OneD, NekDouble> Dtmp0(nq);
2238  Array<OneD, NekDouble> Dtmp1(nq);
2239  Array<OneD, NekDouble> Dtmp2(nq);
2240  Array<OneD, NekDouble> Drv(nq);
2241 
2242  vellc = Array<OneD, NekDouble>(nq, 0.0);
2243 
2244  // m_vellc = \nabla m_vel \cdot tan_i
2245  Array<OneD, NekDouble> tmp(nq);
2246  Array<OneD, NekDouble> vessel(nq);
2247 
2248  for (int j = 0; j < m_shapedim; ++j)
2249  {
2250  Vmath::Zero(nq, velcoeff, 1);
2251  for (int k = 0; k < m_spacedim; ++k)
2252  {
2253  // a_j = tan_j cdot m_vel
2254  Vmath::Vvtvp(nq, &m_movingframes[j][k * nq], 1, &m_velocity[k][0],
2255  1, &velcoeff[0], 1, &velcoeff[0], 1);
2256  }
2257 
2258  // d a_j / d x^k
2259  m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
2260 
2261  for (int k = 0; k < m_spacedim; ++k)
2262  {
2263  // tan_j^k ( d a_j / d x^k )
2264  switch (k)
2265  {
2266  case (0):
2267  {
2268  Vmath::Vvtvp(nq, &Dtmp0[0], 1, &m_movingframes[j][k * nq],
2269  1, &vellc[0], 1, &vellc[0], 1);
2270  }
2271  break;
2272 
2273  case (1):
2274  {
2275  Vmath::Vvtvp(nq, &Dtmp1[0], 1, &m_movingframes[j][k * nq],
2276  1, &vellc[0], 1, &vellc[0], 1);
2277  }
2278  break;
2279 
2280  case (2):
2281  {
2282  Vmath::Vvtvp(nq, &Dtmp2[0], 1, &m_movingframes[j][k * nq],
2283  1, &vellc[0], 1, &vellc[0], 1);
2284  }
2285  break;
2286  }
2287  }
2288  }
2289 }
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:180
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Nektar::MMFSWE::m_velocity
Array< OneD, Array< OneD, NekDouble > > m_velocity
Advection velocity.
Definition: MMFSWE.h:107
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376

◆ ComputeUnstableJetEta()

NekDouble Nektar::MMFSWE::ComputeUnstableJetEta ( const NekDouble  theta)
protected

Definition at line 2737 of file MMFSWE.cpp.

References ComputeUnstableJetuphi(), and m_Omega.

Referenced by UnstableJetFlow().

2738 {
2739  NekDouble uphi, f, dh;
2740 
2741  uphi = ComputeUnstableJetuphi(theta);
2742  f = 2.0 * m_Omega * sin(theta);
2743 
2744  dh = f * uphi + tan(theta) * uphi * uphi;
2745 
2746  return dh;
2747 }
Nektar::MMFSWE::m_Omega
NekDouble m_Omega
Definition: MMFSWE.h:95
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::ComputeUnstableJetuphi
NekDouble ComputeUnstableJetuphi(const NekDouble theta)
Definition: MMFSWE.cpp:2749

◆ ComputeUnstableJetuphi()

NekDouble Nektar::MMFSWE::ComputeUnstableJetuphi ( const NekDouble  theta)
protected

Definition at line 2749 of file MMFSWE.cpp.

References m_en, m_theta0, m_theta1, and m_uthetamax.

Referenced by ComputeUnstableJetEta(), and UnstableJetFlow().

2750 {
2751  NekDouble uphi;
2752 
2753  if ((theta > m_theta0) && (theta < m_theta1))
2754  {
2755  uphi = (m_uthetamax / m_en) *
2756  exp(1.0 / (theta - m_theta0) / (theta - m_theta1));
2757  }
2758 
2759  else
2760  {
2761  uphi = 0.0;
2762  }
2763 
2764  return uphi;
2765 }
Nektar::MMFSWE::m_uthetamax
NekDouble m_uthetamax
Definition: MMFSWE.h:98
Nektar::MMFSWE::m_en
NekDouble m_en
Definition: MMFSWE.h:98
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_theta0
NekDouble m_theta0
Definition: MMFSWE.h:98
Nektar::MMFSWE::m_theta1
NekDouble m_theta1
Definition: MMFSWE.h:98

◆ ComputeVorticity()

void Nektar::MMFSWE::ComputeVorticity ( const Array< OneD, const NekDouble > &  u,
const Array< OneD, const NekDouble > &  v,
Array< OneD, NekDouble > &  Vorticity 
)
protected

Definition at line 2210 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::m_CurlMF, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Vmath::Neg(), Vmath::Vadd(), and Vmath::Vvtvp().

Referenced by Checkpoint_Output_Cartesian(), ComputeEnstrophy(), TestVorticityComputation(), and v_DoSolve().

2213 {
2214  int nq = m_fields[0]->GetTotPoints();
2215 
2216  Array<OneD, NekDouble> tmp(nq);
2217 
2218  Vorticity = Array<OneD, NekDouble>(nq, 0.0);
2219 
2220  m_fields[0]->PhysDirectionalDeriv(m_movingframes[0], v, Vorticity);
2221  Vmath::Vvtvp(nq, &v[0], 1, &m_CurlMF[1][2][0], 1, &Vorticity[0], 1,
2222  &Vorticity[0], 1);
2223 
2224  m_fields[0]->PhysDirectionalDeriv(m_movingframes[1], u, tmp);
2225  Vmath::Neg(nq, tmp, 1);
2226  Vmath::Vvtvp(nq, &u[0], 1, &m_CurlMF[0][2][0], 1, &tmp[0], 1, &tmp[0], 1);
2227 
2228  Vmath::Vadd(nq, tmp, 1, Vorticity, 1, Vorticity, 1);
2229 }
Nektar::SolverUtils::MMFSystem::m_CurlMF
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Definition: MMFSystem.h:190
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:180
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302

◆ ConservativeToPrimitive() [1/2]

void Nektar::MMFSWE::ConservativeToPrimitive ( const Array< OneD, const Array< OneD, NekDouble >> &  physin,
Array< OneD, Array< OneD, NekDouble >> &  physout 
)
protected

Definition at line 2832 of file MMFSWE.cpp.

References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vsub().

2835 {
2836  int nq = GetTotPoints();
2837 
2838  if (physin[0].get() == physout[0].get())
2839  {
2840  // copy indata and work with tmp array
2841  Array<OneD, Array<OneD, NekDouble>> tmp(3);
2842  for (int i = 0; i < 3; ++i)
2843  {
2844  // deep copy
2845  tmp[i] = Array<OneD, NekDouble>(nq);
2846  Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2847  }
2848 
2849  // \eta = h - d
2850  Vmath::Vsub(nq, tmp[0], 1, m_depth, 1, physout[0], 1);
2851 
2852  // u = hu/h
2853  Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
2854 
2855  // v = hv/h
2856  Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
2857  }
2858  else
2859  {
2860  // \eta = h - d
2861  Vmath::Vsub(nq, physin[0], 1, m_depth, 1, physout[0], 1);
2862 
2863  // u = hu/h
2864  Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
2865 
2866  // v = hv/h
2867  Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
2868  }
2869 }
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:244
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064

◆ ConservativeToPrimitive() [2/2]

void Nektar::MMFSWE::ConservativeToPrimitive ( void  )
protected

Definition at line 2912 of file MMFSWE.cpp.

References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Vdiv(), and Vmath::Vsub().

Referenced by DoOdeProjection(), DoOdeRhs(), and v_DoSolve().

2913 {
2914  int nq = GetTotPoints();
2915 
2916  // u = hu/h
2917  Vmath::Vdiv(nq, m_fields[1]->GetPhys(), 1, m_fields[0]->GetPhys(), 1,
2918  m_fields[1]->UpdatePhys(), 1);
2919 
2920  // v = hv/ v
2921  Vmath::Vdiv(nq, m_fields[2]->GetPhys(), 1, m_fields[0]->GetPhys(), 1,
2922  m_fields[2]->UpdatePhys(), 1);
2923 
2924  // \eta = h - d
2925  Vmath::Vsub(nq, m_fields[0]->GetPhys(), 1, m_depth, 1,
2926  m_fields[0]->UpdatePhys(), 1);
2927 }
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:244
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339

◆ create()

static SolverUtils::EquationSystemSharedPtr Nektar::MMFSWE::create ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const SpatialDomains::MeshGraphSharedPtr &  pGraph 
)
inlinestatic

Creates an instance of this class.

Definition at line 67 of file MMFSWE.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

70  {
71  SolverUtils::EquationSystemSharedPtr p =
72  MemoryManager<MMFSWE>::AllocateSharedPtr(pSession, pGraph);
73  p->InitObject();
74  return p;
75  }
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::SolverUtils::EquationSystemSharedPtr
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition: EquationSystem.h:59

◆ DoOdeProjection()

void Nektar::MMFSWE::DoOdeProjection ( const Array< OneD, const Array< OneD, NekDouble >> &  inarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray,
const NekDouble  time 
)
protected

Compute the projection.

Compute the projection for the linear advection equation.

Parameters
inarrayGiven fields.
outarrayCalculated solution.
timeTime.

Definition at line 1356 of file MMFSWE.cpp.

References ASSERTL0, ConservativeToPrimitive(), Nektar::MultiRegions::eDiscontinuous, Nektar::SolverUtils::EquationSystem::m_projectionType, PrimitiveToConservative(), and SetBoundaryConditions().

Referenced by v_InitObject().

1359 {
1360  switch (m_projectionType)
1361  {
1362  case MultiRegions::eDiscontinuous:
1363  {
1364  ConservativeToPrimitive(inarray, outarray);
1365  SetBoundaryConditions(outarray, time);
1366  PrimitiveToConservative(outarray, outarray);
1367  }
1368  break;
1369  default:
1370  ASSERTL0(false, "Unknown projection scheme");
1371  break;
1372  }
1373 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::EquationSystem::m_projectionType
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:387
Nektar::MMFSWE::ConservativeToPrimitive
void ConservativeToPrimitive()
Definition: MMFSWE.cpp:2912
Nektar::MMFSWE::PrimitiveToConservative
void PrimitiveToConservative()
Definition: MMFSWE.cpp:2929
Nektar::MMFSWE::SetBoundaryConditions
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble >> &inarray, NekDouble time)
Definition: MMFSWE.cpp:1375

◆ DoOdeRhs()

void Nektar::MMFSWE::DoOdeRhs ( const Array< OneD, const Array< OneD, NekDouble >> &  inarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray,
const NekDouble  time 
)
protected

Compute the RHS.

Definition at line 424 of file MMFSWE.cpp.

References AddCoriolis(), AddDivForGradient(), AddElevationEffect(), AddRotation(), ConservativeToPrimitive(), Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_AddCoriolis, m_AddRotation, Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Neg(), and WeakDGSWEDirDeriv().

Referenced by v_InitObject().

427 {
428  boost::ignore_unused(time);
429 
430  int i;
431  int nvariables = inarray.num_elements();
432  int ncoeffs = GetNcoeffs();
433  int nq = GetTotPoints();
434 
435  // inarray in physical space
436  Array<OneD, Array<OneD, NekDouble>> physarray(nvariables);
437  Array<OneD, Array<OneD, NekDouble>> modarray(nvariables);
438 
439  for (i = 0; i < nvariables; ++i)
440  {
441  physarray[i] = Array<OneD, NekDouble>(nq);
442  modarray[i] = Array<OneD, NekDouble>(ncoeffs);
443  }
444 
445  // (h, hu, hv) -> (\eta, u, v)
446  ConservativeToPrimitive(inarray, physarray);
447 
448  // Weak Directional Derivative
449  WeakDGSWEDirDeriv(physarray, modarray);
450  AddDivForGradient(physarray, modarray);
451 
452  for (i = 0; i < nvariables; ++i)
453  {
454  Vmath::Neg(ncoeffs, modarray[i], 1);
455  }
456 
457  // coriolis forcing
458  if (m_AddCoriolis)
459  {
460  AddCoriolis(physarray, modarray);
461  }
462 
463  // Bottom Elevation Effect
464  // Add g H \nabla m_depth
465  AddElevationEffect(physarray, modarray);
466 
467  // Add terms concerning the rotation of the moving frame
468  if (m_AddRotation)
469  {
470  AddRotation(physarray, modarray);
471  }
472 
473  for (i = 0; i < nvariables; ++i)
474  {
475  m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
476  m_fields[i]->BwdTrans(modarray[i], outarray[i]);
477  }
478 }
Nektar::MMFSWE::AddDivForGradient
void AddDivForGradient(Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: MMFSWE.cpp:550
Nektar::MMFSWE::ConservativeToPrimitive
void ConservativeToPrimitive()
Definition: MMFSWE.cpp:2912
Nektar::MMFSWE::AddElevationEffect
void AddElevationEffect(Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: MMFSWE.cpp:1244
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::MMFSWE::AddCoriolis
void AddCoriolis(Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: MMFSWE.cpp:1199
Nektar::MMFSWE::m_AddCoriolis
int m_AddCoriolis
Definition: MMFSWE.h:92
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::MMFSWE::WeakDGSWEDirDeriv
void WeakDGSWEDirDeriv(const Array< OneD, Array< OneD, NekDouble >> &InField, Array< OneD, Array< OneD, NekDouble >> &OutField)
Definition: MMFSWE.cpp:480
Nektar::MMFSWE::AddRotation
void AddRotation(Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: MMFSWE.cpp:1273
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::MMFSWE::m_AddRotation
int m_AddRotation
Definition: MMFSWE.h:92
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:646

◆ EvaluateCoriolis()

void Nektar::MMFSWE::EvaluateCoriolis ( void  )
protected

Definition at line 1696 of file MMFSWE.cpp.

References Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, EvaluateCoriolisForZonalFlow(), EvaluateStandardCoriolis(), Nektar::SolverUtils::EquationSystem::GetFunction(), m_coriolis, and m_TestType.

Referenced by v_DoInitialise().

1697 {
1698  switch (m_TestType)
1699  {
1700  case eTestPlane:
1701  {
1702  GetFunction("Coriolis")->Evaluate("f", m_coriolis);
1703  }
1704  break;
1705 
1706  case eTestSteadyZonal:
1707  {
1708  EvaluateCoriolisForZonalFlow(m_coriolis);
1709  }
1710  break;
1711 
1712  case eTestUnsteadyZonal:
1713  case eTestIsolatedMountain:
1714  case eTestUnstableJet:
1715  case eTestRossbyWave:
1716  {
1717  EvaluateStandardCoriolis(m_coriolis);
1718  }
1719  break;
1720 
1721  default:
1722  break;
1723  }
1724 }
Nektar::MMFSWE::m_TestType
TestType m_TestType
Definition: MMFSWE.h:79
Nektar::MMFSWE::m_coriolis
Array< OneD, NekDouble > m_coriolis
Coriolis force.
Definition: MMFSWE.h:90
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:688
Nektar::MMFSWE::EvaluateStandardCoriolis
void EvaluateStandardCoriolis(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1758
Nektar::MMFSWE::EvaluateCoriolisForZonalFlow
void EvaluateCoriolisForZonalFlow(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1726

◆ EvaluateCoriolisForZonalFlow()

void Nektar::MMFSWE::EvaluateCoriolisForZonalFlow ( Array< OneD, NekDouble > &  outarray)
protected

Definition at line 1726 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, and m_Omega.

Referenced by EvaluateCoriolis().

1727 {
1728  int nq = GetTotPoints();
1729  NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1730 
1731  Array<OneD, NekDouble> x(nq);
1732  Array<OneD, NekDouble> y(nq);
1733  Array<OneD, NekDouble> z(nq);
1734 
1735  m_fields[0]->GetCoords(x, y, z);
1736 
1737  NekDouble x0j, x1j, x2j;
1738  NekDouble tmp;
1739 
1740  outarray = Array<OneD, NekDouble>(nq, 0.0);
1741  for (int j = 0; j < nq; ++j)
1742  {
1743  x0j = x[j];
1744  x1j = y[j];
1745  x2j = z[j];
1746 
1747  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1748  cos_theta);
1749 
1750  // H = 2 \Omega *(- \cos \phi \cos \theta \sin \alpha + \sin \theta \cos
1751  // \alpha )
1752  tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
1753  sin_theta * cos(m_alpha);
1754  outarray[j] = 2.0 * m_Omega * tmp;
1755  }
1756 }
Nektar::MMFSWE::m_Omega
NekDouble m_Omega
Definition: MMFSWE.h:95
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_alpha
NekDouble m_alpha
Definition: MMFSWE.h:95
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339

◆ EvaluateStandardCoriolis()

void Nektar::MMFSWE::EvaluateStandardCoriolis ( Array< OneD, NekDouble > &  outarray)
protected

Definition at line 1758 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, and m_Omega.

Referenced by EvaluateCoriolis().

1759 {
1760  int nq = GetTotPoints();
1761 
1762  NekDouble x0j, x1j, x2j;
1763  NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1764 
1765  Array<OneD, NekDouble> x(nq);
1766  Array<OneD, NekDouble> y(nq);
1767  Array<OneD, NekDouble> z(nq);
1768 
1769  m_fields[0]->GetCoords(x, y, z);
1770 
1771  outarray = Array<OneD, NekDouble>(nq, 0.0);
1772  for (int j = 0; j < nq; ++j)
1773  {
1774  x0j = x[j];
1775  x1j = y[j];
1776  x2j = z[j];
1777 
1778  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1779  cos_theta);
1780 
1781  outarray[j] = 2.0 * m_Omega * sin_theta;
1782  }
1783 }
Nektar::MMFSWE::m_Omega
NekDouble m_Omega
Definition: MMFSWE.h:95
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339

◆ EvaluateWaterDepth()

void Nektar::MMFSWE::EvaluateWaterDepth ( void  )
protected

Definition at line 1556 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::AvgAbsInt(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, Vmath::Fill(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, m_Derivdepth, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, m_hs0, m_k2, Nektar::SolverUtils::MMFSystem::m_movingframes, m_Omega, Nektar::SolverUtils::MMFSystem::m_pi, Nektar::SolverUtils::MMFSystem::m_shapedim, m_TestType, Vmath::Vamax(), and Vmath::Zero().

Referenced by v_DoInitialise().

1557 {
1558  int nq = GetTotPoints();
1559 
1560  m_depth = Array<OneD, NekDouble>(nq, 0.0);
1561 
1562  switch (m_TestType)
1563  {
1564  case eTestPlane:
1565  {
1566  Vmath::Fill(nq, 1.0, m_depth, 1);
1567  }
1568  break;
1569 
1570  case eTestUnsteadyZonal:
1571  {
1572  // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1573  Array<OneD, NekDouble> x(nq);
1574  Array<OneD, NekDouble> y(nq);
1575  Array<OneD, NekDouble> z(nq);
1576 
1577  m_fields[0]->GetCoords(x, y, z);
1578 
1579  NekDouble x0j, x1j, x2j;
1580  NekDouble sin_varphi, cos_varphi, sin_theta, cos_theta;
1581 
1582  for (int j = 0; j < nq; ++j)
1583  {
1584  x0j = x[j];
1585  x1j = y[j];
1586  x2j = z[j];
1587 
1588  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1589  sin_theta, cos_theta);
1590  m_depth[j] =
1591  m_H0 - m_k2 -
1592  (0.5 / m_g) * (m_Omega * sin_theta) * (m_Omega * sin_theta);
1593  }
1594  }
1595  break;
1596 
1597  case eTestIsolatedMountain:
1598  {
1599  // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1600  Array<OneD, NekDouble> x(nq);
1601  Array<OneD, NekDouble> y(nq);
1602  Array<OneD, NekDouble> z(nq);
1603 
1604  m_fields[0]->GetCoords(x, y, z);
1605 
1606  int indx = 0;
1607  NekDouble x0j, x1j, x2j, dist2;
1608  NekDouble phi, theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
1609  NekDouble hRad, phic, thetac;
1610  NekDouble Tol = 0.000001;
1611 
1612  hRad = m_pi / 9.0;
1613  phic = -m_pi / 2.0;
1614  thetac = m_pi / 6.0;
1615 
1616  for (int j = 0; j < nq; ++j)
1617  {
1618  x0j = x[j];
1619  x1j = y[j];
1620  x2j = z[j];
1621 
1622  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1623  sin_theta, cos_theta);
1624 
1625  if ((std::abs(sin(phic) - sin_varphi) +
1626  std::abs(sin(thetac) - sin_theta)) < Tol)
1627  {
1628  std::cout << "A point " << j
1629  << " is coincient with the singularity "
1630  << std::endl;
1631  indx = 1;
1632  }
1633 
1634  phi = atan2(sin_varphi, cos_varphi);
1635  theta = atan2(sin_theta, cos_theta);
1636 
1637  // Compute r
1638  dist2 = (phi - phic) * (phi - phic) +
1639  (theta - thetac) * (theta - thetac);
1640 
1641  if (dist2 > hRad * hRad)
1642  {
1643  dist2 = hRad * hRad;
1644  }
1645 
1646  m_depth[j] = m_H0 - m_hs0 * (1.0 - sqrt(dist2) / hRad);
1647  }
1648 
1649  if (!indx)
1650  {
1651  std::cout << "No point is coincident with the singularity point"
1652  << std::endl;
1653  }
1654  }
1655  break;
1656 
1657  case eTestUnstableJet:
1658  {
1659  for (int j = 0; j < nq; ++j)
1660  {
1661  m_depth[j] = m_H0;
1662  }
1663  }
1664  break;
1665 
1666  case eTestSteadyZonal:
1667  case eTestRossbyWave:
1668  {
1669  Vmath::Zero(nq, m_depth, 1);
1670  }
1671  break;
1672 
1673  default:
1674  {
1675  Vmath::Zero(nq, m_depth, 1);
1676  }
1677  break;
1678  }
1679 
1680  // Comptue \nabla m_depth \cdot e^i
1681  m_Derivdepth = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1682 
1683  for (int j = 0; j < m_shapedim; j++)
1684  {
1685  m_Derivdepth[j] = Array<OneD, NekDouble>(nq);
1686  m_fields[0]->PhysDirectionalDeriv(m_movingframes[j], m_depth,
1687  m_Derivdepth[j]);
1688  }
1689 
1690  std::cout << "Water Depth (m_depth) was generated with mag = "
1691  << AvgAbsInt(m_depth) << " with max. deriv = ( "
1692  << Vmath::Vamax(nq, m_Derivdepth[0], 1) << " , "
1693  << Vmath::Vamax(nq, m_Derivdepth[1], 1) << " ) " << std::endl;
1694 }
Nektar::MMFSWE::m_Derivdepth
Array< OneD, Array< OneD, NekDouble > > m_Derivdepth
Definition: MMFSWE.h:87
Nektar::MMFSWE::m_hs0
NekDouble m_hs0
Definition: MMFSWE.h:97
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Nektar::MMFSWE::m_Omega
NekDouble m_Omega
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::MMFSWE::m_TestType
TestType m_TestType
Definition: MMFSWE.h:79
Nektar::MMFSWE::m_H0
NekDouble m_H0
Definition: MMFSWE.h:95
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:180
Nektar::MMFSWE::m_k2
NekDouble m_k2
Definition: MMFSWE.h:96
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Vmath::Vamax
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max. ...
Definition: Vmath.cpp:828
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
Nektar::SolverUtils::MMFSystem::AvgAbsInt
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2336
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376

◆ GetNormalVelocity()

void Nektar::MMFSWE::GetNormalVelocity ( Array< OneD, NekDouble > &  traceVn)
protected

Get the normal velocity.

◆ GetSWEFluxVector()

void Nektar::MMFSWE::GetSWEFluxVector ( const int  i,
const Array< OneD, const Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  flux 
)
protected

Definition at line 580 of file MMFSWE.cpp.

References Blas::Daxpy(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, m_g, Vmath::Vadd(), and Vmath::Vmul().

Referenced by AddDivForGradient(), and WeakDGSWEDirDeriv().

583 {
584  int nq = m_fields[0]->GetTotPoints();
585 
586  switch (i)
587  {
588  // flux function for the h equation = [(\eta + d ) u, (\eta + d) v ]
589  case 0:
590  {
591  // h in flux 1
592  Vmath::Vadd(nq, physfield[0], 1, m_depth, 1, flux[1], 1);
593 
594  // hu in flux 0
595  Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
596 
597  // hv in flux 1
598  Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
599  }
600  break;
601 
602  // flux function for the hu equation = [Hu^2 + 0.5 g H^2, Huv]
603  case 1:
604  {
605  Array<OneD, NekDouble> tmp(nq);
606 
607  // h in tmp
608  Vmath::Vadd(nq, physfield[0], 1, m_depth, 1, tmp, 1);
609 
610  // hu in flux 1
611  Vmath::Vmul(nq, tmp, 1, physfield[1], 1, flux[1], 1);
612 
613  // huu in flux 0
614  Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
615 
616  // hh in tmp
617  Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
618 
619  // huu + 0.5 g hh in flux 0
620  // Daxpy overwrites flux[0] on exit
621  Blas::Daxpy(nq, 0.5 * m_g, tmp, 1, flux[0], 1);
622 
623  // huv in flux 1
624  Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
625  }
626  break;
627 
628  // flux function for the hv equation = [Huv, Hv^2]
629  case 2:
630  {
631  Array<OneD, NekDouble> tmp(nq);
632 
633  // h in tmp
634  Vmath::Vadd(nq, physfield[0], 1, m_depth, 1, tmp, 1);
635 
636  // hv in flux 0
637  Vmath::Vmul(nq, tmp, 1, physfield[2], 1, flux[0], 1);
638 
639  // hvv in flux 1
640  Vmath::Vmul(nq, flux[0], 1, physfield[2], 1, flux[1], 1);
641 
642  // huv in flux 0
643  Vmath::Vmul(nq, flux[0], 1, physfield[1], 1, flux[0], 1);
644 
645  // hh in tmp
646  Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
647 
648  // hvv + 0.5 g hh in flux 1
649  Blas::Daxpy(nq, 0.5 * m_g, tmp, 1, flux[1], 1);
650  }
651  break;
652 
653  // flux function 0.5 g h * h
654  case 3:
655  {
656  Array<OneD, NekDouble> h(nq);
657  flux[0] = Array<OneD, NekDouble>(nq, 0.0);
658 
659  // h in tmp
660  Vmath::Vadd(nq, physfield[0], 1, m_depth, 1, h, 1);
661 
662  // hh in tmp
663  Vmath::Vmul(nq, h, 1, h, 1, h, 1);
664 
665  // 0.5 g hh in flux 0
666  Blas::Daxpy(nq, 0.5 * m_g, h, 1, flux[0], 1);
667  }
668  break;
669 
670  default:
671  break;
672  }
673 }
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Blas::Daxpy
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:110
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ IsolatedMountainFlow()

void Nektar::MMFSWE::IsolatedMountainFlow ( unsigned int  field,
const NekDouble  time,
Array< OneD, NekDouble > &  outfield 
)
protected

Definition at line 2392 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_Omega, Nektar::SolverUtils::EquationSystem::m_spacedim, and m_u0.

Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().

2394 {
2395  boost::ignore_unused(time);
2396 
2397  int nq = GetTotPoints();
2398 
2399  NekDouble uhat, vhat;
2400  NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2401  NekDouble x0j, x1j, x2j;
2402 
2403  Array<OneD, NekDouble> eta(nq, 0.0);
2404  Array<OneD, NekDouble> u(nq, 0.0);
2405  Array<OneD, NekDouble> v(nq, 0.0);
2406 
2407  Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2408  for (int i = 0; i < m_spacedim; ++i)
2409  {
2410  uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2411  }
2412 
2413  Array<OneD, NekDouble> x(nq);
2414  Array<OneD, NekDouble> y(nq);
2415  Array<OneD, NekDouble> z(nq);
2416 
2417  m_fields[0]->GetCoords(x, y, z);
2418 
2419  for (int j = 0; j < nq; ++j)
2420  {
2421  x0j = x[j];
2422  x1j = y[j];
2423  x2j = z[j];
2424 
2425  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2426  cos_theta);
2427 
2428  // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2429  eta[j] = (-1.0 / m_g) * (m_Omega * m_u0 + 0.5 * m_u0 * m_u0) *
2430  sin_theta * sin_theta;
2431 
2432  uhat = m_u0 * cos_theta;
2433  vhat = 0.0;
2434 
2435  uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2436  uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2437  uvec[2][j] = vhat * cos_theta;
2438  }
2439 
2440  // Projection of u onto the tangent plane with conserving the mag. of the
2441  // velocity.
2442  Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2443 
2444  for (int i = 0; i < m_spacedim; ++i)
2445  {
2446  uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2447  }
2448 
2449  // u is projected on the tangent plane with preserving its length
2450  // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2451 
2452  // Change it to the coordinate of moving frames
2453  // CartesianToMovingframes(0,uvecproj,u);
2454  // CartesianToMovingframes(1,uvecproj,v);
2455 
2456  u = CartesianToMovingframes(uvec, 0);
2457  v = CartesianToMovingframes(uvec, 1);
2458 
2459  switch (field)
2460  {
2461  case (0):
2462  {
2463  outfield = eta;
2464  }
2465  break;
2466 
2467  case (1):
2468  {
2469  outfield = u;
2470  }
2471  break;
2472 
2473  case (2):
2474  {
2475  outfield = v;
2476  }
2477  break;
2478  }
2479 }
Nektar::MMFSWE::m_u0
NekDouble m_u0
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_Omega
NekDouble m_Omega
Definition: MMFSWE.h:95
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::SolverUtils::MMFSystem::CartesianToMovingframes
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
Definition: MMFSystem.cpp:774
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339

◆ LaxFriedrichFlux()

void Nektar::MMFSWE::LaxFriedrichFlux ( const int  index,
NekDouble  hL,
NekDouble  uL,
NekDouble  vL,
NekDouble  hR,
NekDouble  uR,
NekDouble  vR,
Array< OneD, NekDouble > &  numfluxF,
Array< OneD, NekDouble > &  numfluxB 
)
protected

Definition at line 1012 of file MMFSWE.cpp.

References ComputeMagAndDot(), m_g, Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd, Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd, and Vmath::Vamax().

Referenced by NumericalSWEFlux().

1016 {
1017  int nvariables = 3;
1018  NekDouble MageF1, MageF2, MageB1, MageB2;
1019  NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1020  NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1021 
1022  NekDouble g = m_g;
1023  NekDouble uRF, vRF, uLB, vLB;
1024  NekDouble velL, velR, lambdaF, lambdaB;
1025 
1026  Array<OneD, NekDouble> EigF(nvariables);
1027  Array<OneD, NekDouble> EigB(nvariables);
1028 
1029  // Compute Magnitude and Dot product of moving frames for the index
1030  ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1031  eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1032 
1033  // Get the velocity in the normal to the edge
1034  velL = uL * m_ncdotMFFwd[0][index] + vL * m_ncdotMFFwd[1][index];
1035  velR = -1.0 * (uR * m_ncdotMFBwd[0][index] + vR * m_ncdotMFBwd[1][index]);
1036 
1037  EigF[0] = velL - sqrt(g * hL);
1038  EigF[1] = velL;
1039  EigF[2] = velL + sqrt(g * hL);
1040 
1041  EigB[0] = velR - sqrt(g * hR);
1042  EigB[1] = velR;
1043  EigB[2] = velR + sqrt(g * hR);
1044 
1045  lambdaF = Vmath::Vamax(nvariables, EigF, 1);
1046  lambdaB = Vmath::Vamax(nvariables, EigB, 1);
1047 
1048  // uRF = uR component in moving frames e^{Fwd}
1049  // vRF = vR component in moving frames e^{Fwd}
1050  uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1051  vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1052 
1053  numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1054  numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1055  numfluxF[2] = 0.5 * lambdaF * (hL - hR);
1056 
1057  numfluxF[3] = 0.5 * (hL * uL * uL * MageF1 + hR * uRF * uRF * MageB1 +
1058  0.5 * g * (hL * hL + hR * hR));
1059  numfluxF[4] = 0.5 * (hL * uL * vL * MageF1 + hR * uRF * vRF * MageB1);
1060  numfluxF[5] = 0.5 * lambdaF * (uL * hL - uRF * hR);
1061 
1062  numfluxF[6] = 0.5 * (hL * uL * vL * MageF2 + hR * uRF * vRF * MageB2);
1063  numfluxF[7] = 0.5 * (hL * vL * vL * MageF2 + hR * vRF * vRF * MageB2 +
1064  0.5 * g * (hL * hL + hR * hR));
1065  numfluxF[8] = 0.5 * lambdaF * (vL * hL - vRF * hR);
1066 
1067  // uLB = uL component in moving frames e^{Bwd}
1068  // vLB = vL component in moving frames e^{Bwd}
1069  uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1070  vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1071 
1072  numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1073  numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1074  numfluxB[2] = 0.5 * lambdaB * (hL - hR);
1075 
1076  numfluxB[3] = 0.5 * (hR * uR * uR * MageB1 + hR * uLB * uLB * MageF1 +
1077  0.5 * g * (hR * hR + hL * hL));
1078  numfluxB[4] = 0.5 * (hR * uR * vR * MageB1 + hR * uLB * vLB * MageF1);
1079  numfluxB[5] = 0.5 * lambdaB * (uLB * hL - uR * hR);
1080 
1081  numfluxB[6] = 0.5 * (hR * uR * vR * MageB2 + hR * uLB * vLB * MageF2);
1082  numfluxB[7] = 0.5 * (hR * vR * vR * MageB2 + hR * vLB * vLB * MageF2 +
1083  0.5 * g * (hR * hR + hL * hL));
1084  numfluxB[8] = 0.5 * lambdaB * (vLB * hL - vR * hR);
1085 }
Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:184
Nektar::MMFSWE::ComputeMagAndDot
void ComputeMagAndDot(const int index, NekDouble &MageF1, NekDouble &MageF2, NekDouble &MageB1, NekDouble &MageB2, NekDouble &eF1_cdot_eB1, NekDouble &eF1_cdot_eB2, NekDouble &eF2_cdot_eB1, NekDouble &eF2_cdot_eB2)
Definition: MMFSWE.cpp:1161
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Vmath::Vamax
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max. ...
Definition: Vmath.cpp:828
Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:183

◆ NumericalSWEFlux()

void Nektar::MMFSWE::NumericalSWEFlux ( Array< OneD, Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  numfluxFwd,
Array< OneD, Array< OneD, NekDouble >> &  numfluxBwd 
)
protected

Definition at line 675 of file MMFSWE.cpp.

References ASSERTL0, AverageFlux(), Nektar::SolverUtils::MMFSystem::CopyBoundaryTrace(), Nektar::SolverUtils::eAverage, Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::eHLLC, Nektar::SolverUtils::eLaxFriedrich, Nektar::SolverUtils::eRusanov, Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), LaxFriedrichFlux(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd, Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd, Nektar::SolverUtils::MMFSystem::m_nperpcdotMFBwd, Nektar::SolverUtils::MMFSystem::m_nperpcdotMFFwd, Nektar::SolverUtils::MMFSystem::m_shapedim, Nektar::SolverUtils::MMFSystem::m_upwindType, RiemannSolverHLLC(), and RusanovFlux().

Referenced by WeakDGSWEDirDeriv().

678 {
679 
680  int i, k;
681  int nTraceNumPoints = GetTraceTotPoints();
682  int nvariables = 3; // only the dependent variables
683 
684  // get temporary arrays
685  Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
686  Array<OneD, Array<OneD, NekDouble>> Bwd(nvariables);
687  Array<OneD, NekDouble> DepthFwd(nTraceNumPoints);
688  Array<OneD, NekDouble> DepthBwd(nTraceNumPoints);
689 
690  for (i = 0; i < nvariables; ++i)
691  {
692  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
693  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
694  }
695 
696  // get the physical values at the trace
697  for (i = 0; i < nvariables; ++i)
698  {
699  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
700 
701  // Copy Fwd to Bwd at boundaries
702  CopyBoundaryTrace(Fwd[i], Bwd[i], SolverUtils::eFwdEQBwd);
703  }
704  m_fields[0]->GetFwdBwdTracePhys(m_depth, DepthFwd, DepthBwd);
705  CopyBoundaryTrace(DepthFwd, DepthBwd, SolverUtils::eFwdEQBwd);
706 
707  // note that we are using the same depth - i.e. the depth is assumed
708  // continuous...
709  switch (m_upwindType)
710  {
711  case SolverUtils::eHLLC:
712  {
713  // rotate the values to the normal direction
714  NekDouble tmpX, tmpY;
715 
716  // Fwd[1] = hu^+ = ( h^1^+ (e^1 \cdot n) + h^2^+ ( e^2 \cdot n) )
717  // Fwd[2] = hv^+ = ( h^1^+ (e^1 \cdot n^{\perp}) + h^2^+ ( e^2 \cdot
718  // n^{\perp}) )
719  // Bwd[1] = hu^- = ( h^1^- (e^1 \cdot n) + h^2^- ( e^2 \cdot n) )
720  // Bwd[2] = hv^- = ( h^1^- (e^1 \cdot n^{\perp}) + h^2^- ( e^2 \cdot
721  // n^{\perp}) )
722 
723  for (k = 0; k < nTraceNumPoints; ++k)
724  {
725  tmpX = Fwd[1][k] * m_ncdotMFFwd[0][k] +
726  Fwd[2][k] * m_ncdotMFFwd[1][k];
727  tmpY = Fwd[1][k] * m_nperpcdotMFFwd[0][k] +
728  Fwd[2][k] * m_nperpcdotMFFwd[1][k];
729  Fwd[1][k] = tmpX;
730  Fwd[2][k] = tmpY;
731 
732  tmpX = Bwd[1][k] * m_ncdotMFBwd[0][k] +
733  Bwd[2][k] * m_ncdotMFBwd[1][k];
734  tmpY = Bwd[1][k] * m_nperpcdotMFBwd[0][k] +
735  Bwd[2][k] * m_nperpcdotMFBwd[1][k];
736  Bwd[1][k] = tmpX;
737  Bwd[2][k] = tmpY;
738  }
739 
740  // Solve the Riemann problem
741  NekDouble denomFwd, denomBwd;
742  Array<OneD, NekDouble> numfluxF(nvariables);
743  Array<OneD, NekDouble> numfluxB(nvariables);
744 
745  NekDouble eF1n, eF2n, eB1n, eB2n;
746  NekDouble eF1t, eF2t, eB1t, eB2t;
747  for (k = 0; k < nTraceNumPoints; ++k)
748  {
749  RiemannSolverHLLC(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
750  Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
751  Bwd[2][k], numfluxF, numfluxB);
752 
753  // uflux = 1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
754  // \cdot n)(e^1 \cdot n^{\perp} ) ] ( e^2 \cdot n^{\perp} huflux
755  // - e^2 \cdot n hvflux
756  // vflux = -1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
757  // \cdot n)(e^1 \cdot n^{\perp} ) ] ( -e^1 \cdot n^{\perp}
758  // huflux + e^1 \cdot n hvflux
759  eF1n = m_ncdotMFFwd[0][k];
760  eF2n = m_ncdotMFFwd[1][k];
761  eB1n = m_ncdotMFBwd[0][k];
762  eB2n = m_ncdotMFBwd[1][k];
763 
764  eF1t = m_nperpcdotMFFwd[0][k];
765  eF2t = m_nperpcdotMFFwd[1][k];
766  eB1t = m_nperpcdotMFBwd[0][k];
767  eB2t = m_nperpcdotMFBwd[1][k];
768 
769  denomFwd = eF1n * eF2t - eF2n * eF1t;
770  denomBwd = eB1n * eB2t - eB2n * eB1t;
771 
772  numfluxFwd[0][k] = numfluxF[0];
773  numfluxFwd[1][k] = (1.0 / denomFwd) *
774  (eF2t * numfluxF[1] - eF2n * numfluxF[2]);
775  numfluxFwd[2][k] =
776  (1.0 / denomFwd) *
777  (-1.0 * eF1t * numfluxF[1] + eF1n * numfluxF[2]);
778 
779  numfluxBwd[0][k] = 1.0 * numfluxB[0];
780  numfluxBwd[1][k] = (1.0 / denomBwd) *
781  (eB2t * numfluxB[1] - eB2n * numfluxB[2]);
782  numfluxBwd[2][k] =
783  (1.0 / denomBwd) *
784  (-1.0 * eB1t * numfluxB[1] + eB1n * numfluxB[2]);
785  }
786  }
787  break;
788 
789  case SolverUtils::eAverage:
790  case SolverUtils::eLaxFriedrich:
791  case SolverUtils::eRusanov:
792  {
793  Array<OneD, NekDouble> numfluxF(nvariables * (m_shapedim + 1));
794  Array<OneD, NekDouble> numfluxB(nvariables * (m_shapedim + 1));
795 
796  for (k = 0; k < nTraceNumPoints; ++k)
797  {
798  if (m_upwindType == SolverUtils::eAverage)
799  {
800  AverageFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
801  Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
802  Bwd[2][k], numfluxF, numfluxB);
803  }
804 
805  else if (m_upwindType == SolverUtils::eLaxFriedrich)
806  {
807  LaxFriedrichFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
808  Fwd[2][k], Bwd[0][k] + DepthFwd[k],
809  Bwd[1][k], Bwd[2][k], numfluxF, numfluxB);
810  }
811 
812  else if (m_upwindType == SolverUtils::eRusanov)
813  {
814  RusanovFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
815  Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
816  Bwd[2][k], numfluxF, numfluxB);
817  }
818 
819  int indx;
820  NekDouble tmpF0, tmpF1, tmpB0, tmpB1;
821  for (i = 0; i < nvariables; ++i)
822  {
823  indx = i * (m_shapedim + 1);
824 
825  tmpF0 = numfluxF[indx] * m_ncdotMFFwd[0][k];
826  tmpF1 = numfluxF[indx + 1] * m_ncdotMFFwd[1][k];
827 
828  tmpB0 = numfluxB[indx] * m_ncdotMFBwd[0][k];
829  tmpB1 = numfluxB[indx + 1] * m_ncdotMFBwd[1][k];
830 
831  numfluxFwd[i][k] = tmpF0 + tmpF1 + numfluxF[indx + 2];
832  numfluxBwd[i][k] = tmpB0 + tmpB1 + numfluxB[indx + 2];
833  }
834  }
835  }
836  break;
837 
838  default:
839  {
840  ASSERTL0(false, "populate switch statement for upwind flux");
841  }
842  break;
843  }
844 }
Nektar::SolverUtils::eAverage
averaged (or centred) flux
Definition: MMFSystem.h:77
Nektar::MMFSWE::AverageFlux
void AverageFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:958
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:184
Nektar::MMFSWE::LaxFriedrichFlux
void LaxFriedrichFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:1012
Nektar::MMFSWE::RusanovFlux
void RusanovFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:1087
Nektar::SolverUtils::eLaxFriedrich
Lax-Friedrich flux.
Definition: MMFSystem.h:78
Nektar::SolverUtils::eHLLC
Harten-Lax-Leer Contact wave flux.
Definition: MMFSystem.h:82
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::SolverUtils::MMFSystem::m_nperpcdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
Definition: MMFSystem.h:186
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:183
Nektar::SolverUtils::MMFSystem::m_nperpcdotMFBwd
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFBwd
Definition: MMFSystem.h:187
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:677
Nektar::MMFSWE::RiemannSolverHLLC
void RiemannSolverHLLC(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:846
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::SolverUtils::MMFSystem::CopyBoundaryTrace
SOLVER_UTILS_EXPORT void CopyBoundaryTrace(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const BoundaryCopyType BDCopyType, const int var=0, const std::string btype="NoUserDefined")
Definition: MMFSystem.cpp:838

◆ PrimitiveToConservative() [1/2]

void Nektar::MMFSWE::PrimitiveToConservative ( const Array< OneD, const Array< OneD, NekDouble >> &  physin,
Array< OneD, Array< OneD, NekDouble >> &  physout 
)
protected

Definition at line 2871 of file MMFSWE.cpp.

References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vmul().

2874 {
2875 
2876  int nq = GetTotPoints();
2877 
2878  if (physin[0].get() == physout[0].get())
2879  {
2880  // copy indata and work with tmp array
2881  Array<OneD, Array<OneD, NekDouble>> tmp(3);
2882  for (int i = 0; i < 3; ++i)
2883  {
2884  // deep copy
2885  tmp[i] = Array<OneD, NekDouble>(nq);
2886  Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2887  }
2888 
2889  // h = \eta + d
2890  Vmath::Vadd(nq, tmp[0], 1, m_depth, 1, physout[0], 1);
2891 
2892  // hu = h * u
2893  Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
2894 
2895  // hv = h * v
2896  Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
2897  }
2898  else
2899  {
2900  // h = \eta + d
2901  Vmath::Vadd(nq, physin[0], 1, m_depth, 1, physout[0], 1);
2902 
2903  // hu = h * u
2904  Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
2905 
2906  // hv = h * v
2907  Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
2908  }
2909 }
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ PrimitiveToConservative() [2/2]

void Nektar::MMFSWE::PrimitiveToConservative ( void  )
protected

Definition at line 2929 of file MMFSWE.cpp.

References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Vadd(), and Vmath::Vmul().

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

2930 {
2931  int nq = GetTotPoints();
2932 
2933  // h = \eta + d
2934  Vmath::Vadd(nq, m_fields[0]->GetPhys(), 1, m_depth, 1,
2935  m_fields[0]->UpdatePhys(), 1);
2936 
2937  // hu = h * u
2938  Vmath::Vmul(nq, m_fields[0]->GetPhys(), 1, m_fields[1]->GetPhys(), 1,
2939  m_fields[1]->UpdatePhys(), 1);
2940 
2941  // hv = h * v
2942  Vmath::Vmul(nq, m_fields[0]->GetPhys(), 1, m_fields[2]->GetPhys(), 1,
2943  m_fields[2]->UpdatePhys(), 1);
2944 }
Nektar::MMFSWE::m_depth
Array< OneD, NekDouble > m_depth
Still water depth.
Definition: MMFSWE.h:86
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ RiemannSolverHLLC()

void Nektar::MMFSWE::RiemannSolverHLLC ( const int  index,
NekDouble  hL,
NekDouble  uL,
NekDouble  vL,
NekDouble  hR,
NekDouble  uR,
NekDouble  vR,
Array< OneD, NekDouble > &  numfluxF,
Array< OneD, NekDouble > &  numfluxB 
)
protected

Definition at line 846 of file MMFSWE.cpp.

References Computehhuhvflux(), and m_g.

Referenced by NumericalSWEFlux().

850 {
851  boost::ignore_unused(index);
852 
853  NekDouble g = m_g;
854 
855  NekDouble cL = sqrt(g * hL);
856  NekDouble cR = sqrt(g * hR);
857 
858  NekDouble uRF, vRF, uLB, vLB;
859  NekDouble hstarF, hstarB;
860 
861  // Temporary assignments
862  uRF = uR;
863  vRF = vR;
864  uLB = uL;
865  vLB = vL;
866 
867  // the two-rarefaction wave assumption
868  hstarF = 0.5 * (cL + cR) + 0.25 * (uL - uRF);
869  hstarF *= hstarF;
870  hstarF *= (1.0 / g);
871 
872  hstarB = 0.5 * (cL + cR) + 0.25 * (uLB - uR);
873  hstarB *= hstarB;
874  hstarB *= (1.0 / g);
875 
876  NekDouble hfluxF, hufluxF, hvfluxF;
877  NekDouble hfluxB, hufluxB, hvfluxB;
878  Computehhuhvflux(hL, uL, vL, hR, uRF, vRF, hstarF, hfluxF, hufluxF,
879  hvfluxF);
880  Computehhuhvflux(hL, uLB, vLB, hR, uR, vR, hstarB, hfluxB, hufluxB,
881  hvfluxB);
882 
883  numfluxF[0] = hfluxF;
884  numfluxF[1] = hufluxF;
885  numfluxF[2] = hvfluxF;
886 
887  numfluxB[0] = hfluxB;
888  numfluxB[1] = hufluxB;
889  numfluxB[2] = hvfluxB;
890 }
Nektar::MMFSWE::Computehhuhvflux
void Computehhuhvflux(NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, NekDouble hstar, NekDouble &hflux, NekDouble &huflux, NekDouble &hvflux)
Definition: MMFSWE.cpp:892
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94

◆ RossbyWave()

void Nektar::MMFSWE::RossbyWave ( unsigned int  field,
Array< OneD, NekDouble > &  outfield 
)
protected

Definition at line 2592 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_angfreq, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, m_K, m_Omega, Nektar::SolverUtils::MMFSystem::m_pi, m_RossbyDisturbance, and Nektar::SolverUtils::EquationSystem::m_spacedim.

Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().

2593 {
2594  int nq = GetTotPoints();
2595  NekDouble uhat, vhat;
2596  NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2597  NekDouble x0j, x1j, x2j;
2598  NekDouble Ath, Bth, Cth, tmp;
2599 
2600  Array<OneD, NekDouble> eta(nq, 0.0);
2601  Array<OneD, NekDouble> u(nq, 0.0);
2602  Array<OneD, NekDouble> v(nq, 0.0);
2603 
2604  Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2605  for (int i = 0; i < m_spacedim; ++i)
2606  {
2607  uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2608  }
2609 
2610  Array<OneD, NekDouble> x(nq);
2611  Array<OneD, NekDouble> y(nq);
2612  Array<OneD, NekDouble> z(nq);
2613 
2614  m_fields[0]->GetCoords(x, y, z);
2615 
2616  NekDouble R = 4.0;
2617  NekDouble cos2theta, cosRtheta, cos6theta, cos2Rtheta, cosRm1theta;
2618  NekDouble cos2phi, cos4phi, sin4phi, cos8phi;
2619 
2620  // disturbancees of Rossby-Haurwitz Wave
2621  NekDouble x0d, y0d, z0d, phi0, theta0;
2622  //NekDouble rad_earth = 6.37122 * 1000000;
2623 
2624  phi0 = 40.0 * m_pi / 180.0;
2625  theta0 = 50.0 * m_pi / 180.0;
2626 
2627  x0d = cos(phi0) * cos(theta0);
2628  y0d = sin(phi0) * cos(theta0);
2629  z0d = sin(theta0);
2630 
2631  for (int j = 0; j < nq; ++j)
2632  {
2633  x0j = x[j];
2634  x1j = y[j];
2635  x2j = z[j];
2636 
2637  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2638  cos_theta);
2639 
2640  // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2641  // \sin \alpha + \sin \theta \cos \alpha )^2
2642 
2643  // tmp = cos^{2R} \theta, R = 4;
2644  cos2theta = cos_theta * cos_theta;
2645  cosRm1theta = cos_theta * cos2theta;
2646  cosRtheta = cos2theta * cos2theta;
2647  cos6theta = cos2theta * cosRtheta;
2648  cos2Rtheta = cosRtheta * cosRtheta;
2649 
2650  tmp = (0.5 * m_angfreq) * (2.0 * m_Omega + m_angfreq);
2651  Ath = tmp * cos2theta +
2652  0.25 * m_K * m_K * cos6theta *
2653  ((R + 1.0) * cosRtheta + (2 * R * R - R - 2.0) * cos2theta -
2654  2 * R * R);
2655 
2656  tmp = (2.0 * m_K) * (m_Omega + m_angfreq) / ((R + 1.0) * (R + 2.0));
2657  Bth = tmp * cosRtheta *
2658  ((R * R + 2 * R + 2) - (R + 1.0) * (R + 1.0) * cos2theta);
2659 
2660  Cth =
2661  0.25 * m_K * m_K * cos2Rtheta * ((R + 1.0) * cos2theta - (R + 2.0));
2662 
2663  // cos (2 \phi) = 2 * cos \phi * cos \phi - 1.0
2664  cos2phi = 2.0 * cos_varphi * cos_varphi - 1.0;
2665  cos4phi = 2.0 * cos2phi * cos2phi - 1.0;
2666  cos8phi = 2.0 * cos4phi * cos4phi - 1.0;
2667 
2668  // sin (2 \phi) = 2 * cos \phi * sin \phi
2669  sin4phi = 4.0 * sin_varphi * cos_varphi * cos2phi;
2670 
2671  eta[j] = m_H0 + (1.0 / m_g) * (Ath + Bth * cos4phi + Cth * cos8phi);
2672 
2673  // disturbances is added
2674  if (m_RossbyDisturbance)
2675  {
2676  eta[j] = eta[j] *
2677  (1.0 + (1.0 / 40.0) * (x0j * x0d + x1j * y0d + x2j * z0d));
2678  }
2679 
2680  // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2681  // || e^2 ||^2
2682  uhat = m_angfreq * cos_theta +
2683  m_K * cosRm1theta * (R * sin_theta * sin_theta - cos2theta) *
2684  cos4phi;
2685  vhat = -1.0 * m_K * R * cosRm1theta * sin_theta * sin4phi;
2686 
2687  uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2688  uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2689  uvec[2][j] = vhat * cos_theta;
2690  }
2691 
2692  // NekDouble etamin, etaaver;
2693  // etamin = Vmath::Vmin(nq, eta, 1)*rad_earth;
2694  // etaaver = Vmath::Vsum(nq, eta, 1)/nq*rad_earth;
2695 
2696  // Projection of u onto the tangent plane with conserving the mag. of the
2697  // velocity.
2698  Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2699 
2700  for (int i = 0; i < m_spacedim; ++i)
2701  {
2702  uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2703  }
2704 
2705  // u is projected on the tangent plane with preserving its length
2706  // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2707 
2708  // Change it to the coordinate of moving frames
2709  // CartesianToMovingframes(0,uvecproj,u);
2710  // CartesianToMovingframes(1,uvecproj,v);
2711 
2712  u = CartesianToMovingframes(uvec, 0);
2713  v = CartesianToMovingframes(uvec, 1);
2714 
2715  switch (field)
2716  {
2717  case (0):
2718  {
2719  outfield = eta;
2720  }
2721  break;
2722 
2723  case (1):
2724  {
2725  outfield = u;
2726  }
2727  break;
2728 
2729  case (2):
2730  {
2731  outfield = v;
2732  }
2733  break;
2734  }
2735 }
Nektar::MMFSWE::m_Omega
NekDouble m_Omega
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_H0
NekDouble m_H0
Definition: MMFSWE.h:95
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::SolverUtils::MMFSystem::CartesianToMovingframes
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
Definition: MMFSystem.cpp:774
Nektar::MMFSWE::m_angfreq
NekDouble m_angfreq
Definition: MMFSWE.h:99
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Nektar::MMFSWE::m_K
NekDouble m_K
Definition: MMFSWE.h:99
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::MMFSWE::m_RossbyDisturbance
int m_RossbyDisturbance
Definition: MMFSWE.h:275

◆ RusanovFlux()

void Nektar::MMFSWE::RusanovFlux ( const int  index,
NekDouble  hL,
NekDouble  uL,
NekDouble  vL,
NekDouble  hR,
NekDouble  uR,
NekDouble  vR,
Array< OneD, NekDouble > &  numfluxF,
Array< OneD, NekDouble > &  numfluxB 
)
protected

Definition at line 1087 of file MMFSWE.cpp.

References ComputeMagAndDot(), m_g, Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd, and Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd.

Referenced by NumericalSWEFlux().

1091 {
1092  int nvariables = 3;
1093  NekDouble MageF1, MageF2, MageB1, MageB2;
1094  NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1095  NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1096 
1097  NekDouble g = m_g;
1098  NekDouble uRF, vRF, uLB, vLB;
1099  NekDouble velL, velR;
1100 
1101  Array<OneD, NekDouble> EigF(nvariables);
1102  Array<OneD, NekDouble> EigB(nvariables);
1103 
1104  // Compute Magnitude and Dot product of moving frames for the index
1105  ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1106  eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1107 
1108  // Get the velocity in the normal to the edge
1109  velL = uL * m_ncdotMFFwd[0][index] + vL * m_ncdotMFFwd[1][index];
1110  velR = -1.0 * (uR * m_ncdotMFBwd[0][index] + vR * m_ncdotMFBwd[1][index]);
1111 
1112  NekDouble SL, SR;
1113  SL = fabs(velL) + sqrt(g * hL);
1114  SR = fabs(velR) + sqrt(g * hR);
1115 
1116  NekDouble S;
1117  if (SL > SR)
1118  S = SL;
1119  else
1120  S = SR;
1121 
1122  // uRF = uR component in moving frames e^{Fwd}
1123  // vRF = vR component in moving frames e^{Fwd}
1124  uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1125  vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1126 
1127  numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1128  numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1129  numfluxF[2] = 0.5 * S * (hL - hR);
1130 
1131  numfluxF[3] =
1132  0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
1133  numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1134  numfluxF[5] = 0.5 * S * (uL * hL - uRF * hR);
1135 
1136  numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1137  numfluxF[7] =
1138  0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1139  numfluxF[8] = 0.5 * S * (vL * hL - vRF * hR);
1140 
1141  // uLB = uL component in moving frames e^{Bwd}
1142  // vLB = vL component in moving frames e^{Bwd}
1143  uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1144  vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1145 
1146  numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1147  numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1148  numfluxB[2] = 0.5 * S * (hL - hR);
1149 
1150  numfluxB[3] =
1151  0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1152  numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1153  numfluxB[5] = 0.5 * S * (uLB * hL - uR * hR);
1154 
1155  numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1156  numfluxB[7] =
1157  0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1158  numfluxB[8] = 0.5 * S * (vLB * hL - vR * hR);
1159 }
Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:184
Nektar::MMFSWE::ComputeMagAndDot
void ComputeMagAndDot(const int index, NekDouble &MageF1, NekDouble &MageF2, NekDouble &MageB1, NekDouble &MageB2, NekDouble &eF1_cdot_eB1, NekDouble &eF1_cdot_eB2, NekDouble &eF2_cdot_eB1, NekDouble &eF2_cdot_eB2)
Definition: MMFSWE.cpp:1161
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:183

◆ SetBoundaryConditions()

void Nektar::MMFSWE::SetBoundaryConditions ( Array< OneD, Array< OneD, NekDouble >> &  inarray,
NekDouble  time 
)
protected

Definition at line 1375 of file MMFSWE.cpp.

References ASSERTL0, Nektar::SolverUtils::EquationSystem::m_expdim, Nektar::SolverUtils::EquationSystem::m_fields, and WallBoundary2D().

Referenced by DoOdeProjection().

1377 {
1378 
1379  int nvariables = m_fields.num_elements();
1380  int cnt = 0;
1381 
1382  // loop over Boundary Regions
1383  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
1384  {
1385 
1386  // Zonal Boundary Condition
1387  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() == "eMG")
1388  {
1389  if (m_expdim == 1)
1390  {
1391  ASSERTL0(false, "Illegal dimension");
1392  }
1393  else if (m_expdim == 2)
1394  {
1395  // ZonalBoundary2D(n,cnt,inarray);
1396  }
1397  }
1398 
1399  // Wall Boundary Condition
1400  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() == "eWall")
1401  {
1402  if (m_expdim == 1)
1403  {
1404  ASSERTL0(false, "Illegal dimension");
1405  }
1406  else if (m_expdim == 2)
1407  {
1408  WallBoundary2D(n, cnt, inarray);
1409  }
1410  }
1411 
1412  // Time Dependent Boundary Condition (specified in meshfile)
1413  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
1414  "eTimeDependent")
1415  {
1416  for (int i = 0; i < nvariables; ++i)
1417  {
1418  m_fields[i]->EvaluateBoundaryConditions(time);
1419  }
1420  }
1421  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
1422  }
1423 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::EquationSystem::m_expdim
int m_expdim
Expansion dimension.
Definition: EquationSystem.h:367
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::MMFSWE::WallBoundary2D
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble >> &physarray)
Definition: MMFSWE.cpp:1425

◆ SteadyZonalFlow()

void Nektar::MMFSWE::SteadyZonalFlow ( unsigned int  field,
Array< OneD, NekDouble > &  outfield 
)
protected

Definition at line 2049 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_H0, m_Hvar, Nektar::SolverUtils::EquationSystem::m_spacedim, and m_u0.

Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().

2051 {
2052  int nq = GetTotPoints();
2053  NekDouble uhat, vhat;
2054  NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2055  NekDouble x0j, x1j, x2j, tmp;
2056 
2057  Array<OneD, NekDouble> eta(nq, 0.0);
2058  Array<OneD, NekDouble> u(nq, 0.0);
2059  Array<OneD, NekDouble> v(nq, 0.0);
2060 
2061  Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2062  for (int i = 0; i < m_spacedim; ++i)
2063  {
2064  uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2065  }
2066 
2067  Array<OneD, NekDouble> x(nq);
2068  Array<OneD, NekDouble> y(nq);
2069  Array<OneD, NekDouble> z(nq);
2070 
2071  m_fields[0]->GetCoords(x, y, z);
2072 
2073  for (int j = 0; j < nq; ++j)
2074  {
2075  x0j = x[j];
2076  x1j = y[j];
2077  x2j = z[j];
2078 
2079  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2080  cos_theta);
2081 
2082  // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2083  // \sin \alpha + \sin \theta \cos \alpha )^2
2084  tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
2085  sin_theta * cos(m_alpha);
2086  eta[j] = m_H0 - m_Hvar * tmp * tmp;
2087 
2088  // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2089  // || e^2 ||^2
2090  uhat = m_u0 * (cos_theta * cos(m_alpha) +
2091  sin_theta * cos_varphi * sin(m_alpha));
2092  vhat = -1.0 * m_u0 * sin_varphi * sin(m_alpha);
2093 
2094  uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2095  uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2096  uvec[2][j] = vhat * cos_theta;
2097  }
2098 
2099  // Projection of u onto the tangent plane with conserving the mag. of the
2100  // velocity.
2101  Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2102 
2103  for (int i = 0; i < m_spacedim; ++i)
2104  {
2105  uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2106  }
2107 
2108  // u is projected on the tangent plane with preserving its length
2109  // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2110 
2111  // Change it to the coordinate of moving frames
2112  // CartesianToMovingframes(0,uvecproj,u);
2113  // CartesianToMovingframes(1,uvecproj,v);
2114 
2115  u = CartesianToMovingframes(uvec, 0);
2116  v = CartesianToMovingframes(uvec, 1);
2117 
2118  switch (field)
2119  {
2120  case (0):
2121  {
2122  outfield = eta;
2123  }
2124  break;
2125 
2126  case (1):
2127  {
2128  outfield = u;
2129  }
2130  break;
2131 
2132  case (2):
2133  {
2134  outfield = v;
2135  }
2136  break;
2137  }
2138 }
Nektar::MMFSWE::m_u0
NekDouble m_u0
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_H0
NekDouble m_H0
Definition: MMFSWE.h:95
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::SolverUtils::MMFSystem::CartesianToMovingframes
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
Definition: MMFSystem.cpp:774
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_alpha
NekDouble m_alpha
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_Hvar
NekDouble m_Hvar
Definition: MMFSWE.h:95
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339

◆ TestSWE2Dproblem()

void Nektar::MMFSWE::TestSWE2Dproblem ( const NekDouble  time,
unsigned int  field,
Array< OneD, NekDouble > &  outfield 
)
private

Definition at line 1984 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::EquationSystem::m_fields, and Nektar::SolverUtils::EquationSystem::m_spacedim.

Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().

1986 {
1987  boost::ignore_unused(time);
1988 
1989  int nq = m_fields[0]->GetNpoints();
1990 
1991  Array<OneD, NekDouble> x0(nq);
1992  Array<OneD, NekDouble> x1(nq);
1993  Array<OneD, NekDouble> x2(nq);
1994 
1995  m_fields[0]->GetCoords(x0, x1, x2);
1996 
1997  Array<OneD, NekDouble> eta0(nq);
1998  Array<OneD, NekDouble> u0(nq);
1999  Array<OneD, NekDouble> v0(nq);
2000 
2001  Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2002 
2003  for (int i = 0; i < m_spacedim; ++i)
2004  {
2005  uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2006  }
2007 
2008  for (int i = 0; i < nq; ++i)
2009  {
2010  eta0[i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2011  (1.0 / cosh(0.395 * x0[i]))) *
2012  (3.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2013  exp(-0.5 * x1[i] * x1[i]);
2014  uvec[0][i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2015  (1.0 / cosh(0.395 * x0[i]))) *
2016  (-9.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2017  exp(-0.5 * x1[i] * x1[i]);
2018  uvec[1][i] = (-2.0 * 0.395 * tanh(0.395 * x0[i])) *
2019  (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2020  (1.0 / cosh(0.395 * x0[i]))) *
2021  (2.0 * x1[i]) * exp(-0.5 * x1[i] * x1[i]);
2022  }
2023 
2024  u0 = CartesianToMovingframes(uvec, 0);
2025  v0 = CartesianToMovingframes(uvec, 1);
2026 
2027  switch (field)
2028  {
2029  case (0):
2030  {
2031  outfield = eta0;
2032  }
2033  break;
2034 
2035  case (1):
2036  {
2037  outfield = u0;
2038  }
2039  break;
2040 
2041  case (2):
2042  {
2043  outfield = v0;
2044  }
2045  break;
2046  }
2047 }
Nektar::SolverUtils::MMFSystem::CartesianToMovingframes
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
Definition: MMFSystem.cpp:774
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339

◆ TestVorticityComputation()

void Nektar::MMFSWE::TestVorticityComputation ( void  )
protected

Definition at line 2946 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::AvgAbsInt(), Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), ComputeVorticity(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vamax(), and Vmath::Vsub().

Referenced by v_InitObject().

2947 {
2948  // Construct beta
2949  int i, k;
2950  int n = 1, m = 1;
2951  int nq = m_fields[0]->GetTotPoints();
2952 
2953  NekDouble alpha, beta_theta, beta_phi;
2954 
2955  NekDouble xp, yp, zp, Re;
2956  NekDouble theta, phi, sin_theta, cos_theta, sin_varphi, cos_varphi;
2957  NekDouble cosntheta3;
2958 
2959  NekDouble thetax, thetay, thetaz;
2960  NekDouble phix, phiy, phiz;
2961 
2962  Array<OneD, NekDouble> x(nq);
2963  Array<OneD, NekDouble> y(nq);
2964  Array<OneD, NekDouble> z(nq);
2965 
2966  Array<OneD, NekDouble> u(nq);
2967  Array<OneD, NekDouble> v(nq);
2968 
2969  Array<OneD, NekDouble> vorticitycompt(nq);
2970  Array<OneD, NekDouble> vorticityexact(nq);
2971 
2972  m_fields[0]->GetCoords(x, y, z);
2973 
2974  Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2975  for (i = 0; i < m_spacedim; ++i)
2976  {
2977  uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2978  }
2979 
2980  // Generate SphericalCoords
2981  for (k = 0; k < nq; ++k)
2982  {
2983  xp = x[k];
2984  yp = y[k];
2985  zp = z[k];
2986 
2987  Re = sqrt(xp * xp + yp * yp + zp * zp);
2988 
2989  CartesianToSpherical(xp, yp, zp, sin_varphi, cos_varphi, sin_theta,
2990  cos_theta);
2991 
2992  alpha = sin_varphi;
2993 
2994  theta = atan2(sin_theta, cos_theta);
2995  phi = atan2(sin_varphi, cos_varphi);
2996 
2997  cosntheta3 = cos(n * theta) * cos(n * theta) * cos(n * theta);
2998 
2999  beta_theta = -4.0 * n * cosntheta3 * cos(m * phi) * sin(n * theta) / Re;
3000  beta_phi = -m * cosntheta3 * sin(m * phi) / Re;
3001 
3002  thetax = -1.0 * cos_varphi * sin_theta;
3003  thetay = -1.0 * sin_varphi * sin_theta;
3004  thetaz = cos_theta;
3005 
3006  phix = -1.0 * sin_varphi;
3007  phiy = cos_varphi;
3008  phiz = 0.0;
3009 
3010  uvec[0][k] = alpha * (beta_theta * thetax + beta_phi * phix);
3011  uvec[1][k] = alpha * (beta_theta * thetay + beta_phi * phiy);
3012  uvec[2][k] = alpha * (beta_theta * thetaz + beta_phi * phiz);
3013 
3014  vorticityexact[k] = -4.0 * n / Re / Re * cos_theta * cos_theta *
3015  cos_varphi * cos(m * phi) * sin(n * theta);
3016  }
3017 
3018  u = CartesianToMovingframes(uvec, 0);
3019  v = CartesianToMovingframes(uvec, 1);
3020 
3021  std::cout << "chi migi1" << std::endl;
3022 
3023  ComputeVorticity(u, v, vorticitycompt);
3024  /*for (int k=0; k < nq; k++)
3025  {
3026 
3027  std::cout << "vorticitycompt[ " << k << "]"<< "\t"<<vorticitycompt[k]<<
3028  std::endl;
3029 
3030  }*/
3031 
3032  Vmath::Vsub(nq, vorticityexact, 1, vorticitycompt, 1, vorticitycompt, 1);
3033 
3034  std::cout << "Vorticity: L2 error = " << AvgAbsInt(vorticitycompt)
3035  << ", Linf error = " << Vmath::Vamax(nq, vorticitycompt, 1)
3036  << std::endl;
3037 }
Nektar::MMFSWE::ComputeVorticity
void ComputeVorticity(const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
Definition: MMFSWE.cpp:2210
Nektar::SolverUtils::MMFSystem::CartesianToMovingframes
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
Definition: MMFSystem.cpp:774
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vamax
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max. ...
Definition: Vmath.cpp:828
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Nektar::SolverUtils::MMFSystem::AvgAbsInt
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2336
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339

◆ UnstableJetFlow()

void Nektar::MMFSWE::UnstableJetFlow ( unsigned int  field,
const NekDouble  time,
Array< OneD, NekDouble > &  outfield 
)
protected

Definition at line 2481 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), ComputeUnstableJetEta(), ComputeUnstableJetuphi(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_hbar, Nektar::SolverUtils::MMFSystem::m_pi, m_PurturbedJet, and Nektar::SolverUtils::EquationSystem::m_spacedim.

Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().

2483 {
2484  boost::ignore_unused(time);
2485 
2486  int nq = GetTotPoints();
2487 
2488  NekDouble uhat, vhat;
2489  NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2490  NekDouble x0j, x1j, x2j;
2491  NekDouble Ttheta, Tphi;
2492 
2493  Array<OneD, NekDouble> eta(nq, 0.0);
2494  Array<OneD, NekDouble> u(nq, 0.0);
2495  Array<OneD, NekDouble> v(nq, 0.0);
2496 
2497  Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2498  for (int i = 0; i < m_spacedim; ++i)
2499  {
2500  uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2501  }
2502 
2503  Array<OneD, NekDouble> x(nq);
2504  Array<OneD, NekDouble> y(nq);
2505  Array<OneD, NekDouble> z(nq);
2506 
2507  m_fields[0]->GetCoords(x, y, z);
2508 
2509  int Nint = 1000;
2510  NekDouble dth, intj;
2511  for (int j = 0; j < nq; ++j)
2512  {
2513  x0j = x[j];
2514  x1j = y[j];
2515  x2j = z[j];
2516 
2517  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2518  cos_theta);
2519 
2520  Ttheta = atan2(sin_theta, cos_theta);
2521  Tphi = atan2(sin_varphi, cos_varphi);
2522 
2523  uhat = ComputeUnstableJetuphi(Ttheta);
2524  vhat = 0.0;
2525 
2526  // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2527  dth = Ttheta / Nint;
2528  eta[j] = dth * 0.5 *
2529  (ComputeUnstableJetEta(0.0) + ComputeUnstableJetEta(Ttheta));
2530  for (int i = 1; i < Nint - 1; i++)
2531  {
2532  intj = i * dth;
2533  eta[j] = eta[j] + dth * ComputeUnstableJetEta(intj);
2534  }
2535  eta[j] = (-1.0 / m_g) * eta[j];
2536 
2537  // Add perturbation
2538  if (m_PurturbedJet)
2539  {
2540  eta[j] =
2541  eta[j] +
2542  m_hbar * cos_theta * exp(-9.0 * Tphi * Tphi) *
2543  exp(-225.0 * (m_pi / 4.0 - Ttheta) * (m_pi / 4.0 - Ttheta));
2544  }
2545 
2546  uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2547  uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2548  uvec[2][j] = vhat * cos_theta;
2549  }
2550 
2551  // Projection of u onto the tangent plane with conserving the mag. of the
2552  // velocity.
2553  Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2554 
2555  for (int i = 0; i < m_spacedim; ++i)
2556  {
2557  uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2558  }
2559 
2560  // u is projected on the tangent plane with preserving its length
2561  // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2562 
2563  // Change it to the coordinate of moving frames
2564  // CartesianToMovingframes(0,uvecproj,u);
2565  // CartesianToMovingframes(1,uvecproj,v);
2566 
2567  u = CartesianToMovingframes(uvec, 0);
2568  v = CartesianToMovingframes(uvec, 1);
2569 
2570  switch (field)
2571  {
2572  case (0):
2573  {
2574  outfield = eta;
2575  }
2576  break;
2577 
2578  case (1):
2579  {
2580  outfield = u;
2581  }
2582  break;
2583 
2584  case (2):
2585  {
2586  outfield = v;
2587  }
2588  break;
2589  }
2590 }
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::SolverUtils::MMFSystem::CartesianToMovingframes
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
Definition: MMFSystem.cpp:774
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Nektar::MMFSWE::m_hbar
NekDouble m_hbar
Definition: MMFSWE.h:98
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_PurturbedJet
int m_PurturbedJet
Definition: MMFSWE.h:276
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::MMFSWE::ComputeUnstableJetuphi
NekDouble ComputeUnstableJetuphi(const NekDouble theta)
Definition: MMFSWE.cpp:2749
Nektar::MMFSWE::ComputeUnstableJetEta
NekDouble ComputeUnstableJetEta(const NekDouble theta)
Definition: MMFSWE.cpp:2737

◆ UnsteadyZonalFlow()

void Nektar::MMFSWE::UnsteadyZonalFlow ( unsigned int  field,
const NekDouble  time,
Array< OneD, NekDouble > &  outfield 
)
protected

Definition at line 2291 of file MMFSWE.cpp.

References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_Omega, Nektar::SolverUtils::EquationSystem::m_spacedim, and m_u0.

Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().

2293 {
2294  int nq = GetTotPoints();
2295  NekDouble uhat, vhat;
2296  NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2297  NekDouble x0j, x1j, x2j, tmp;
2298 
2299  NekDouble TR, Ttheta;
2300 
2301  Array<OneD, NekDouble> eta(nq, 0.0);
2302  Array<OneD, NekDouble> u(nq, 0.0);
2303  Array<OneD, NekDouble> v(nq, 0.0);
2304 
2305  Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2306  for (int i = 0; i < m_spacedim; ++i)
2307  {
2308  uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2309  }
2310 
2311  Array<OneD, NekDouble> x(nq);
2312  Array<OneD, NekDouble> y(nq);
2313  Array<OneD, NekDouble> z(nq);
2314 
2315  m_fields[0]->GetCoords(x, y, z);
2316 
2317  for (int j = 0; j < nq; ++j)
2318  {
2319  x0j = x[j];
2320  x1j = y[j];
2321  x2j = z[j];
2322 
2323  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2324  cos_theta);
2325 
2326  // \eta = ( - ( u_0 ( - T_R sin \alpha \cos \theta + \cos \alpha \sin
2327  // \theta ) + \Omega \sin \theta )^2 + \Omega \sin \theta )^2 + (\Omega
2328  // \sin \theta )^2
2329  TR =
2330  cos_varphi * cos(m_Omega * time) - sin_varphi * sin(m_Omega * time);
2331  Ttheta =
2332  sin_varphi * cos(m_Omega * time) + cos_varphi * sin(m_Omega * time);
2333  tmp = -1.0 * TR * sin(m_alpha) * cos_theta + cos(m_alpha) * sin_theta;
2334 
2335  eta[j] = -1.0 * (m_u0 * tmp + m_Omega * sin_theta) *
2336  (m_u0 * tmp + m_Omega * sin_theta) +
2337  m_Omega * m_Omega * sin_theta * sin_theta;
2338  eta[j] = 0.5 * eta[j] / m_g;
2339 
2340  // u = u_0*(TR*\sin \alpha * \sin \theta + \cos \alpha * \cos \theta
2341  // v = - u_0 Ttheta * \sin \alpha
2342  uhat =
2343  m_u0 * (TR * sin(m_alpha) * sin_theta + cos(m_alpha) * cos_theta);
2344  vhat = -1.0 * m_u0 * Ttheta * sin(m_alpha);
2345 
2346  uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2347  uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2348  uvec[2][j] = vhat * cos_theta;
2349  }
2350 
2351  // Projection of u onto the tangent plane with conserving the mag. of the
2352  // velocity.
2353  Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2354 
2355  for (int i = 0; i < m_spacedim; ++i)
2356  {
2357  uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2358  }
2359 
2360  // u is projected on the tangent plane with preserving its length
2361  // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2362 
2363  // Change it to the coordinate of moving frames
2364  // CartesianToMovingframes(0,uvecproj,u);
2365  // CartesianToMovingframes(1,uvecproj,v);
2366 
2367  u = CartesianToMovingframes(uvec, 0);
2368  v = CartesianToMovingframes(uvec, 1);
2369 
2370  switch (field)
2371  {
2372  case (0):
2373  {
2374  outfield = eta;
2375  }
2376  break;
2377 
2378  case (1):
2379  {
2380  outfield = u;
2381  }
2382  break;
2383 
2384  case (2):
2385  {
2386  outfield = v;
2387  }
2388  break;
2389  }
2390 }
Nektar::MMFSWE::m_u0
NekDouble m_u0
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_Omega
NekDouble m_Omega
Definition: MMFSWE.h:95
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::SolverUtils::MMFSystem::CartesianToMovingframes
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
Definition: MMFSystem.cpp:774
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:365
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Nektar::MMFSWE::m_alpha
NekDouble m_alpha
Definition: MMFSWE.h:95
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339

◆ v_DoInitialise()

void Nektar::MMFSWE::v_DoInitialise ( void  )
protectedvirtual

Sets up initial conditions.

Sets the initial conditions.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1539 of file MMFSWE.cpp.

References EvaluateCoriolis(), EvaluateWaterDepth(), Nektar::SolverUtils::EquationSystem::m_fields, PrimitiveToConservative(), and Nektar::SolverUtils::EquationSystem::SetInitialConditions().

1540 {
1541  // Compute m_depth and m_Derivdepth
1542  EvaluateWaterDepth();
1543  EvaluateCoriolis();
1544  SetInitialConditions();
1545  PrimitiveToConservative();
1546 
1547  // transfer the initial conditions to modal values
1548  for(int i = 0; i < m_fields.num_elements(); ++i)
1549  {
1550  m_fields[i]->SetPhysState(true);
1551  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),m_fields[i]->UpdateCoeffs());
1552  }
1553 
1554 }
Nektar::MMFSWE::EvaluateCoriolis
void EvaluateCoriolis(void)
Definition: MMFSWE.cpp:1696
Nektar::MMFSWE::EvaluateWaterDepth
void EvaluateWaterDepth(void)
Definition: MMFSWE.cpp:1556
Nektar::MMFSWE::PrimitiveToConservative
void PrimitiveToConservative()
Definition: MMFSWE.cpp:2929
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
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:620

◆ v_DoSolve()

void Nektar::MMFSWE::v_DoSolve ( void  )
protectedvirtual

Solves an unsteady problem.

Initialises the time integration scheme (as specified in the session file), and perform the time integration.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 248 of file MMFSWE.cpp.

References ASSERTL0, Nektar::SolverUtils::EquationSystem::Checkpoint_Output(), ComputeEnergy(), ComputeEnstrophy(), ComputeMass(), ComputeVorticity(), ConservativeToPrimitive(), Nektar::NekConstants::kNekZeroTol, Nektar::SolverUtils::UnsteadySystem::m_cflSafetyFactor, Nektar::SolverUtils::EquationSystem::m_checksteps, Nektar::SolverUtils::EquationSystem::m_checktime, m_Energy0, m_Enstrophy0, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_fintime, Nektar::SolverUtils::UnsteadySystem::m_infosteps, Nektar::SolverUtils::UnsteadySystem::m_intScheme, Nektar::SolverUtils::UnsteadySystem::m_intSoln, Nektar::SolverUtils::UnsteadySystem::m_intVariables, m_Mass0, 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_Vorticity0, PrimitiveToConservative(), Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), and Nektar::LibUtilities::Timer::TimePerTest().

249 {
250  ASSERTL0(m_intScheme != 0, "No time integration scheme.");
251 
252  int i, nchk = 1;
253  int nvariables = 0;
254  int nfields = m_fields.num_elements();
255  int nq = m_fields[0]->GetNpoints();
256 
257  if (m_intVariables.empty())
258  {
259  for (i = 0; i < nfields; ++i)
260  {
261  m_intVariables.push_back(i);
262  }
263  nvariables = nfields;
264  }
265  else
266  {
267  nvariables = m_intVariables.size();
268  }
269 
270  // Set up wrapper to fields data storage.
271  Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
272  Array<OneD, Array<OneD, NekDouble>> tmp(nvariables);
273 
274  // Order storage to list time-integrated fields first.
275  for (i = 0; i < nvariables; ++i)
276  {
277  fields[i] = m_fields[m_intVariables[i]]->GetPhys();
278  m_fields[m_intVariables[i]]->SetPhysState(false);
279  }
280 
281  // Initialise time integration scheme
282  m_intSoln =
283  m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
284 
285  // Check uniqueness of checkpoint output
286  ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
287  (m_checktime > 0.0 && m_checksteps == 0) ||
288  (m_checktime == 0.0 && m_checksteps > 0),
289  "Only one of IO_CheckTime and IO_CheckSteps "
290  "should be set!");
291 
292  LibUtilities::Timer timer;
293  bool doCheckTime = false;
294  int step = 0;
295  NekDouble intTime = 0.0;
296  NekDouble cpuTime = 0.0;
297  NekDouble elapsed = 0.0;
298 
299  NekDouble Mass = 0.0, Energy = 0.0, Enstrophy = 0.0, Vorticity = 0.0;
300  Array<OneD, NekDouble> zeta(nq);
301  Array<OneD, Array<OneD, NekDouble>> fieldsprimitive(nvariables);
302  for (int i = 0; i < nvariables; ++i)
303  {
304  fieldsprimitive[i] = Array<OneD, NekDouble>(nq);
305  }
306 
307  while (step < m_steps || m_time < m_fintime - NekConstants::kNekZeroTol)
308  {
309  timer.Start();
310  fields = m_intScheme->TimeIntegrate(step, m_timestep, m_intSoln, m_ode);
311  timer.Stop();
312 
313  m_time += m_timestep;
314  elapsed = timer.TimePerTest(1);
315  intTime += elapsed;
316  cpuTime += elapsed;
317 
318  // Write out status information
319  if (m_session->GetComm()->GetRank() == 0 && !((step + 1) % m_infosteps))
320  {
321  std::cout << "Steps: " << std::setw(8) << std::left << step + 1 << " "
322  << "Time: " << std::setw(12) << std::left << m_time;
323 
324  std::stringstream ss;
325  ss << cpuTime << "s";
326  std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
327  << std::endl;
328 
329  // Printout Mass, Energy, Enstrophy
330  ConservativeToPrimitive(fields, fieldsprimitive);
331 
332  // Vorticity zeta
333  ComputeVorticity(fieldsprimitive[1], fieldsprimitive[2], zeta);
334  Vorticity =
335  std::abs(m_fields[0]->PhysIntegral(zeta) - m_Vorticity0);
336 
337  // Masss = h^*
338  Mass = (ComputeMass(fieldsprimitive[0]) - m_Mass0) / m_Mass0;
339 
340  // Energy = 0.5*( h^*(u^2 + v^2) + g ( h^2 - h_s^s ) )
341  Energy = (ComputeEnergy(fieldsprimitive[0], fieldsprimitive[1],
342  fieldsprimitive[2]) -
343  m_Energy0) /
344  m_Energy0;
345 
346  // Enstrophy = 0.5/h^* ( \mathbf{k} \cdot (\nabla \times \mathbf{v}
347  // ) + f )^2
348  Enstrophy =
349  (ComputeEnstrophy(fieldsprimitive[0], fieldsprimitive[1],
350  fieldsprimitive[2]) -
351  m_Enstrophy0) /
352  m_Enstrophy0;
353 
354  std::cout << "dMass = " << std::setw(8) << std::left << Mass << " "
355  << ", dEnergy = " << std::setw(8) << std::left << Energy
356  << " "
357  << ", dEnstrophy = " << std::setw(8) << std::left
358  << Enstrophy << " "
359  << ", dVorticity = " << std::setw(8) << std::left
360  << Vorticity << std::endl
361  << std::endl;
362 
363  cpuTime = 0.0;
364  }
365 
366  // (h, hu, hv) -> (\eta, u, v)
367  ConservativeToPrimitive();
368 
369  // Transform data into coefficient space
370  for (i = 0; i < nvariables; ++i)
371  {
372  m_fields[m_intVariables[i]]->SetPhys(fields[i]);
373  m_fields[m_intVariables[i]]->FwdTrans_IterPerExp(
374  fields[i], m_fields[m_intVariables[i]]->UpdateCoeffs());
375  m_fields[m_intVariables[i]]->SetPhysState(false);
376  }
377 
378  // Write out checkpoint files
379  if ((m_checksteps && step && !((step + 1) % m_checksteps)) ||
380  doCheckTime)
381  {
382  Checkpoint_Output(nchk++);
383  doCheckTime = false;
384 
385  // (\eta, u, v) -> (h, hu, hv)
386  PrimitiveToConservative();
387  }
388 
389  // Step advance
390  ++step;
391  }
392 
393  // Print out summary statistics
394  if (m_session->GetComm()->GetRank() == 0)
395  {
396  if (m_cflSafetyFactor > 0.0)
397  {
398  std::cout << "CFL safety factor : " << m_cflSafetyFactor << std::endl
399  << "CFL time-step : " << m_timestep << std::endl;
400  }
401 
402  if (m_session->GetSolverInfo("Driver") != "SteadyState")
403  {
404  std::cout << "Time-integration : " << intTime << "s" << std::endl;
405  }
406  }
407 
408  for (i = 0; i < nvariables; ++i)
409  {
410  m_fields[m_intVariables[i]]->SetPhys(fields[i]);
411  m_fields[m_intVariables[i]]->SetPhysState(true);
412  }
413 
414  // (h, hu, hv) -> (\eta, u, v)
415  ConservativeToPrimitive();
416 
417  for (i = 0; i < nvariables; ++i)
418  {
419  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
420  m_fields[i]->UpdateCoeffs());
421  }
422 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::MMFSWE::m_Vorticity0
NekDouble m_Vorticity0
Definition: MMFSWE.h:101
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:347
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:71
Nektar::SolverUtils::EquationSystem::m_timestep
NekDouble m_timestep
Time step size.
Definition: EquationSystem.h:353
Nektar::SolverUtils::EquationSystem::Checkpoint_Output
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition: EquationSystem.cpp:1086
Nektar::MMFSWE::m_Mass0
NekDouble m_Mass0
Definition: MMFSWE.h:101
Nektar::SolverUtils::EquationSystem::m_checktime
NekDouble m_checktime
Time between checkpoints.
Definition: EquationSystem.h:357
Nektar::MMFSWE::ConservativeToPrimitive
void ConservativeToPrimitive()
Definition: MMFSWE.cpp:2912
Nektar::MMFSWE::PrimitiveToConservative
void PrimitiveToConservative()
Definition: MMFSWE.cpp:2929
Nektar::SolverUtils::EquationSystem::m_checksteps
int m_checksteps
Number of steps between checkpoints.
Definition: EquationSystem.h:363
Nektar::SolverUtils::EquationSystem::m_fintime
NekDouble m_fintime
Finish time of the simulation.
Definition: EquationSystem.h:351
Nektar::MMFSWE::ComputeVorticity
void ComputeVorticity(const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
Definition: MMFSWE.cpp:2210
Nektar::MMFSWE::ComputeMass
NekDouble ComputeMass(const Array< OneD, const NekDouble > &eta)
Definition: MMFSWE.cpp:2140
Nektar::MMFSWE::ComputeEnstrophy
NekDouble ComputeEnstrophy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2179
Nektar::SolverUtils::EquationSystem::m_steps
int m_steps
Number of steps to take.
Definition: EquationSystem.h:361
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:50
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_Energy0
NekDouble m_Energy0
Definition: MMFSWE.h:101
Nektar::SolverUtils::UnsteadySystem::m_cflSafetyFactor
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
Definition: UnsteadySystem.h:58
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:333
Nektar::SolverUtils::UnsteadySystem::m_intScheme
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
Wrapper to the time integration scheme.
Definition: UnsteadySystem.h:69
Nektar::SolverUtils::UnsteadySystem::m_infosteps
int m_infosteps
Number of time steps between outputting status information.
Definition: UnsteadySystem.h:62
Nektar::MMFSWE::m_Enstrophy0
NekDouble m_Enstrophy0
Definition: MMFSWE.h:101
Nektar::MMFSWE::ComputeEnergy
NekDouble ComputeEnergy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2150
Nektar::SolverUtils::UnsteadySystem::m_intSoln
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
Definition: UnsteadySystem.h:73

◆ v_EvaluateExactSolution()

void Nektar::MMFSWE::v_EvaluateExactSolution ( unsigned int  field,
Array< OneD, NekDouble > &  outfield,
const NekDouble  time 
)
protectedvirtual

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 3228 of file MMFSWE.cpp.

References Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, IsolatedMountainFlow(), m_TestType, RossbyWave(), SteadyZonalFlow(), TestSWE2Dproblem(), UnstableJetFlow(), and UnsteadyZonalFlow().

Referenced by v_L2Error().

3231 {
3232  switch (m_TestType)
3233  {
3234  case eTestSteadyZonal:
3235  {
3236  SteadyZonalFlow(field, outfield);
3237  }
3238  break;
3239 
3240  case eTestUnsteadyZonal:
3241  {
3242  UnsteadyZonalFlow(field, time, outfield);
3243  }
3244  break;
3245 
3246  case eTestIsolatedMountain:
3247  {
3248  IsolatedMountainFlow(field, time, outfield);
3249  }
3250  break;
3251 
3252  case eTestUnstableJet:
3253  {
3254  UnstableJetFlow(field, time, outfield);
3255  }
3256  break;
3257 
3258  case eTestRossbyWave:
3259  {
3260  RossbyWave(field, outfield);
3261  }
3262  break;
3263 
3264  case eTestPlane:
3265  {
3266  TestSWE2Dproblem(time, field, outfield);
3267  }
3268  break;
3269 
3270  default:
3271  break;
3272  }
3273 }
Nektar::MMFSWE::RossbyWave
void RossbyWave(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2592
Nektar::MMFSWE::IsolatedMountainFlow
void IsolatedMountainFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2392
Nektar::MMFSWE::m_TestType
TestType m_TestType
Definition: MMFSWE.h:79
Nektar::MMFSWE::UnstableJetFlow
void UnstableJetFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2481
Nektar::MMFSWE::SteadyZonalFlow
void SteadyZonalFlow(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2049
Nektar::MMFSWE::TestSWE2Dproblem
void TestSWE2Dproblem(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:1984
Nektar::MMFSWE::UnsteadyZonalFlow
void UnsteadyZonalFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2291

◆ v_GenerateSummary()

void Nektar::MMFSWE::v_GenerateSummary ( SolverUtils::SummaryList &  s)
protectedvirtual

Print Summary.

Reimplemented from Nektar::SolverUtils::MMFSystem.

Definition at line 3275 of file MMFSWE.cpp.

References Nektar::SolverUtils::AddSummaryItem(), Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, m_alpha, m_PurturbedJet, m_RossbyDisturbance, m_TestType, Nektar::TestTypeMap, and Nektar::SolverUtils::MMFSystem::v_GenerateSummary().

3276 {
3277  MMFSystem::v_GenerateSummary(s);
3278  SolverUtils::AddSummaryItem(s, "TestType", TestTypeMap[m_TestType]);
3279 
3280  switch (m_TestType)
3281  {
3282  case eTestSteadyZonal:
3283  {
3284  SolverUtils::AddSummaryItem(s, "Rotation Angle", m_alpha);
3285  }
3286  break;
3287 
3288  case eTestRossbyWave:
3289  {
3290  SolverUtils::AddSummaryItem(s, "RossbyDistrubance",
3291  m_RossbyDisturbance);
3292  }
3293  break;
3294 
3295  case eTestUnstableJet:
3296  {
3297  SolverUtils::AddSummaryItem(s, "PurturbedJet", m_PurturbedJet);
3298  }
3299  break;
3300 
3301  default:
3302  break;
3303  }
3304 }
Nektar::TestTypeMap
const char *const TestTypeMap[]
Definition: MMFDiffusion.h:58
Nektar::MMFSWE::m_TestType
TestType m_TestType
Definition: MMFSWE.h:79
Nektar::SolverUtils::AddSummaryItem
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:49
Nektar::MMFSWE::m_PurturbedJet
int m_PurturbedJet
Definition: MMFSWE.h:276
Nektar::MMFSWE::m_alpha
NekDouble m_alpha
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_RossbyDisturbance
int m_RossbyDisturbance
Definition: MMFSWE.h:275
Nektar::SolverUtils::MMFSystem::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
Definition: MMFSystem.cpp:2492

◆ v_InitObject()

void Nektar::MMFSWE::v_InitObject ( )
protectedvirtual

Initialise the object.

Initialisation object for the unsteady linear advection equation.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 61 of file MMFSWE.cpp.

References ASSERTL0, Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineProjection(), DoOdeProjection(), DoOdeRhs(), Nektar::SolverUtils::eSphere, Nektar::eTestIsolatedMountain, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, m_AddCoriolis, m_AddRotation, m_alpha, m_angfreq, m_en, Nektar::SolverUtils::UnsteadySystem::m_explicitAdvection, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, m_hbar, m_hs0, m_Hvar, m_K, m_k2, Nektar::SolverUtils::UnsteadySystem::m_ode, m_Omega, Nektar::SolverUtils::MMFSystem::m_pi, m_PurturbedJet, m_RossbyDisturbance, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::MMFSystem::m_surfaceType, m_TestType, m_theta0, m_theta1, m_u0, m_uthetamax, Nektar::SolverUtils::MMFSystem::MMFInitObject(), Nektar::SIZE_TestType, Nektar::TestTypeMap, TestVorticityComputation(), and Nektar::SolverUtils::UnsteadySystem::v_InitObject().

62 {
63  // Call to the initialisation object
64  UnsteadySystem::v_InitObject();
65 
66  int nq = m_fields[0]->GetNpoints();
67  int shapedim = m_fields[0]->GetShapeDimension();
68  Array<OneD, Array<OneD, NekDouble>> Anisotropy(shapedim);
69  for (int j = 0; j < shapedim; ++j)
70  {
71  Anisotropy[j] = Array<OneD, NekDouble>(nq, 1.0);
72  }
73 
74  MMFSystem::MMFInitObject(Anisotropy);
75 
76  // Load acceleration of gravity
77  m_session->LoadParameter("Gravity", m_g, 9.81);
78 
79  // Add Coriolois effects
80  m_session->LoadParameter("AddCoriolis", m_AddCoriolis, 1);
81 
82  // Add Rotation of the sphere along the pole
83  m_session->LoadParameter("AddRotation", m_AddRotation, 1);
84 
85  m_session->LoadParameter("AddRossbyDisturbance", m_RossbyDisturbance, 0);
86  m_session->LoadParameter("PurturbedJet", m_PurturbedJet, 0);
87 
88  // Define TestType
89  ASSERTL0(m_session->DefinesSolverInfo("TESTTYPE"),
90  "No TESTTYPE defined in session.");
91  std::string TestTypeStr = m_session->GetSolverInfo("TESTTYPE");
92  for (int i = 0; i < (int)SIZE_TestType; ++i)
93  {
94  if (boost::iequals(TestTypeMap[i], TestTypeStr))
95  {
96  m_TestType = (TestType)i;
97  break;
98  }
99  }
100 
101  // Variable Setting for each test problem
102  NekDouble gms = 9.80616;
103 
104  switch (m_TestType)
105  {
106  case eTestSteadyZonal:
107  {
108  NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
109  NekDouble rad_earth = 6.37122 * 1000000;
110  NekDouble Omegams;
111 
112  // Nondimensionalized coeffs.
113  m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
114 
115  m_session->LoadParameter("RotationAngle", m_alpha, 0.0);
116  m_session->LoadParameter("u0", m_u0, 2.0 * m_pi / 12.0);
117  m_session->LoadParameter("Omega", Omegams, 7.292 * 0.00001);
118  m_Omega = Omegams * SecondToDay;
119 
120  m_session->LoadParameter("H0", m_H0, 2.94 * 10000);
121  m_H0 = m_H0 / (rad_earth * gms);
122 
123  m_Hvar = (1.0 / m_g) * (m_Omega * m_u0 + 0.5 * m_u0 * m_u0);
124  }
125  break;
126 
127  case eTestUnsteadyZonal:
128  {
129  NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
130  NekDouble rad_earth = 6.37122 * 1000000;
131  NekDouble Omegams;
132 
133  // Nondimensionalized coeffs.
134  m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
135 
136  m_session->LoadParameter("RotationAngle", m_alpha, 0.0);
137  m_session->LoadParameter("u0", m_u0, 2.0 * m_pi / 12.0);
138  m_session->LoadParameter("Omega", Omegams, 7.292 * 0.00001);
139  m_Omega = Omegams * SecondToDay;
140 
141  m_H0 = 133681.0 / (rad_earth * gms); // m^2 / s^2
142  m_k2 = 10.0 / (rad_earth * gms); // m^2 / s^2
143  }
144  break;
145 
146  case eTestIsolatedMountain:
147  {
148  NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
149  NekDouble rad_earth = 6.37122 * 1000000;
150  NekDouble Omegams;
151 
152  // Nondimensionalized coeffs.
153  m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
154 
155  m_session->LoadParameter("RotationAngle", m_alpha, 0.0);
156 
157  m_session->LoadParameter("u0", m_u0, 20.0);
158  m_u0 = m_u0 * SecondToDay / rad_earth;
159 
160  m_session->LoadParameter("Omega", Omegams, 7.292 * 0.00001);
161  m_Omega = Omegams * SecondToDay;
162 
163  m_H0 = 5960.0 / rad_earth;
164  m_hs0 = 2000.0 / rad_earth;
165  }
166  break;
167 
168  case eTestUnstableJet:
169  {
170  NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
171  NekDouble rad_earth = 6.37122 * 1000000;
172  NekDouble Omegams;
173 
174  // Nondimensionalized coeffs.
175  m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
176 
177  m_session->LoadParameter("u0", m_u0, 80.0);
178  m_u0 = m_u0 * SecondToDay / rad_earth;
179 
180  m_session->LoadParameter("Omega", Omegams, 7.292 * 0.00001);
181  m_Omega = Omegams * SecondToDay;
182 
183  m_session->LoadParameter("H0", m_H0, 10000.0);
184  m_H0 = m_H0 / rad_earth;
185 
186  m_uthetamax = 80 * SecondToDay / rad_earth;
187  m_theta0 = m_pi / 7.0;
188  m_theta1 = m_pi / 2.0 - m_theta0;
189  m_en = exp(-4.0 / (m_theta1 - m_theta0) / (m_theta1 - m_theta0));
190  m_hbar = 120.0 / rad_earth;
191 
192  std::cout << "m_theta0 = " << m_theta0 << ", m_theta1 = " << m_theta1
193  << ", m_en = " << m_en << ", m_hbar = " << m_hbar
194  << std::endl;
195  }
196  break;
197 
198  case eTestRossbyWave:
199  {
200  NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
201  NekDouble rad_earth = 6.37122 * 1000000;
202  NekDouble Omegams;
203 
204  // Nondimensionalized coeffs.
205  m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
206 
207  m_session->LoadParameter("Omega", Omegams, 7.292 * 0.00001);
208  m_Omega = Omegams * SecondToDay;
209 
210  m_session->LoadParameter("H0", m_H0, 8000.0);
211  m_H0 = m_H0 / rad_earth;
212 
213  m_angfreq = 7.848 * 0.000001 * SecondToDay;
214  m_K = 7.848 * 0.000001 * SecondToDay;
215  }
216  break;
217 
218  default:
219  break;
220  }
221 
222  // TestVorticityComputation
223  if (m_surfaceType == SolverUtils::eSphere)
224  {
225  TestVorticityComputation();
226  }
227 
228  // If explicit it computes RHS and PROJECTION for the time integration
229  if (m_explicitAdvection)
230  {
231  m_ode.DefineOdeRhs(&MMFSWE::DoOdeRhs, this);
232  m_ode.DefineProjection(&MMFSWE::DoOdeProjection, this);
233  }
234  // Otherwise it gives an error (no implicit integration)
235  else
236  {
237  ASSERTL0(false, "Implicit unsteady Advection not set up.");
238  }
239 }
Nektar::MMFSWE::TestVorticityComputation
void TestVorticityComputation()
Definition: MMFSWE.cpp:2946
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::TestTypeMap
const char *const TestTypeMap[]
Definition: MMFDiffusion.h:58
Nektar::MMFSWE::DoOdeProjection
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the projection.
Definition: MMFSWE.cpp:1356
Nektar::MMFSWE::m_uthetamax
NekDouble m_uthetamax
Definition: MMFSWE.h:98
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:71
Nektar::MMFSWE::m_hs0
NekDouble m_hs0
Definition: MMFSWE.h:97
Nektar::MMFSWE::m_u0
NekDouble m_u0
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_Omega
NekDouble m_Omega
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_en
NekDouble m_en
Definition: MMFSWE.h:98
Nektar::MMFSWE::m_TestType
TestType m_TestType
Definition: MMFSWE.h:79
Nektar::MMFSWE::m_H0
NekDouble m_H0
Definition: MMFSWE.h:95
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineProjection
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationScheme.h:215
Nektar::SolverUtils::UnsteadySystem::m_explicitAdvection
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
Definition: UnsteadySystem.h:79
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationScheme.h:191
Nektar::MMFSWE::m_AddCoriolis
int m_AddCoriolis
Definition: MMFSWE.h:92
Nektar::MMFSWE::m_angfreq
NekDouble m_angfreq
Definition: MMFSWE.h:99
Nektar::MMFSWE::m_k2
NekDouble m_k2
Definition: MMFSWE.h:96
Nektar::MMFSWE::m_hbar
NekDouble m_hbar
Definition: MMFSWE.h:98
Nektar::SIZE_TestType
Length of enum list.
Definition: MMFDiffusion.h:55
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::MMFSWE::m_PurturbedJet
int m_PurturbedJet
Definition: MMFSWE.h:276
Nektar::SolverUtils::MMFSystem::MMFInitObject
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:53
Nektar::MMFSWE::m_g
NekDouble m_g
Definition: MMFSWE.h:94
Nektar::MMFSWE::m_K
NekDouble m_K
Definition: MMFSWE.h:99
Nektar::SolverUtils::UnsteadySystem::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:80
Nektar::MMFSWE::m_alpha
NekDouble m_alpha
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_Hvar
NekDouble m_Hvar
Definition: MMFSWE.h:95
Nektar::MMFSWE::m_theta0
NekDouble m_theta0
Definition: MMFSWE.h:98
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:333
Nektar::MMFSWE::m_AddRotation
int m_AddRotation
Definition: MMFSWE.h:92
Nektar::MMFSWE::DoOdeRhs
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the RHS.
Definition: MMFSWE.cpp:424
Nektar::MMFSWE::m_RossbyDisturbance
int m_RossbyDisturbance
Definition: MMFSWE.h:275
Nektar::MMFSWE::m_theta1
NekDouble m_theta1
Definition: MMFSWE.h:98

◆ v_L2Error()

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

Virtual function for the L_2 error computation between fields and a given exact solution.

Compute the error in the L2-norm.

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 3039 of file MMFSWE.cpp.

References Nektar::SolverUtils::EquationSystem::ErrorExtraPoints(), Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::SolverUtils::EquationSystem::m_comm, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError, Nektar::SolverUtils::EquationSystem::m_time, Nektar::LibUtilities::ReduceSum, v_EvaluateExactSolution(), Vmath::Vabs(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

3042 {
3043  boost::ignore_unused(exactsoln);
3044 
3045  int nq = m_fields[field]->GetNpoints();
3046  NekDouble L2error = -1.0;
3047 
3048  if (m_NumQuadPointsError == 0)
3049  {
3050  if (m_fields[field]->GetPhysState() == false)
3051  {
3052  m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
3053  m_fields[field]->UpdatePhys());
3054  }
3055 
3056  switch (field)
3057  {
3058  case (0):
3059  {
3060  // I(h) = I (h - h_exact ) / I (h_exact)
3061  Array<OneD, NekDouble> exactsolution(nq);
3062  v_EvaluateExactSolution(0, exactsolution, m_time);
3063 
3064  // exactsoln = u - u_T so that L2 compute u_T
3065  NekDouble L2exact = m_fields[0]->PhysIntegral(exactsolution);
3066 
3067  Vmath::Vsub(nq, &(m_fields[0]->GetPhys())[0], 1,
3068  &exactsolution[0], 1, &exactsolution[0], 1);
3069  Vmath::Vabs(nq, exactsolution, 1, exactsolution, 1);
3070 
3071  L2error = (m_fields[0]->PhysIntegral(exactsolution)) / L2exact;
3072  }
3073  break;
3074 
3075  case (1):
3076  {
3077  // I2 (u) = I( (u - u_ext)^2 + (v - v_ext)^2 )^{1/2} / I(
3078  // u_ext^2 + v_ext^2 )^{1/2}
3079  Array<OneD, NekDouble> exactu(nq);
3080  Array<OneD, NekDouble> exactv(nq);
3081  Array<OneD, NekDouble> tmp(nq);
3082 
3083  // L2exact = \int (\sqrt{exactu*exactu+exactv*exactv})
3084  NekDouble L2exact;
3085  v_EvaluateExactSolution(1, exactu, m_time);
3086  v_EvaluateExactSolution(2, exactv, m_time);
3087  Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3088  Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3089  Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3090 
3091  L2exact = m_fields[1]->PhysIntegral(tmp);
3092 
3093  // L2exact = \int
3094  // (\sqrt{(u-exactu)*(u-exactu)+(v-exactv)*(v-exactv)})
3095  Vmath::Vsub(nq, &(m_fields[1]->GetPhys())[0], 1, &exactu[0], 1,
3096  &exactu[0], 1);
3097  Vmath::Vsub(nq, &(m_fields[2]->GetPhys())[0], 1, &exactv[0], 1,
3098  &exactv[0], 1);
3099  Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3100  Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3101  Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3102 
3103  L2error = (m_fields[1]->PhysIntegral(tmp)) / L2exact;
3104  }
3105  break;
3106 
3107  case (2):
3108  {
3109  L2error = 0.0;
3110  }
3111  break;
3112 
3113  default:
3114  break;
3115  }
3116 
3117  if (Normalised == true)
3118  {
3119  Array<OneD, NekDouble> one(m_fields[field]->GetNpoints(), 1.0);
3120 
3121  NekDouble Vol = m_fields[field]->PhysIntegral(one);
3122  m_comm->AllReduce(Vol, LibUtilities::ReduceSum);
3123 
3124  L2error = sqrt(L2error * L2error / Vol);
3125  }
3126  }
3127  else
3128  {
3129  Array<OneD, NekDouble> L2INF(2);
3130  L2INF = ErrorExtraPoints(field);
3131  L2error = L2INF[0];
3132  }
3133 
3134  return L2error;
3135 }
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:347
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:411
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Nektar::SolverUtils::EquationSystem::ErrorExtraPoints
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].
Definition: EquationSystem.cpp:857
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:427
Nektar::SolverUtils::EquationSystem::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: EquationSystem.h:331
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:712
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Definition: EquationSystem.h:396
Nektar::MMFSWE::v_EvaluateExactSolution
virtual void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Definition: MMFSWE.cpp:3228
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ v_LinfError()

NekDouble Nektar::MMFSWE::v_LinfError ( unsigned int  field,
const Array< OneD, NekDouble > &  exactsoln 
)
protectedvirtual

Virtual function for the L_inf error computation 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 3137 of file MMFSWE.cpp.

References Nektar::SolverUtils::EquationSystem::EvaluateExactSolution(), Nektar::SolverUtils::EquationSystem::LinfError(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_time, Vmath::Vadd(), Vmath::Vamax(), Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vsub().

3139 {
3140  boost::ignore_unused(exactsoln);
3141 
3142  NekDouble LinfError = -1.0;
3143 
3144  if (m_fields[field]->GetPhysState() == false)
3145  {
3146  m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
3147  m_fields[field]->UpdatePhys());
3148  }
3149 
3150  int nq = m_fields[field]->GetNpoints();
3151 
3152  // Obtain \vec{u} in cartesian coordinate
3153  Array<OneD, NekDouble> x(nq);
3154  Array<OneD, NekDouble> y(nq);
3155  Array<OneD, NekDouble> z(nq);
3156 
3157  m_fields[0]->GetCoords(x, y, z);
3158 
3159  switch (field)
3160  {
3161  case (0):
3162  {
3163  NekDouble Letaint;
3164 
3165  Array<OneD, NekDouble> exactsolution(nq);
3166 
3167  EvaluateExactSolution(field, exactsolution, m_time);
3168  LinfError = m_fields[field]->Linf(m_fields[field]->GetPhys(),
3169  exactsolution);
3170 
3171  Letaint = Vmath::Vamax(nq, exactsolution, 1);
3172 
3173  Vmath::Vsub(nq, &(m_fields[0]->GetPhys())[0], 1, &exactsolution[0],
3174  1, &exactsolution[0], 1);
3175 
3176  LinfError = fabs(LinfError / Letaint);
3177  }
3178  break;
3179 
3180  case (1):
3181  {
3182  Array<OneD, NekDouble> exactu(nq);
3183  Array<OneD, NekDouble> exactv(nq);
3184  Array<OneD, NekDouble> tmpu(nq);
3185  Array<OneD, NekDouble> tmpv(nq);
3186  Array<OneD, NekDouble> Lerr(nq);
3187  Array<OneD, NekDouble> uT(nq);
3188 
3189  EvaluateExactSolution(1, exactu, m_time);
3190  EvaluateExactSolution(2, exactv, m_time);
3191 
3192  // Compute max[sqrt{(u-uex)^2 + (v-vex)^2}]
3193  Vmath::Vcopy(nq, &(m_fields[1]->UpdatePhys())[0], 1, &tmpu[0], 1);
3194  Vmath::Vcopy(nq, &(m_fields[2]->UpdatePhys())[0], 1, &tmpv[0], 1);
3195 
3196  Vmath::Vsub(nq, &exactu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3197  Vmath::Vsub(nq, &exactv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3198 
3199  Vmath::Vmul(nq, &tmpu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3200  Vmath::Vmul(nq, &tmpv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3201 
3202  Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &Lerr[0], 1);
3203  Vmath::Vsqrt(nq, &Lerr[0], 1, &Lerr[0], 1);
3204 
3205  // uT = max[sqrt( u_T^2 + v_T^2 ) ]
3206  Vmath::Vmul(nq, &exactu[0], 1, &exactu[0], 1, &tmpu[0], 1);
3207  Vmath::Vmul(nq, &exactv[0], 1, &exactv[0], 1, &tmpv[0], 1);
3208  Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &uT[0], 1);
3209  Vmath::Vsqrt(nq, &uT[0], 1, &uT[0], 1);
3210 
3211  LinfError = Vmath::Vamax(nq, Lerr, 1) / Vmath::Vamax(nq, uT, 1);
3212  }
3213  break;
3214 
3215  case (2):
3216  {
3217  LinfError = 0.0;
3218  }
3219  break;
3220 
3221  default:
3222  break;
3223  }
3224 
3225  return LinfError;
3226 }
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:347
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:411
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vamax
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max. ...
Definition: Vmath.cpp:828
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::SolverUtils::EquationSystem::LinfError
SOLVER_UTILS_EXPORT NekDouble LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Linf error computation.
Definition: EquationSystem.h:554
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:628
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ v_SetInitialConditions()

void Nektar::MMFSWE::v_SetInitialConditions ( const NekDouble  initialtime,
bool  dumpInitialConditions,
const int  domain 
)
protectedvirtual

Set the physical fields based on a restart file, or a function describing the initial condition given in the session.

Parameters
initialtimeTime at which to evaluate the function.
dumpInitialConditionsWrite the initial condition to file?

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1785 of file MMFSWE.cpp.

References Checkpoint_Output_Cartesian(), ComputeEnergy(), ComputeEnstrophy(), ComputeMass(), Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, Nektar::SolverUtils::EquationSystem::GetTotPoints(), IsolatedMountainFlow(), m_Energy0, m_Enstrophy0, Nektar::SolverUtils::EquationSystem::m_fields, m_Mass0, Nektar::SolverUtils::EquationSystem::m_sessionName, m_TestType, m_Vorticity0, RossbyWave(), SteadyZonalFlow(), TestSWE2Dproblem(), UnstableJetFlow(), UnsteadyZonalFlow(), and Nektar::SolverUtils::EquationSystem::WriteFld().

1788 {
1789  boost::ignore_unused(domain);
1790 
1791  int nq = GetTotPoints();
1792 
1793  switch (m_TestType)
1794  {
1795  case eTestPlane:
1796  {
1797  Array<OneD, NekDouble> eta0(nq);
1798  Array<OneD, NekDouble> u0(nq);
1799  Array<OneD, NekDouble> v0(nq);
1800 
1801  TestSWE2Dproblem(initialtime, 0, eta0);
1802  m_fields[0]->SetPhys(eta0);
1803 
1804  TestSWE2Dproblem(initialtime, 1, u0);
1805  m_fields[1]->SetPhys(u0);
1806 
1807  TestSWE2Dproblem(initialtime, 2, v0);
1808  m_fields[2]->SetPhys(v0);
1809 
1810  // forward transform to fill the modal coeffs
1811  for (int i = 0; i < m_fields.num_elements(); ++i)
1812  {
1813  m_fields[i]->SetPhysState(true);
1814  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1815  m_fields[i]->UpdateCoeffs());
1816  }
1817  }
1818  break;
1819 
1820  case eTestSteadyZonal:
1821  {
1822  Array<OneD, NekDouble> eta0(nq);
1823  Array<OneD, NekDouble> u0(nq);
1824  Array<OneD, NekDouble> v0(nq);
1825  Array<OneD, NekDouble> zeta0(nq);
1826 
1827  SteadyZonalFlow(0, eta0);
1828  m_fields[0]->SetPhys(eta0);
1829 
1830  SteadyZonalFlow(1, u0);
1831  m_fields[1]->SetPhys(u0);
1832 
1833  SteadyZonalFlow(2, v0);
1834  m_fields[2]->SetPhys(v0);
1835 
1836  // ComputeVorticity(u0, v0, zeta0);
1837  m_Vorticity0 = m_fields[0]->PhysIntegral(zeta0);
1838 
1839  m_Mass0 = ComputeMass(eta0);
1840  m_Energy0 = ComputeEnergy(eta0, u0, v0);
1841  m_Enstrophy0 = ComputeEnstrophy(eta0, u0, v0);
1842 
1843  // forward transform to fill the modal coeffs
1844  for (int i = 0; i < m_fields.num_elements(); ++i)
1845  {
1846  m_fields[i]->SetPhysState(true);
1847  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1848  m_fields[i]->UpdateCoeffs());
1849  }
1850  }
1851  break;
1852 
1853  case eTestUnsteadyZonal:
1854  {
1855  Array<OneD, NekDouble> eta0(nq);
1856  Array<OneD, NekDouble> u0(nq);
1857  Array<OneD, NekDouble> v0(nq);
1858 
1859  UnsteadyZonalFlow(0, initialtime, eta0);
1860  m_fields[0]->SetPhys(eta0);
1861 
1862  UnsteadyZonalFlow(1, initialtime, u0);
1863  m_fields[1]->SetPhys(u0);
1864 
1865  UnsteadyZonalFlow(2, initialtime, v0);
1866  m_fields[2]->SetPhys(v0);
1867 
1868  m_Mass0 = ComputeMass(eta0);
1869  m_Energy0 = ComputeEnergy(eta0, u0, v0);
1870  m_Enstrophy0 = ComputeEnstrophy(eta0, u0, v0);
1871 
1872  // forward transform to fill the modal coeffs
1873  for (int i = 0; i < m_fields.num_elements(); ++i)
1874  {
1875  m_fields[i]->SetPhysState(true);
1876  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1877  m_fields[i]->UpdateCoeffs());
1878  }
1879  }
1880  break;
1881 
1882  case eTestIsolatedMountain:
1883  {
1884  Array<OneD, NekDouble> eta0(nq);
1885  Array<OneD, NekDouble> u0(nq);
1886  Array<OneD, NekDouble> v0(nq);
1887 
1888  IsolatedMountainFlow(0, initialtime, eta0);
1889  m_fields[0]->SetPhys(eta0);
1890 
1891  IsolatedMountainFlow(1, initialtime, u0);
1892  m_fields[1]->SetPhys(u0);
1893 
1894  IsolatedMountainFlow(2, initialtime, v0);
1895  m_fields[2]->SetPhys(v0);
1896 
1897  m_Mass0 = ComputeMass(eta0);
1898  m_Energy0 = ComputeEnergy(eta0, u0, v0);
1899  m_Enstrophy0 = ComputeEnstrophy(eta0, u0, v0);
1900 
1901  // forward transform to fill the modal coeffs
1902  for (int i = 0; i < m_fields.num_elements(); ++i)
1903  {
1904  m_fields[i]->SetPhysState(true);
1905  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1906  m_fields[i]->UpdateCoeffs());
1907  }
1908  }
1909  break;
1910 
1911  case eTestUnstableJet:
1912  {
1913  Array<OneD, NekDouble> eta0(nq);
1914  Array<OneD, NekDouble> u0(nq);
1915  Array<OneD, NekDouble> v0(nq);
1916 
1917  UnstableJetFlow(0, initialtime, eta0);
1918  m_fields[0]->SetPhys(eta0);
1919 
1920  UnstableJetFlow(1, initialtime, u0);
1921  m_fields[1]->SetPhys(u0);
1922 
1923  UnstableJetFlow(2, initialtime, v0);
1924  m_fields[2]->SetPhys(v0);
1925 
1926  m_Mass0 = ComputeMass(eta0);
1927  m_Energy0 = ComputeEnergy(eta0, u0, v0);
1928  m_Enstrophy0 = ComputeEnstrophy(eta0, u0, v0);
1929 
1930  // forward transform to fill the modal coeffs
1931  for (int i = 0; i < m_fields.num_elements(); ++i)
1932  {
1933  m_fields[i]->SetPhysState(true);
1934  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1935  m_fields[i]->UpdateCoeffs());
1936  }
1937  }
1938  break;
1939 
1940  case eTestRossbyWave:
1941  {
1942  Array<OneD, NekDouble> eta0(nq);
1943  Array<OneD, NekDouble> u0(nq);
1944  Array<OneD, NekDouble> v0(nq);
1945 
1946  RossbyWave(0, eta0);
1947  m_fields[0]->SetPhys(eta0);
1948 
1949  RossbyWave(1, u0);
1950  m_fields[1]->SetPhys(u0);
1951 
1952  RossbyWave(2, v0);
1953  m_fields[2]->SetPhys(v0);
1954 
1955  m_Mass0 = ComputeMass(eta0);
1956  m_Energy0 = ComputeEnergy(eta0, u0, v0);
1957  m_Enstrophy0 = ComputeEnstrophy(eta0, u0, v0);
1958 
1959  // forward transform to fill the modal coeffs
1960  for (int i = 0; i < m_fields.num_elements(); ++i)
1961  {
1962  m_fields[i]->SetPhysState(true);
1963  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1964  m_fields[i]->UpdateCoeffs());
1965  }
1966  }
1967  break;
1968 
1969  default:
1970  break;
1971  }
1972 
1973  if (dumpInitialConditions)
1974  {
1975  // dump initial conditions to file
1976  std::string outname = m_sessionName + "_initial.chk";
1977  WriteFld(outname);
1978 
1979  outname = m_sessionName + "_initialCART.chk";
1980  Checkpoint_Output_Cartesian(outname);
1981  }
1982 }
Nektar::MMFSWE::m_Vorticity0
NekDouble m_Vorticity0
Definition: MMFSWE.h:101
Nektar::MMFSWE::RossbyWave
void RossbyWave(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2592
Nektar::MMFSWE::m_Mass0
NekDouble m_Mass0
Definition: MMFSWE.h:101
Nektar::MMFSWE::IsolatedMountainFlow
void IsolatedMountainFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2392
Nektar::MMFSWE::m_TestType
TestType m_TestType
Definition: MMFSWE.h:79
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:345
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:702
Nektar::MMFSWE::UnstableJetFlow
void UnstableJetFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2481
Nektar::MMFSWE::ComputeMass
NekDouble ComputeMass(const Array< OneD, const NekDouble > &eta)
Definition: MMFSWE.cpp:2140
Nektar::MMFSWE::ComputeEnstrophy
NekDouble ComputeEnstrophy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2179
Nektar::MMFSWE::Checkpoint_Output_Cartesian
void Checkpoint_Output_Cartesian(std::string outname)
Definition: MMFSWE.cpp:2767
Nektar::MMFSWE::m_Energy0
NekDouble m_Energy0
Definition: MMFSWE.h:101
Nektar::MMFSWE::SteadyZonalFlow
void SteadyZonalFlow(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2049
Nektar::MMFSWE::TestSWE2Dproblem
void TestSWE2Dproblem(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:1984
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1124
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::MMFSWE::m_Enstrophy0
NekDouble m_Enstrophy0
Definition: MMFSWE.h:101
Nektar::MMFSWE::ComputeEnergy
NekDouble ComputeEnergy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2150
Nektar::MMFSWE::UnsteadyZonalFlow
void UnsteadyZonalFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2291

◆ WallBoundary2D()

void Nektar::MMFSWE::WallBoundary2D ( int  bcRegion,
int  cnt,
Array< OneD, Array< OneD, NekDouble >> &  physarray 
)
protected

Definition at line 1425 of file MMFSWE.cpp.

References ASSERTL0, Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_expdim, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd, Nektar::SolverUtils::MMFSystem::m_nperpcdotMFFwd, Vmath::Neg(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vvtvm(), and Vmath::Vvtvp().

Referenced by SetBoundaryConditions().

1427 {
1428 
1429  int i;
1430  int nTraceNumPoints = GetTraceTotPoints();
1431  int nvariables = physarray.num_elements();
1432 
1433  // get physical values of the forward trace
1434  Array<OneD, Array<OneD, NekDouble>> Fwd0(nvariables);
1435  for (i = 0; i < nvariables; ++i)
1436  {
1437  Fwd0[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1438  m_fields[i]->ExtractTracePhys(physarray[i], Fwd0[i]);
1439  }
1440 
1441  Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
1442  for (i = 0; i < nvariables; ++i)
1443  {
1444  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1445  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1446  }
1447 
1448  // Adjust the physical values of the trace to take
1449  // user defined boundaries into account
1450  int e, id1, id2, npts;
1451 
1452  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
1453  ++e)
1454  {
1455  npts = m_fields[0]
1456  ->GetBndCondExpansions()[bcRegion]
1457  ->GetExp(e)
1458  ->GetTotPoints();
1459  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
1460  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
1461  m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(cnt +
1462  e));
1463 
1464  switch (m_expdim)
1465  {
1466  case 1:
1467  {
1468  // negate the forward flux
1469  Vmath::Neg(npts, &Fwd[1][id2], 1);
1470  }
1471  break;
1472  case 2:
1473  {
1474  Array<OneD, NekDouble> tmp_n(npts);
1475  Array<OneD, NekDouble> tmp_t(npts);
1476 
1477  Vmath::Vmul(npts, &Fwd[1][id2], 1, &m_ncdotMFFwd[0][id2], 1,
1478  &tmp_n[0], 1);
1479  Vmath::Vvtvp(npts, &Fwd[2][id2], 1, &m_ncdotMFFwd[1][id2], 1,
1480  &tmp_n[0], 1, &tmp_n[0], 1);
1481 
1482  Vmath::Vmul(npts, &Fwd[1][id2], 1, &m_nperpcdotMFFwd[0][id2], 1,
1483  &tmp_t[0], 1);
1484  Vmath::Vvtvp(npts, &Fwd[2][id2], 1, &m_nperpcdotMFFwd[1][id2],
1485  1, &tmp_t[0], 1, &tmp_t[0], 1);
1486 
1487  // negate the normal flux
1488  Vmath::Neg(npts, tmp_n, 1);
1489 
1490  Array<OneD, NekDouble> denom(npts);
1491  Array<OneD, NekDouble> tmp_u(npts);
1492  Array<OneD, NekDouble> tmp_v(npts);
1493 
1494  // denom = (e^1 \cdot n ) (e^2 \cdot t) - (e^2 \cdot n ) (e^1
1495  // \cdot t)
1496  Vmath::Vmul(npts, &m_ncdotMFFwd[1][id2], 1,
1497  &m_nperpcdotMFFwd[0][id2], 1, &denom[0], 1);
1498  Vmath::Vvtvm(npts, &m_ncdotMFFwd[0][id2], 1,
1499  &m_nperpcdotMFFwd[1][id2], 1, &denom[0], 1,
1500  &denom[0], 1);
1501 
1502  Vmath::Vmul(npts, &m_ncdotMFFwd[1][id2], 1, &tmp_t[0], 1,
1503  &tmp_u[0], 1);
1504  Vmath::Vvtvm(npts, &m_nperpcdotMFFwd[1][id2], 1, &tmp_n[0], 1,
1505  &tmp_u[0], 1, &tmp_u[0], 1);
1506  Vmath::Vdiv(npts, &tmp_u[0], 1, &denom[0], 1, &tmp_u[0], 1);
1507 
1508  Vmath::Vcopy(npts, &tmp_u[0], 1, &Fwd[1][id2], 1);
1509 
1510  Vmath::Vmul(npts, &m_nperpcdotMFFwd[0][id2], 1, &tmp_n[0], 1,
1511  &tmp_v[0], 1);
1512  Vmath::Vvtvm(npts, &m_ncdotMFFwd[0][id2], 1, &tmp_t[0], 1,
1513  &tmp_v[0], 1, &tmp_v[0], 1);
1514  Vmath::Vdiv(npts, &tmp_v[0], 1, &denom[0], 1, &tmp_v[0], 1);
1515 
1516  Vmath::Vcopy(npts, &tmp_v[0], 1, &Fwd[2][id2], 1);
1517  }
1518  break;
1519 
1520  default:
1521  ASSERTL0(false, "Illegal expansion dimension");
1522  }
1523 
1524  // copy boundary adjusted values into the boundary expansion
1525  for (i = 0; i < nvariables; ++i)
1526  {
1527  Vmath::Vcopy(npts, &Fwd[i][id2], 1,
1528  &(m_fields[i]
1529  ->GetBndCondExpansions()[bcRegion]
1530  ->UpdatePhys())[id1],
1531  1);
1532  }
1533  }
1534 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::EquationSystem::m_expdim
int m_expdim
Expansion dimension.
Definition: EquationSystem.h:367
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:244
Nektar::SolverUtils::MMFSystem::m_nperpcdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
Definition: MMFSystem.h:186
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:183
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:677
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Vmath::Vvtvm
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector plus vector): z = w*x - y
Definition: Vmath.cpp:468
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ WeakDGSWEDirDeriv()

void Nektar::MMFSWE::WeakDGSWEDirDeriv ( const Array< OneD, Array< OneD, NekDouble >> &  InField,
Array< OneD, Array< OneD, NekDouble >> &  OutField 
)
protected

Definition at line 480 of file MMFSWE.cpp.

References Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetNpoints(), GetSWEFluxVector(), Nektar::SolverUtils::EquationSystem::GetTraceNpoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Neg(), NumericalSWEFlux(), and Vmath::Vadd().

Referenced by DoOdeRhs().

483 {
484  int i;
485  int nq = GetNpoints();
486  int ncoeffs = GetNcoeffs();
487  int nTracePointsTot = GetTraceNpoints();
488  int nvariables = m_fields.num_elements();
489 
490  Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
491  Array<OneD, Array<OneD, NekDouble>> physfield(nvariables);
492 
493  for (i = 0; i < m_shapedim; ++i)
494  {
495  fluxvector[i] = Array<OneD, NekDouble>(nq);
496  }
497 
498  // InField is Primitive
499  for (i = 0; i < nvariables; ++i)
500  {
501  physfield[i] = InField[i];
502  }
503 
504  // Get the ith component of the flux vector in (physical space)
505  // fluxvector[0] = component for e^1 cdot \nabla \varphi
506  // fluxvector[1] = component for e^2 cdot \nabla \varphi
507  Array<OneD, NekDouble> fluxtmp(nq);
508  Array<OneD, NekDouble> tmp(ncoeffs);
509 
510  // Compute Divergence Components
511  for (i = 0; i < nvariables; ++i)
512  {
513  GetSWEFluxVector(i, physfield, fluxvector);
514 
515  OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
516  for (int j = 0; j < m_shapedim; ++j)
517  {
518  // Directional derivation with respect to the j'th moving frame
519  // tmp_j = ( \nabla \phi, fluxvector[j] \mathbf{e}^j )
520  m_fields[i]->IProductWRTDirectionalDerivBase(m_movingframes[j],
521  fluxvector[j], tmp);
522  Vmath::Vadd(ncoeffs, &tmp[0], 1, &OutField[i][0], 1,
523  &OutField[i][0], 1);
524  }
525  }
526 
527  // Numerical Flux
528  Array<OneD, Array<OneD, NekDouble>> numfluxFwd(nvariables);
529  Array<OneD, Array<OneD, NekDouble>> numfluxBwd(nvariables);
530 
531  for (i = 0; i < nvariables; ++i)
532  {
533  numfluxFwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
534  numfluxBwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
535  }
536 
537  NumericalSWEFlux(physfield, numfluxFwd, numfluxBwd);
538 
539  // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
540  for (i = 0; i < nvariables; ++i)
541  {
542  Vmath::Neg(ncoeffs, OutField[i], 1);
543  m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
544  OutField[i]);
545  m_fields[i]->SetPhysState(false);
546  }
547 }
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:180
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::MMFSWE::GetSWEFluxVector
void GetSWEFluxVector(const int i, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
Definition: MMFSWE.cpp:580
Nektar::MMFSWE::NumericalSWEFlux
void NumericalSWEFlux(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
Definition: MMFSWE.cpp:675
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:712
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:339
Nektar::SolverUtils::EquationSystem::GetTraceNpoints
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:682
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:646
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302

Friends And Related Function Documentation

◆ MemoryManager< MMFSWE >

friend class MemoryManager< MMFSWE >
friend

Definition at line 64 of file MMFSWE.h.

Member Data Documentation

◆ className

std::string Nektar::MMFSWE::className
static
Initial value:
=
SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
"MMFSWE", MMFSWE::create, "MMFSWE equation.")

Name of class.

Definition at line 77 of file MMFSWE.h.

◆ m_AddCoriolis

int Nektar::MMFSWE::m_AddCoriolis
protected

Definition at line 92 of file MMFSWE.h.

Referenced by DoOdeRhs(), and v_InitObject().

◆ m_AddRotation

int Nektar::MMFSWE::m_AddRotation
protected

Definition at line 92 of file MMFSWE.h.

Referenced by DoOdeRhs(), and v_InitObject().

◆ m_alpha

NekDouble Nektar::MMFSWE::m_alpha
protected

Definition at line 95 of file MMFSWE.h.

Referenced by EvaluateCoriolisForZonalFlow(), SteadyZonalFlow(), UnsteadyZonalFlow(), v_GenerateSummary(), and v_InitObject().

◆ m_angfreq

NekDouble Nektar::MMFSWE::m_angfreq
protected

Definition at line 99 of file MMFSWE.h.

Referenced by RossbyWave(), and v_InitObject().

◆ m_coriolis

Array<OneD, NekDouble> Nektar::MMFSWE::m_coriolis
protected

Coriolis force.

Definition at line 90 of file MMFSWE.h.

Referenced by AddCoriolis(), ComputeEnstrophy(), and EvaluateCoriolis().

◆ m_depth

Array<OneD, NekDouble> Nektar::MMFSWE::m_depth
protected

Still water depth.

Definition at line 86 of file MMFSWE.h.

Referenced by AddCoriolis(), AddElevationEffect(), AddRotation(), Checkpoint_Output_Cartesian(), ComputeEnergy(), ComputeEnstrophy(), ComputeMass(), ConservativeToPrimitive(), EvaluateWaterDepth(), GetSWEFluxVector(), NumericalSWEFlux(), and PrimitiveToConservative().

◆ m_Derivdepth

Array<OneD, Array<OneD, NekDouble> > Nektar::MMFSWE::m_Derivdepth
protected

Definition at line 87 of file MMFSWE.h.

Referenced by AddElevationEffect(), and EvaluateWaterDepth().

◆ m_en

NekDouble Nektar::MMFSWE::m_en
protected

Definition at line 98 of file MMFSWE.h.

Referenced by ComputeUnstableJetuphi(), and v_InitObject().

◆ m_Energy0

NekDouble Nektar::MMFSWE::m_Energy0
protected

Definition at line 101 of file MMFSWE.h.

Referenced by v_DoSolve(), and v_SetInitialConditions().

◆ m_Enstrophy0

NekDouble Nektar::MMFSWE::m_Enstrophy0
protected

Definition at line 101 of file MMFSWE.h.

Referenced by v_DoSolve(), and v_SetInitialConditions().

◆ m_g

NekDouble Nektar::MMFSWE::m_g
protected

Definition at line 94 of file MMFSWE.h.

Referenced by AddElevationEffect(), AverageFlux(), ComputeEnergy(), Computehhuhvflux(), EvaluateWaterDepth(), GetSWEFluxVector(), IsolatedMountainFlow(), LaxFriedrichFlux(), RiemannSolverHLLC(), RossbyWave(), RusanovFlux(), UnstableJetFlow(), UnsteadyZonalFlow(), and v_InitObject().

◆ m_H0

NekDouble Nektar::MMFSWE::m_H0
protected

Definition at line 95 of file MMFSWE.h.

Referenced by Checkpoint_Output_Cartesian(), ComputeEnergy(), EvaluateWaterDepth(), RossbyWave(), SteadyZonalFlow(), and v_InitObject().

◆ m_hbar

NekDouble Nektar::MMFSWE::m_hbar
protected

Definition at line 98 of file MMFSWE.h.

Referenced by UnstableJetFlow(), and v_InitObject().

◆ m_hs0

NekDouble Nektar::MMFSWE::m_hs0
protected

Definition at line 97 of file MMFSWE.h.

Referenced by EvaluateWaterDepth(), and v_InitObject().

◆ m_Hvar

NekDouble Nektar::MMFSWE::m_Hvar
protected

Definition at line 95 of file MMFSWE.h.

Referenced by SteadyZonalFlow(), and v_InitObject().

◆ m_K

NekDouble Nektar::MMFSWE::m_K
protected

Definition at line 99 of file MMFSWE.h.

Referenced by RossbyWave(), and v_InitObject().

◆ m_k2

NekDouble Nektar::MMFSWE::m_k2
protected

Definition at line 96 of file MMFSWE.h.

Referenced by EvaluateWaterDepth(), and v_InitObject().

◆ m_Mass0

NekDouble Nektar::MMFSWE::m_Mass0
protected

Definition at line 101 of file MMFSWE.h.

Referenced by v_DoSolve(), and v_SetInitialConditions().

◆ m_Omega

NekDouble Nektar::MMFSWE::m_Omega
protected

Definition at line 95 of file MMFSWE.h.

Referenced by ComputeUnstableJetEta(), EvaluateCoriolisForZonalFlow(), EvaluateStandardCoriolis(), EvaluateWaterDepth(), IsolatedMountainFlow(), RossbyWave(), UnsteadyZonalFlow(), and v_InitObject().

◆ m_planeNumber

int Nektar::MMFSWE::m_planeNumber
protected

Definition at line 115 of file MMFSWE.h.

Referenced by MMFSWE().

◆ m_primitive

bool Nektar::MMFSWE::m_primitive
protected

Indicates if variables are primitive or conservative.

Definition at line 104 of file MMFSWE.h.

◆ m_PurturbedJet

int Nektar::MMFSWE::m_PurturbedJet
private

Definition at line 276 of file MMFSWE.h.

Referenced by UnstableJetFlow(), v_GenerateSummary(), and v_InitObject().

◆ m_RossbyDisturbance

int Nektar::MMFSWE::m_RossbyDisturbance
private

Definition at line 275 of file MMFSWE.h.

Referenced by RossbyWave(), v_GenerateSummary(), and v_InitObject().

◆ m_TestType

TestType Nektar::MMFSWE::m_TestType

Definition at line 79 of file MMFSWE.h.

Referenced by EvaluateCoriolis(), EvaluateWaterDepth(), v_EvaluateExactSolution(), v_GenerateSummary(), v_InitObject(), and v_SetInitialConditions().

◆ m_theta0

NekDouble Nektar::MMFSWE::m_theta0
protected

Definition at line 98 of file MMFSWE.h.

Referenced by ComputeUnstableJetuphi(), and v_InitObject().

◆ m_theta1

NekDouble Nektar::MMFSWE::m_theta1
protected

Definition at line 98 of file MMFSWE.h.

Referenced by ComputeUnstableJetuphi(), and v_InitObject().

◆ m_traceVn

Array<OneD, NekDouble> Nektar::MMFSWE::m_traceVn
protected

Definition at line 108 of file MMFSWE.h.

◆ m_u0

NekDouble Nektar::MMFSWE::m_u0
protected

Definition at line 95 of file MMFSWE.h.

Referenced by IsolatedMountainFlow(), SteadyZonalFlow(), UnsteadyZonalFlow(), and v_InitObject().

◆ m_uthetamax

NekDouble Nektar::MMFSWE::m_uthetamax
protected

Definition at line 98 of file MMFSWE.h.

Referenced by ComputeUnstableJetuphi(), and v_InitObject().

◆ m_veldotMF

Array<OneD, Array<OneD, NekDouble> > Nektar::MMFSWE::m_veldotMF
protected

Definition at line 110 of file MMFSWE.h.

◆ m_vellc

Array<OneD, NekDouble> Nektar::MMFSWE::m_vellc
protected

Definition at line 111 of file MMFSWE.h.

◆ m_velocity

Array<OneD, Array<OneD, NekDouble> > Nektar::MMFSWE::m_velocity
protected

Advection velocity.

Definition at line 107 of file MMFSWE.h.

Referenced by ComputeNablaCdotVelocity().

◆ m_Vorticity0

NekDouble Nektar::MMFSWE::m_Vorticity0
protected

Definition at line 101 of file MMFSWE.h.

Referenced by v_DoSolve(), and v_SetInitialConditions().

Generated on Wed Jan 20 2021 03:33:54 for Nektar++ by   doxygen 1.8.13