Nektar++
Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
Nektar::SolverUtils::MMFSystem Class Reference

A base class for PDEs which include an advection component. More...

#include <MMFSystem.h>

Inheritance diagram for Nektar::SolverUtils::MMFSystem:
[legend]

Public Member Functions

SOLVER_UTILS_EXPORT MMFSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 
virtual SOLVER_UTILS_EXPORT ~MMFSystem ()
 
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &s) override
 Print a summary of time stepping parameters. More...
 
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...
 
SOLVER_UTILS_EXPORT void SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void SetUpTraceNormals (void)
 
SOLVER_UTILS_EXPORT void InitObject (bool DeclareField=true)
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise ()
 Perform any initialisation necessary before solving the problem. More...
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem. More...
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space. More...
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space. More...
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve. More...
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name. More...
 
template<class T >
std::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name. More...
 
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession ()
 Get Session name. More...
 
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure ()
 Get pressure field if available. More...
 
SOLVER_UTILS_EXPORT void ExtraFldOutput (std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
 
SOLVER_UTILS_EXPORT void PrintSummary (std::ostream &out)
 Print a summary of parameters and solver characteristics. More...
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda. More...
 
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction (std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
 Get a SessionFunction by name. More...
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields. More...
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the L2 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDoubleErrorExtraPoints (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::FieldMetaDataMapUpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetFinalTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, const Array< OneD, const NekDouble > &input)
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT int GetInfoSteps ()
 
SOLVER_UTILS_EXPORT void SetInfoSteps (int num)
 
SOLVER_UTILS_EXPORT int GetPararealIterationNumber ()
 
SOLVER_UTILS_EXPORT void SetPararealIterationNumber (int num)
 
SOLVER_UTILS_EXPORT bool GetUseInitialCondition ()
 
SOLVER_UTILS_EXPORT void SetUseInitialCondition (bool num)
 
SOLVER_UTILS_EXPORT Array< OneD, const Array< OneD, NekDouble > > GetTraceNormals ()
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp ()
 Virtual function to identify if operator is negated in DoSolve. More...
 
SOLVER_UTILS_EXPORT bool ParallelInTime ()
 Check if solver use Parallel-in-Time. More...
 

Public Attributes

NekDouble m_pi
 
int m_shapedim
 
SurfaceType m_surfaceType
 
UpwindType m_upwindType
 
TestMaxwellType m_TestMaxwellType
 
PolType m_PolType
 
IncType m_IncType
 
Array< OneD, NekDoublem_MMFfactors
 
- Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 
NekDouble m_cflNonAcoustic
 
NekDouble m_CFLGrowth
 CFL growth rate. More...
 
NekDouble m_CFLEnd
 maximun cfl in cfl growth More...
 

Protected Member Functions

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, NekDoubleCartesianToMovingframes (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, NekDoubleGetIncidentField (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...
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. More...
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoSolve () override
 Solves an unsteady problem. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise () override
 Sets up initial conditions. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_GetTimeStep (const Array< OneD, const Array< OneD, NekDouble >> &inarray)
 Return the timestep to be used for the next step in the time-marching loop. More...
 
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_RequireFwdTrans ()
 
virtual SOLVER_UTILS_EXPORT void v_SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time, int &nchk)
 
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff (const Array< OneD, Array< OneD, NekDouble >> vel, StdRegions::VarCoeffMap &varCoeffMap)
 Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h time velocity. More...
 
virtual SOLVER_UTILS_EXPORT bool v_UpdateTimeStepCheck ()
 
SOLVER_UTILS_EXPORT void DoDummyProjection (const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
 Perform dummy projection. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Virtual function for the L_inf error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the L_2 error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys ()
 Virtual function for transformation to physical space. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff ()
 Virtual function for transformation to coefficient space. More...
 
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution (unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 
virtual SOLVER_UTILS_EXPORT void v_Output (void)
 
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure (void)
 
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
 

Protected Attributes

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, NekDoublem_MFlength
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
int m_nanSteps
 
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
NekDouble m_epsilon
 
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used. More...
 
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used. More...
 
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used. More...
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
NekDouble m_steadyStateRes = 1.0
 
NekDouble m_steadyStateRes0 = 1.0
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check. More...
 
std::ofstream m_errFile
 
std::vector< int > m_intVariables
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
NekDouble m_filterTimeWarning
 Number of time steps between outputting status information. More...
 
NekDouble m_TimeIntegLambda = 0.0
 coefff of spacial derivatives(rhs or m_F in GLM) in calculating the residual of the whole equation(used in unsteady time integrations) More...
 
bool m_flagImplicitItsStatistics
 
bool m_flagImplicitSolver = false
 
Array< OneD, NekDoublem_magnitdEstimat
 estimate the magnitude of each conserved varibles More...
 
Array< OneD, NekDoublem_locTimeStep
 local time step(notice only for jfnk other see m_cflSafetyFactor) More...
 
NekDouble m_inArrayNorm = -1.0
 
int m_TotLinItePerStep = 0
 
int m_StagesPerStep = 1
 
bool m_flagUpdatePreconMat
 
int m_maxLinItePerNewton
 
int m_TotNewtonIts = 0
 
int m_TotLinIts = 0
 
int m_TotImpStages = 0
 
bool m_CalcPhysicalAV = true
 flag to update artificial viscosity More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
bool m_verbose
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtrm_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_fields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_timestepMax = -1.0
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
NekDouble m_lastCheckTime
 
NekDouble m_TimeIncrementFactor
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_infosteps
 Number of time steps between outputting status information. More...
 
int m_pararealIter
 Number of parareal time iteration. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_useInitialCondition
 Flag to determine if IC are used. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
Array< OneD, NekDoublem_movingFrameVelsxyz
 Moving frame of reference velocities. More...
 
Array< OneD, NekDoublem_movingFrameTheta
 Moving frame of reference angles with respect to the. More...
 
boost::numeric::ublas::matrix< NekDoublem_movingFrameProjMat
 Projection matrix for transformation between inertial and moving. More...
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error. More...
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous) More...
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous) More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
int m_HomoDirec
 number of homogenous directions More...
 

Additional Inherited Members

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

Detailed Description

A base class for PDEs which include an advection component.

Definition at line 146 of file MMFSystem.h.

Constructor & Destructor Documentation

◆ MMFSystem()

Nektar::SolverUtils::MMFSystem::MMFSystem ( const LibUtilities::SessionReaderSharedPtr pSession,
const SpatialDomains::MeshGraphSharedPtr pGraph 
)

Definition at line 43 of file MMFSystem.cpp.

45  : UnsteadySystem(pSession, pGraph)
46 {
47 }
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.

◆ ~MMFSystem()

Nektar::SolverUtils::MMFSystem::~MMFSystem ( )
virtual

Definition at line 49 of file MMFSystem.cpp.

50 {
51 }

Member Function Documentation

◆ AbsIntegral()

NekDouble Nektar::SolverUtils::MMFSystem::AbsIntegral ( const Array< OneD, const NekDouble > &  inarray)
protected

Definition at line 2330 of file MMFSystem.cpp.

2331 {
2332  int nq = m_fields[0]->GetNpoints();
2333  Array<OneD, NekDouble> tmp(nq);
2334 
2335  if (inarray.size() != nq)
2336  {
2337  ASSERTL0(false, "AbsIntegral Error: Vector size is not correct");
2338  }
2339 
2340  Vmath::Vabs(nq, inarray, 1, tmp, 1);
2341  return m_fields[0]->Integral(tmp);
2342 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:553

References ASSERTL0, Nektar::SolverUtils::EquationSystem::m_fields, and Vmath::Vabs().

◆ AdddedtMaxwell()

void Nektar::SolverUtils::MMFSystem::AdddedtMaxwell ( const Array< OneD, const Array< OneD, NekDouble >> &  physarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray 
)
protected

Definition at line 1370 of file MMFSystem.cpp.

1373 {
1374  int nq = GetTotPoints();
1375 
1376  // m_dedxi_cdot_e[m][j][n][] = de^m / d \xi^j \cdot e^n
1377  Array<OneD, NekDouble> dedtej(nq);
1378  Array<OneD, NekDouble> de1dtcdotej(nq);
1379  Array<OneD, NekDouble> de2dtcdotej(nq);
1380 
1381  Array<OneD, NekDouble> normH(nq);
1382  Array<OneD, NekDouble> NormedH1(nq, 0.0);
1383  Array<OneD, NekDouble> NormednegH2(nq, 0.0);
1384 
1385  Vmath::Vmul(nq, &physarray[0][0], 1, &physarray[0][0], 1, &normH[0], 1);
1386  Vmath::Vvtvp(nq, &physarray[1][0], 1, &physarray[1][0], 1, &normH[0], 1,
1387  &normH[0], 1);
1388  Vmath::Vsqrt(nq, normH, 1, normH, 1);
1389 
1390  NekDouble Tol = 0.001;
1391  for (int i = 0; i < nq; ++i)
1392  {
1393  if (normH[i] > Tol)
1394  {
1395  NormedH1[i] = physarray[0][i] / normH[i];
1396  NormednegH2[i] = -1.0 * physarray[1][i] / normH[i];
1397  }
1398  }
1399 
1400  for (int j = 0; j < m_shapedim; ++j)
1401  {
1402  // Compute de1 / dt \cdot ej = (-H2 de^1/d\xi1 \cdot e^j + H1 de^1/d\xi2
1403  // \cdot e^j) / sqrt{ H1^2 + H2^2 }
1404  Vmath::Vmul(nq, &NormednegH2[0], 1, &m_dedxi_cdot_e[0][0][j][0], 1,
1405  &de1dtcdotej[0], 1);
1406  Vmath::Vvtvp(nq, &NormedH1[0], 1, &m_dedxi_cdot_e[0][1][j][0], 1,
1407  &de1dtcdotej[0], 1, &de1dtcdotej[0], 1);
1408 
1409  // Compute de2 / dt \cdot ej = (-H2 de2/d\xi1 \cdot e^j + H1 de2/d\xi2
1410  // \cdot e^j) / sqrt{ H1^2 + H2^2 }
1411  Vmath::Vmul(nq, &NormednegH2[0], 1, &m_dedxi_cdot_e[1][0][j][0], 1,
1412  &de2dtcdotej[0], 1);
1413  Vmath::Vvtvp(nq, &NormedH1[0], 1, &m_dedxi_cdot_e[1][1][j][0], 1,
1414  &de2dtcdotej[0], 1, &de2dtcdotej[0], 1);
1415 
1416  // Add dedt component: (H1 (de1/dt) + H2 (de2/dt) ) \cdot ej
1417  Vmath::Vmul(nq, &physarray[0][0], 1, &de1dtcdotej[0], 1, &dedtej[0], 1);
1418  Vmath::Vvtvp(nq, &physarray[1][0], 1, &de2dtcdotej[0], 1, &dedtej[0], 1,
1419  &dedtej[0], 1);
1420 
1421  Vmath::Neg(nq, dedtej, 1);
1422 
1423  switch (m_PolType)
1424  {
1426  {
1427  if (j == 0)
1428  {
1429  Vmath::Vmul(nq, m_muvec[0], 1, dedtej, 1, dedtej, 1);
1430  }
1431 
1432  else if (j == 1)
1433  {
1434  Vmath::Vmul(nq, m_muvec[1], 1, dedtej, 1, dedtej, 1);
1435  }
1436  }
1437  break;
1438 
1440  {
1441  if (j == 0)
1442  {
1443  Vmath::Vmul(nq, m_epsvec[0], 1, dedtej, 1, dedtej, 1);
1444  }
1445 
1446  else if (j == 1)
1447  {
1448  Vmath::Vmul(nq, m_epsvec[1], 1, dedtej, 1, dedtej, 1);
1449  }
1450  }
1451  break;
1452 
1453  default:
1454  break;
1455  }
1456 
1457  Vmath::Vadd(nq, &dedtej[0], 1, &outarray[j][0], 1, &outarray[j][0], 1);
1458  }
1459 }
SOLVER_UTILS_EXPORT int GetTotPoints()
Array< OneD, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > > m_dedxi_cdot_e
Definition: MMFSystem.h:220
Array< OneD, Array< OneD, NekDouble > > m_muvec
Definition: MMFSystem.h:213
Array< OneD, Array< OneD, NekDouble > > m_epsvec
Definition: MMFSystem.h:212
double NekDouble
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
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:209
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
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:574
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:359

References Nektar::SolverUtils::eTransElectric, Nektar::SolverUtils::eTransMagnetic, Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_dedxi_cdot_e, m_epsvec, m_muvec, m_PolType, m_shapedim, Vmath::Neg(), Vmath::Vadd(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vvtvp().

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

◆ AverageMaxwellFlux1D()

void Nektar::SolverUtils::MMFSystem::AverageMaxwellFlux1D ( Array< OneD, Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  numfluxFwd,
Array< OneD, Array< OneD, NekDouble >> &  numfluxBwd 
)
protected

Definition at line 1568 of file MMFSystem.cpp.

1572 {
1573  int i;
1574  int nTraceNumPoints = GetTraceTotPoints();
1575  int nvar = 2;
1576 
1577  // get temporary arrays
1578  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1579  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1580 
1581  for (i = 0; i < nvar; ++i)
1582  {
1583  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1584  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1585  }
1586 
1587  // get the physical values at the trace from the dependent variables
1588  for (i = 0; i < nvar; ++i)
1589  {
1590  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1591  }
1592 
1593  // E = 0 at the boundaries
1594  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1595 
1596  // d H / d n = 0 at the boundaries
1597  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1598 
1599  for (i = 0; i < nTraceNumPoints; ++i)
1600  {
1601  numfluxFwd[0][i] = 0.5 * m_traceNormals[0][i] * (Fwd[1][i] + Bwd[1][i]);
1602  numfluxFwd[1][i] = 0.5 * m_traceNormals[0][i] * (Fwd[0][i] + Bwd[0][i]);
1603 
1604  numfluxBwd[0][i] = 0.5 * m_traceNormals[0][i] * (Fwd[1][i] + Bwd[1][i]);
1605  numfluxBwd[1][i] = 0.5 * m_traceNormals[0][i] * (Fwd[0][i] + Bwd[0][i]);
1606  }
1607 }
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
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

References CopyBoundaryTrace(), Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::eFwdEQNegBwd, Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, and Nektar::SolverUtils::EquationSystem::m_traceNormals.

Referenced by NumericalMaxwellFlux().

◆ AvgAbsInt()

NekDouble Nektar::SolverUtils::MMFSystem::AvgAbsInt ( const Array< OneD, const NekDouble > &  inarray)
protected

Definition at line 2313 of file MMFSystem.cpp.

2314 {
2315  int nq = m_fields[0]->GetNpoints();
2316  Array<OneD, NekDouble> Ones(nq, 1.0);
2317  Array<OneD, NekDouble> tmp(nq);
2318 
2319  if (inarray.size() != nq)
2320  {
2321  ASSERTL0(false, "AvgAbsInt Error: Vector size is not correct");
2322  }
2323 
2324  NekDouble jac = m_fields[0]->Integral(Ones);
2325 
2326  Vmath::Vabs(nq, inarray, 1, tmp, 1);
2327  return (m_fields[0]->Integral(tmp)) / jac;
2328 }

References ASSERTL0, Nektar::SolverUtils::EquationSystem::m_fields, and Vmath::Vabs().

Referenced by CheckMovingFrames(), Nektar::MMFSWE::EvaluateWaterDepth(), and Nektar::MMFSWE::TestVorticityComputation().

◆ AvgInt()

NekDouble Nektar::SolverUtils::MMFSystem::AvgInt ( const Array< OneD, const NekDouble > &  inarray)
protected

Definition at line 2298 of file MMFSystem.cpp.

2299 {
2300  int nq = m_fields[0]->GetNpoints();
2301  Array<OneD, NekDouble> Ones(nq, 1.0);
2302 
2303  if (inarray.size() != nq)
2304  {
2305  ASSERTL0(false, "AvgInt Error: Vector size is not correct");
2306  }
2307 
2308  NekDouble jac = m_fields[0]->Integral(Ones);
2309 
2310  return (m_fields[0]->Integral(inarray)) / jac;
2311 }

References ASSERTL0, and Nektar::SolverUtils::EquationSystem::m_fields.

Referenced by ComputeDivCurlMF(), and Nektar::MMFAdvection::v_InitObject().

◆ BubbleSort()

void Nektar::SolverUtils::MMFSystem::BubbleSort ( Array< OneD, NekDouble > &  refarray,
Array< OneD, NekDouble > &  sortarray 
)
protected

Definition at line 2374 of file MMFSystem.cpp.

2376 {
2377  int nq = refarray.size();
2378 
2379  bool swapped = true;
2380  int j = 0;
2381  NekDouble tmp;
2382 
2383  while (swapped)
2384  {
2385  swapped = false;
2386  j++;
2387  for (int i = 0; i < nq - j; i++)
2388  {
2389  if (refarray[i] > refarray[i + 1])
2390  {
2391  tmp = refarray[i];
2392  refarray[i] = refarray[i + 1];
2393  refarray[i + 1] = tmp;
2394 
2395  tmp = sortarray[i];
2396  sortarray[i] = sortarray[i + 1];
2397  sortarray[i + 1] = tmp;
2398 
2399  swapped = true;
2400  }
2401  }
2402  }
2403 }

◆ CartesianToMovingframes()

Array< OneD, NekDouble > Nektar::SolverUtils::MMFSystem::CartesianToMovingframes ( const Array< OneD, const Array< OneD, NekDouble >> &  uvec,
unsigned int  field 
)
protected

Definition at line 774 of file MMFSystem.cpp.

776 {
777  int nq = m_fields[0]->GetNpoints();
778 
779  Array<OneD, NekDouble> outarray(nq, 0.0);
780 
781  // new u0 = ( [u v] \cdot e^1 )/ |e^1|^2
782  Vmath::Vmul(nq, &m_movingframes[field][0], 1, &uvec[0][0], 1, &outarray[0],
783  1);
784  Vmath::Vvtvp(nq, &m_movingframes[field][nq], 1, &uvec[1][0], 1,
785  &outarray[0], 1, &outarray[0], 1);
786  Vmath::Vvtvp(nq, &m_movingframes[field][2 * nq], 1, &uvec[2][0], 1,
787  &outarray[0], 1, &outarray[0], 1);
788 
789  return outarray;
790 }
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:186

References Nektar::SolverUtils::EquationSystem::m_fields, m_movingframes, Vmath::Vmul(), and Vmath::Vvtvp().

Referenced by Nektar::MMFSWE::IsolatedMountainFlow(), Nektar::MMFSWE::RossbyWave(), Nektar::MMFSWE::SteadyZonalFlow(), Nektar::MMFMaxwell::TestMaxwellSphere(), Nektar::MMFSWE::TestSWE2Dproblem(), Nektar::MMFSWE::TestVorticityComputation(), Nektar::MMFSWE::UnstableJetFlow(), and Nektar::MMFSWE::UnsteadyZonalFlow().

◆ CartesianToSpherical()

void Nektar::SolverUtils::MMFSystem::CartesianToSpherical ( const NekDouble  x0j,
const NekDouble  x1j,
const NekDouble  x2j,
NekDouble sin_varphi,
NekDouble cos_varphi,
NekDouble sin_theta,
NekDouble cos_theta 
)
protected

Definition at line 795 of file MMFSystem.cpp.

799 {
800  NekDouble radius;
801  NekDouble radxy;
802  NekDouble Tol = 0.0000001;
803 
804  NekDouble m_Xscale = 1.0;
805  NekDouble m_Yscale = 1.0;
806  NekDouble m_Zscale = 1.0;
807 
808  radius = sqrt(x0j * x0j / (m_Xscale * m_Xscale) +
809  x1j * x1j / (m_Yscale * m_Yscale) +
810  x2j * x2j / (m_Zscale * m_Zscale));
811  radxy = sqrt(x0j * x0j / (m_Xscale * m_Xscale) +
812  x1j * x1j / (m_Yscale * m_Yscale));
813 
814  if (radxy > Tol)
815  {
816  sin_varphi = x1j / (radxy * m_Yscale);
817  cos_varphi = x0j / (radxy * m_Xscale);
818  }
819 
820  else
821  {
822  sin_varphi = 0.0;
823  if (x2j > 0)
824  {
825  cos_varphi = 1.0;
826  }
827 
828  else
829  {
830  cos_varphi = -1.0;
831  }
832  }
833 
834  sin_theta = x2j / (radius * m_Zscale);
835  cos_theta = radxy / radius;
836 }
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References tinysimd::sqrt().

Referenced by Nektar::MMFAdvection::EvaluateAdvectionVelocity(), Nektar::MMFMaxwell::EvaluateCoriolis(), Nektar::MMFSWE::EvaluateCoriolisForZonalFlow(), Nektar::MMFSWE::EvaluateStandardCoriolis(), Nektar::MMFSWE::EvaluateWaterDepth(), Nektar::MMFSWE::IsolatedMountainFlow(), Nektar::MMFSWE::RossbyWave(), Nektar::MMFSWE::SteadyZonalFlow(), Nektar::MMFMaxwell::TestMaxwellSphere(), Nektar::MMFSWE::TestVorticityComputation(), Nektar::MMFSWE::UnstableJetFlow(), and Nektar::MMFSWE::UnsteadyZonalFlow().

◆ CheckMovingFrames()

void Nektar::SolverUtils::MMFSystem::CheckMovingFrames ( const Array< OneD, const Array< OneD, NekDouble >> &  movingframes)
protected

Definition at line 233 of file MMFSystem.cpp.

235 {
236  NekDouble t1x, t1y, t1z, t2x, t2y, t2z, t3x, t3y, t3z;
237  NekDouble dot12 = 0.0, dot23 = 0.0, dot31 = 0.0;
238  NekDouble Tol = 0.0001;
239 
240  int nq = m_fields[0]->GetNpoints();
241 
242  for (int i = 0; i < nq; ++i)
243  {
244  t1x = movingframes[0][i];
245  t1y = movingframes[0][i + nq];
246  t1z = movingframes[0][i + 2 * nq];
247 
248  t2x = movingframes[1][i];
249  t2y = movingframes[1][i + nq];
250  t2z = movingframes[1][i + 2 * nq];
251 
252  t3x = movingframes[2][i];
253  t3y = movingframes[2][i + nq];
254  t3z = movingframes[2][i + 2 * nq];
255 
256  dot12 = t1x * t2x + t1y * t2y + t1z * t2z;
257  dot23 = t2x * t3x + t2y * t3y + t2z * t3z;
258  dot31 = t3x * t1x + t3y * t1y + t3z * t1z;
259  }
260 
261  std::cout << "======================================================"
262  << std::endl;
263  std::cout << "======================================================"
264  << std::endl;
265  std::cout << "*** The first moving frame is alinged along"
266  << SpatialDomains::GeomMMFMap[m_MMFdir] << std::endl;
267 
268  Array<OneD, NekDouble> tmpx(nq), tmpy(nq), tmpz(nq);
269 
270  Vmath::Vcopy(nq, &movingframes[0][0], 1, &tmpx[0], 1);
271  Vmath::Vcopy(nq, &movingframes[0][nq], 1, &tmpy[0], 1);
272  Vmath::Vcopy(nq, &movingframes[0][2 * nq], 1, &tmpz[0], 1);
273  std::cout << nq << " , "
274  << "*** Avg MF1 = ( " << AvgAbsInt(tmpx) << " , "
275  << AvgAbsInt(tmpy) << " , " << AvgAbsInt(tmpz) << " ) "
276  << std::endl;
277 
278  Vmath::Vcopy(nq, &movingframes[1][0], 1, &tmpx[0], 1);
279  Vmath::Vcopy(nq, &movingframes[1][nq], 1, &tmpy[0], 1);
280  Vmath::Vcopy(nq, &movingframes[1][2 * nq], 1, &tmpz[0], 1);
281  std::cout << "*** Avg MF2 = ( " << AvgAbsInt(tmpx) << " , "
282  << AvgAbsInt(tmpy) << " , " << AvgAbsInt(tmpz) << " ) "
283  << std::endl;
284 
285  if (m_shapedim == 3)
286  {
287  Vmath::Vcopy(nq, &movingframes[2][0], 1, &tmpx[0], 1);
288  Vmath::Vcopy(nq, &movingframes[2][nq], 1, &tmpy[0], 1);
289  Vmath::Vcopy(nq, &movingframes[2][2 * nq], 1, &tmpz[0], 1);
290  std::cout << "*** Avg MF3 = ( " << AvgAbsInt(tmpx) << " , "
291  << AvgAbsInt(tmpy) << " , " << AvgAbsInt(tmpz) << " ) "
292  << std::endl;
293  }
294 
295  if ((fabs(dot12) + fabs(dot23) + fabs(dot31)) < Tol)
296  {
297  std::cout << "*** Moving frames are Orthogonal" << std::endl;
298  }
299 
300  else
301  {
302  std::cout << "*** Moving frames are NOT Orthogonal" << std::endl;
303  }
304 
305  Array<OneD, NekDouble> tmp;
306  for (int j = 0; j < m_shapedim; ++j)
307  {
308  tmp = Array<OneD, NekDouble>(nq, 0.0);
309  for (int k = 0; k < m_spacedim; ++k)
310  {
311  Vmath::Vvtvp(nq, &movingframes[j][k * nq], 1,
312  &movingframes[j][k * nq], 1, &tmp[0], 1, &tmp[0], 1);
313  }
314  Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
315  std::cout << "*** Avg. Magnitude of MF" << j << " = " << AvgAbsInt(tmp)
316  << std::endl;
317  }
318 }
int m_spacedim
Spatial dimension (>= expansion dim).
SpatialDomains::GeomMMF m_MMFdir
Definition: MMFSystem.h:222
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2313
const char *const GeomMMFMap[]
Session file names associated with tangent principle directions.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255

References AvgAbsInt(), Nektar::SpatialDomains::GeomMMFMap, Nektar::SolverUtils::EquationSystem::m_fields, m_MMFdir, m_shapedim, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vcopy(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by SetUpMovingFrames().

◆ ComputeCurl()

void Nektar::SolverUtils::MMFSystem::ComputeCurl ( const Array< OneD, const Array< OneD, NekDouble >> &  inarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray 
)
protected

Definition at line 734 of file MMFSystem.cpp.

738 {
739  int nq = inarray[0].size();
740 
741  Array<OneD, NekDouble> tmpx, tmpy, tmpz;
742  Array<OneD, NekDouble> Dtmpzdx, Dtmpydx, Dtmpxdy, Dtmpzdy, Dtmpxdz, Dtmpydz;
743 
744  tmpx = Array<OneD, NekDouble>(nq);
745  tmpy = Array<OneD, NekDouble>(nq);
746  tmpz = Array<OneD, NekDouble>(nq);
747 
748  Dtmpzdx = Array<OneD, NekDouble>(nq);
749  Dtmpydx = Array<OneD, NekDouble>(nq);
750  Dtmpxdy = Array<OneD, NekDouble>(nq);
751  Dtmpzdy = Array<OneD, NekDouble>(nq);
752  Dtmpxdz = Array<OneD, NekDouble>(nq);
753  Dtmpydz = Array<OneD, NekDouble>(nq);
754 
755  for (int k = 0; k < m_spacedim; ++k)
756  {
757  Vmath::Vcopy(nq, &inarray[0][0], 1, &tmpx[0], 1);
758  Vmath::Vcopy(nq, &inarray[1][0], 1, &tmpy[0], 1);
759  Vmath::Vcopy(nq, &inarray[2][0], 1, &tmpz[0], 1);
760 
761  m_fields[0]->PhysDeriv(0, tmpz, Dtmpzdx);
762  m_fields[0]->PhysDeriv(0, tmpy, Dtmpydx);
763  m_fields[0]->PhysDeriv(1, tmpx, Dtmpxdy);
764  m_fields[0]->PhysDeriv(1, tmpz, Dtmpzdy);
765  m_fields[0]->PhysDeriv(2, tmpx, Dtmpxdz);
766  m_fields[0]->PhysDeriv(2, tmpy, Dtmpydz);
767 
768  Vmath::Vsub(nq, &Dtmpzdy[0], 1, &Dtmpydz[0], 1, &outarray[0][0], 1);
769  Vmath::Vsub(nq, &Dtmpxdz[0], 1, &Dtmpzdx[0], 1, &outarray[1][0], 1);
770  Vmath::Vsub(nq, &Dtmpydx[0], 1, &Dtmpxdy[0], 1, &outarray[2][0], 1);
771  }
772 }
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:419

References Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vcopy(), and Vmath::Vsub().

Referenced by ComputeDivCurlMF().

◆ Computedemdxicdote()

void Nektar::SolverUtils::MMFSystem::Computedemdxicdote ( )
protected

Definition at line 1280 of file MMFSystem.cpp.

1281 {
1282  int MFdim = 3;
1283  int nq = GetTotPoints();
1284 
1285  // Initialization
1286  m_dedxi_cdot_e =
1287  Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble>>>>(MFdim);
1288  for (int indm = 0; indm < MFdim; ++indm)
1289  {
1290  m_dedxi_cdot_e[indm] =
1291  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
1292  for (int indj = 0; indj < MFdim; ++indj)
1293  {
1294  m_dedxi_cdot_e[indm][indj] =
1295  Array<OneD, Array<OneD, NekDouble>>(MFdim);
1296  for (int indn = 0; indn < MFdim; ++indn)
1297  {
1298  m_dedxi_cdot_e[indm][indj][indn] =
1299  Array<OneD, NekDouble>(nq, 0.0);
1300  }
1301  }
1302  }
1303 
1304  Array<OneD, NekDouble> tmp(nq);
1305  Array<OneD, NekDouble> tmpderiv(nq);
1306  Array<OneD, NekDouble> dedt;
1307  for (int indm = 0; indm < MFdim; ++indm)
1308  {
1309  for (int indj = 0; indj < MFdim; ++indj)
1310  {
1311  for (int indn = 0; indn < MFdim; ++indn)
1312  {
1313  dedt = Array<OneD, NekDouble>(nq, 0.0);
1314  for (int k = 0; k < m_spacedim; ++k)
1315  {
1316  // Compute d e^m / d \xi_j cdot e^n
1317  Vmath::Vcopy(nq, &m_movingframes[indm][k * nq], 1, &tmp[0],
1318  1);
1319  m_fields[0]->PhysDirectionalDeriv(m_movingframes[indj], tmp,
1320  tmpderiv);
1321 
1322  Vmath::Vvtvp(nq, &tmpderiv[0], 1,
1323  &m_movingframes[indn][k * nq], 1, &dedt[0], 1,
1324  &dedt[0], 1);
1325  }
1326 
1327  Vmath::Vcopy(nq, &dedt[0], 1,
1328  &m_dedxi_cdot_e[indm][indj][indn][0], 1);
1329  }
1330  }
1331  }
1332 
1333  int indx = 0;
1334  std::cout << "*** m_dedxi_cdot_e[0]/dxi1 = ( "
1335  << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1336  << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1337  << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1338  << std::endl;
1339  std::cout << "*** m_dedxi_cdot_e[0]/dxi2 = ( "
1340  << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1341  << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1342  << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1343  << std::endl;
1344 
1345  indx = 1;
1346  std::cout << "*** m_dedxi_cdot_e[1]/dxi1 = ( "
1347  << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1348  << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1349  << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1350  << std::endl;
1351  std::cout << "*** m_dedxi_cdot_e[1]/dxi2 = ( "
1352  << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1353  << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1354  << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1355  << std::endl;
1356 
1357  indx = 2;
1358  std::cout << "*** m_dedxi_cdot_e[2]/dxi1 = ( "
1359  << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1360  << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1361  << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1362  << std::endl;
1363  std::cout << "*** m_dedxi_cdot_e[2]/dxi2 = ( "
1364  << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1365  << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1366  << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1367  << std::endl;
1368 }
SOLVER_UTILS_EXPORT NekDouble RootMeanSquare(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2344

References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_dedxi_cdot_e, Nektar::SolverUtils::EquationSystem::m_fields, m_movingframes, Nektar::SolverUtils::EquationSystem::m_spacedim, RootMeanSquare(), Vmath::Vcopy(), and Vmath::Vvtvp().

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

◆ ComputeDivCurlMF()

void Nektar::SolverUtils::MMFSystem::ComputeDivCurlMF ( )
protected

Definition at line 445 of file MMFSystem.cpp.

446 {
447  int nq = m_fields[0]->GetNpoints();
448  int MMFdim = 3;
449 
450  Array<OneD, NekDouble> tmp(nq);
451  Array<OneD, NekDouble> Dtmp(nq);
452 
453  m_DivMF = Array<OneD, Array<OneD, NekDouble>>(MMFdim);
454  for (int j = 0; j < MMFdim; ++j)
455  {
456  m_DivMF[j] = Array<OneD, NekDouble>(nq, 0.0);
457  for (int k = 0; k < m_spacedim; ++k)
458  {
459  Vmath::Vcopy(nq, &m_movingframes[j][k * nq], 1, &tmp[0], 1);
460 
461  m_fields[0]->PhysDeriv(k, tmp, Dtmp);
462  Vmath::Vadd(nq, &Dtmp[0], 1, &m_DivMF[j][0], 1, &m_DivMF[j][0], 1);
463  }
464  }
465 
466  std::cout << "*** Divergence of MF1 = " << AvgInt(m_DivMF[0])
467  << ", MF2 = " << AvgInt(m_DivMF[1])
468  << ", MF3 = " << AvgInt(m_DivMF[2]) << std::endl;
469 
470  // Compute Curl of MF: CurlMF[i][j] = (\nabla \times e^i) cdot e^j
471  m_CurlMF = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MMFdim);
472  for (int i = 0; i < MMFdim; ++i)
473  {
474  m_CurlMF[i] = Array<OneD, Array<OneD, NekDouble>>(MMFdim);
475 
476  for (int j = 0; j < MMFdim; ++j)
477  {
478  m_CurlMF[i][j] = Array<OneD, NekDouble>(nq, 0.0);
479  }
480  }
481 
482  Array<OneD, Array<OneD, NekDouble>> MFtmp(m_spacedim);
483  Array<OneD, Array<OneD, NekDouble>> CurlMFtmp(m_spacedim);
484 
485  for (int i = 0; i < m_spacedim; ++i)
486  {
487  MFtmp[i] = Array<OneD, NekDouble>(nq);
488  CurlMFtmp[i] = Array<OneD, NekDouble>(nq);
489  }
490 
491  for (int dir = 0; dir < MMFdim; dir++)
492  {
493  for (int i = 0; i < m_spacedim; ++i)
494  {
495  Vmath::Vcopy(nq, &m_movingframes[dir][i * nq], 1, &MFtmp[i][0], 1);
496  }
497 
498  ComputeCurl(MFtmp, CurlMFtmp);
499 
500  for (int j = 0; j < MMFdim; ++j)
501  {
502  for (int i = 0; i < m_spacedim; ++i)
503  {
504  Vmath::Vvtvp(nq, &m_movingframes[j][i * nq], 1,
505  &CurlMFtmp[i][0], 1, &m_CurlMF[dir][j][0], 1,
506  &m_CurlMF[dir][j][0], 1);
507  }
508  }
509  }
510 
511  std::cout << "*** Curl of MF1 = ( " << AvgInt(m_CurlMF[0][0]) << " , "
512  << AvgInt(m_CurlMF[0][1]) << " , " << AvgInt(m_CurlMF[0][2])
513  << " ) " << std::endl;
514  std::cout << "*** Curl of MF2 = ( " << AvgInt(m_CurlMF[1][0]) << " , "
515  << AvgInt(m_CurlMF[1][1]) << " , " << AvgInt(m_CurlMF[1][2])
516  << " ) " << std::endl;
517  std::cout << "*** Curl of MF3 = ( " << AvgInt(m_CurlMF[2][0]) << " , "
518  << AvgInt(m_CurlMF[2][1]) << " , " << AvgInt(m_CurlMF[2][2])
519  << " ) " << std::endl;
520 }
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Definition: MMFSystem.h:195
SOLVER_UTILS_EXPORT NekDouble AvgInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2298
SOLVER_UTILS_EXPORT void ComputeCurl(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: MMFSystem.cpp:734
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Definition: MMFSystem.h:196

References AvgInt(), ComputeCurl(), m_CurlMF, m_DivMF, Nektar::SolverUtils::EquationSystem::m_fields, m_movingframes, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vvtvp().

Referenced by MMFInitObject().

◆ ComputeMFtrace()

void Nektar::SolverUtils::MMFSystem::ComputeMFtrace ( )
protected

Definition at line 522 of file MMFSystem.cpp.

523 {
524  int MFdim = 3;
525 
526  int nq = m_fields[0]->GetNpoints();
527  int nTraceNumPoints = GetTraceTotPoints();
528 
529  Array<OneD, NekDouble> tmp(nq);
530  Array<OneD, NekDouble> Fwdtmp(nq);
531  Array<OneD, NekDouble> Bwdtmp(nq);
532 
533  m_MFtraceFwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
534  m_MFtraceBwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
535 
536  for (int j = 0; j < MFdim; ++j)
537  {
538  m_MFtraceFwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
539  m_MFtraceBwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
540  }
541 
542  // m_MFtraceFwd[0] = e^1_{Fwd}, m_MFtraceFwd[1] = e^2_{Fwd}
543  for (int j = 0; j < MFdim; ++j)
544  {
545  for (int i = 0; i < m_spacedim; ++i)
546  {
547  m_MFtraceFwd[j][i] = Array<OneD, NekDouble>(nTraceNumPoints);
548  m_MFtraceBwd[j][i] = Array<OneD, NekDouble>(nTraceNumPoints);
549 
550  Vmath::Vcopy(nq, &m_movingframes[j][i * nq], 1, &tmp[0], 1);
551 
552  m_fields[0]->GetFwdBwdTracePhys(tmp, Fwdtmp, Bwdtmp);
553 
554  CopyBoundaryTrace(Fwdtmp, Bwdtmp, eFwdEQBwd);
555 
556  Vmath::Vcopy(nTraceNumPoints, &Fwdtmp[0], 1, &m_MFtraceFwd[j][i][0],
557  1);
558  Vmath::Vcopy(nTraceNumPoints, &Bwdtmp[0], 1, &m_MFtraceBwd[j][i][0],
559  1);
560  }
561  }
562 
563  std::cout << "*** MFtraceFwd = ( " << VectorAvgMagnitude(m_MFtraceFwd[0])
564  << " , " << VectorAvgMagnitude(m_MFtraceFwd[1]) << " , "
565  << VectorAvgMagnitude(m_MFtraceFwd[2]) << " ) " << std::endl;
566  std::cout << "*** MFtraceBwd = ( " << VectorAvgMagnitude(m_MFtraceBwd[0])
567  << " , " << VectorAvgMagnitude(m_MFtraceBwd[1]) << " , "
568  << VectorAvgMagnitude(m_MFtraceBwd[2]) << " ) " << std::endl;
569 }
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
Definition: MMFSystem.h:199
SOLVER_UTILS_EXPORT NekDouble VectorAvgMagnitude(const Array< OneD, const Array< OneD, NekDouble >> &inarray)
Definition: MMFSystem.cpp:2358
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
Definition: MMFSystem.h:200

References CopyBoundaryTrace(), Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_MFtraceBwd, m_MFtraceFwd, m_movingframes, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vcopy(), and VectorAvgMagnitude().

Referenced by MMFInitObject().

◆ ComputencdotMF()

void Nektar::SolverUtils::MMFSystem::ComputencdotMF ( )
protected

Definition at line 320 of file MMFSystem.cpp.

321 {
322  int nq = m_fields[0]->GetNpoints();
323  int nTracePointsTot = GetTraceNpoints();
324 
325  // Compute MFjFwd and MFjBwd
326  Array<OneD, NekDouble> tmp(nq);
327 
328  Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtraceFwd(m_shapedim);
329  Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtraceBwd(m_shapedim);
330  Array<OneD, Array<OneD, NekDouble>> SurfaceNormalFwd;
331  Array<OneD, Array<OneD, NekDouble>> SurfaceNormalBwd;
332 
333  for (int j = 0; j < m_shapedim; ++j)
334  {
335  MFtraceFwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
336  MFtraceBwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
337 
338  SurfaceNormalFwd = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
339  SurfaceNormalBwd = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
340 
341  for (int k = 0; k < m_spacedim; ++k)
342  {
343  MFtraceFwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
344  MFtraceBwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
345 
346  SurfaceNormalFwd[k] = Array<OneD, NekDouble>(nTracePointsTot);
347  SurfaceNormalBwd[k] = Array<OneD, NekDouble>(nTracePointsTot);
348 
349  Vmath::Vcopy(nq, &m_movingframes[j][k * nq], 1, &tmp[0], 1);
350 
351  m_fields[0]->GetFwdBwdTracePhys(tmp, MFtraceFwd[j][k],
352  MFtraceBwd[j][k]);
353 
354  CopyBoundaryTrace(MFtraceFwd[j][k], MFtraceBwd[j][k], eFwdEQBwd);
355  }
356  }
357 
358  VectorCrossProd(MFtraceFwd[0], MFtraceFwd[1], SurfaceNormalFwd);
359  VectorCrossProd(MFtraceBwd[0], MFtraceBwd[1], SurfaceNormalBwd);
360 
361  // Compute n \times e^i
362  m_ncdotMFFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
363  m_ncdotMFBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
364  for (int j = 0; j < m_shapedim; ++j)
365  {
366  m_ncdotMFFwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
367  m_ncdotMFBwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
368 
369  VectorDotProd(m_traceNormals, MFtraceFwd[j], m_ncdotMFFwd[j]);
370  VectorDotProd(m_traceNormals, MFtraceBwd[j], m_ncdotMFBwd[j]);
371  }
372 
373  // Compute n^{\perp} \times e^i
374  Array<OneD, Array<OneD, NekDouble>> SurfaceNormal(m_spacedim);
375  Array<OneD, Array<OneD, NekDouble>> Tracevector(m_spacedim);
376  for (int k = 0; k < m_spacedim; k++)
377  {
378  SurfaceNormal[k] = Array<OneD, NekDouble>(nTracePointsTot);
379  Tracevector[k] = Array<OneD, NekDouble>(nTracePointsTot);
380 
381  Vmath::Vcopy(nTracePointsTot, &m_MFtraceFwd[2][k][0], 1,
382  &SurfaceNormal[k][0], 1);
383  }
384 
385  VectorCrossProd(m_traceNormals, SurfaceNormal, Tracevector);
386 
387  m_nperpcdotMFFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
388  m_nperpcdotMFBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
389  for (int j = 0; j < m_shapedim; ++j)
390  {
391  m_nperpcdotMFFwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
392  m_nperpcdotMFBwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
393 
394  VectorDotProd(Tracevector, MFtraceFwd[j], m_nperpcdotMFFwd[j]);
395  VectorDotProd(Tracevector, MFtraceBwd[j], m_nperpcdotMFBwd[j]);
396  }
397 
398  if (m_shapedim == 2)
399  {
400  std::cout << "*** m_ncdotMFFwd = ( " << RootMeanSquare(m_ncdotMFFwd[0])
401  << " , " << RootMeanSquare(m_ncdotMFFwd[1]) << " ) "
402  << std::endl;
403  std::cout << "*** m_ncdotMFBwd = ( " << RootMeanSquare(m_ncdotMFBwd[0])
404  << " , " << RootMeanSquare(m_ncdotMFBwd[1]) << " ) "
405  << std::endl;
406 
407  std::cout << "*** m_nperpcdotMFFwd = ( "
408  << RootMeanSquare(m_nperpcdotMFFwd[0]) << " , "
409  << RootMeanSquare(m_nperpcdotMFFwd[1]) << " ) " << std::endl;
410  std::cout << "*** m_nperpcdotMFBwd = ( "
411  << RootMeanSquare(m_nperpcdotMFBwd[0]) << " , "
412  << RootMeanSquare(m_nperpcdotMFBwd[1]) << " ) " << std::endl;
413  }
414 
415  else if (m_shapedim == 3)
416  {
417  std::cout << "*** m_ncdotMFFwd = ( "
418  << Vmath::Vsum(nTracePointsTot, m_ncdotMFFwd[0], 1) << " , "
419  << Vmath::Vsum(nTracePointsTot, m_ncdotMFFwd[1], 1) << " , "
420  << Vmath::Vsum(nTracePointsTot, m_ncdotMFFwd[2], 1) << " ) "
421  << std::endl;
422  std::cout << "*** m_ncdotMFBwd = ( "
423  << Vmath::Vsum(nTracePointsTot, m_ncdotMFBwd[0], 1) << " , "
424  << Vmath::Vsum(nTracePointsTot, m_ncdotMFBwd[1], 1) << " , "
425  << Vmath::Vsum(nTracePointsTot, m_ncdotMFBwd[2], 1) << " ) "
426  << std::endl;
427 
428  std::cout << "*** m_nperpcdotMFFwd = ( "
429  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFFwd[0], 1)
430  << " , "
431  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFFwd[1], 1)
432  << " , "
433  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFFwd[2], 1)
434  << " ) " << std::endl;
435  std::cout << "*** m_nperpcdotMFBwd = ( "
436  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFBwd[0], 1)
437  << " , "
438  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFBwd[1], 1)
439  << " , "
440  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFBwd[2], 1)
441  << " ) " << std::endl;
442  }
443 }
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
Definition: MMFSystem.h:192
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)
Definition: MMFSystem.cpp:699
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFBwd
Definition: MMFSystem.h:193
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)
Definition: MMFSystem.cpp:676
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:190
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:189
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:895

References CopyBoundaryTrace(), Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::EquationSystem::GetTraceNpoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_MFtraceFwd, m_movingframes, m_ncdotMFBwd, m_ncdotMFFwd, m_nperpcdotMFBwd, m_nperpcdotMFFwd, m_shapedim, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_traceNormals, RootMeanSquare(), Vmath::Vcopy(), VectorCrossProd(), VectorDotProd(), and Vmath::Vsum().

Referenced by MMFInitObject(), and Nektar::MMFAdvection::v_InitObject().

◆ ComputeNtimesF12()

void Nektar::SolverUtils::MMFSystem::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 
)
protected

Definition at line 951 of file MMFSystem.cpp.

959 {
960  int nTraceNumPoints = GetTraceTotPoints();
961 
962  NekDouble tmpFwd, tmpBwd, Aver1, Aver2, HFwdk, HBwdk;
963 
964  Array<OneD, NekDouble> z1HAver(m_spacedim);
965  Array<OneD, NekDouble> z2HAver(m_spacedim);
966  Array<OneD, NekDouble> n1e1(m_spacedim);
967  Array<OneD, NekDouble> n2e2(m_spacedim);
968 
969  Array<OneD, NekDouble> n1e1_times_z1HAver(m_spacedim);
970  Array<OneD, NekDouble> n2e2_times_z2HAver(m_spacedim);
971 
972  for (int i = 0; i < nTraceNumPoints; ++i)
973  {
974  Aver1 = 0.5 * (im1Fwd[i] + im1Bwd[i]);
975  Aver2 = 0.5 * (im2Fwd[i] + im2Bwd[i]);
976 
977  for (int k = 0; k < m_spacedim; k++)
978  {
979  // Compute \vec{HFwd} and \vec{HBwd}
980  HFwdk = Fwd[0][i] * m_MFtraceFwd[0][k][i] +
981  Fwd[1][i] * m_MFtraceFwd[1][k][i];
982  HBwdk = Bwd[0][i] * m_MFtraceBwd[0][k][i] +
983  Bwd[1][i] * m_MFtraceBwd[1][k][i];
984 
985  // Compute z_i {{ \vec{H} }}
986  z1HAver[k] = 0.5 * (im1Fwd[i] * HFwdk + im1Bwd[i] * HBwdk);
987  z2HAver[k] = 0.5 * (im2Fwd[i] * HFwdk + im2Bwd[i] * HBwdk);
988 
989  // Choose e^i for the one in anisotropy region
990  n1e1[k] = m_ncdotMFFwd[0][i] * m_MFtraceFwd[0][k][i];
991  n2e2[k] = m_ncdotMFFwd[1][i] * m_MFtraceFwd[1][k][i];
992  }
993 
994  // Compute n1e1 \times z1HAver and n2e2 \times z2HAver
995  VectorCrossProd(n1e1, z1HAver, n1e1_times_z1HAver);
996  VectorCrossProd(n2e2, z2HAver, n2e2_times_z2HAver);
997 
998  // e^3 \cdot ( n1e1 \times z1HAver + n2e2 \times z2HAver)
999  tmpFwd = 0.0;
1000  tmpBwd = 0.0;
1001  for (int k = 0; k < m_spacedim; k++)
1002  {
1003  tmpFwd += m_MFtraceFwd[2][k][i] * (n1e1_times_z1HAver[k] / Aver1 +
1004  n2e2_times_z2HAver[k] / Aver2);
1005  tmpBwd += m_MFtraceBwd[2][k][i] * (n1e1_times_z1HAver[k] / Aver1 +
1006  n2e2_times_z2HAver[k] / Aver2);
1007  }
1008 
1009  outarrayFwd[i] = tmpFwd;
1010  outarrayBwd[i] = tmpBwd;
1011  }
1012 }

References Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), m_MFtraceBwd, m_MFtraceFwd, m_ncdotMFFwd, Nektar::SolverUtils::EquationSystem::m_spacedim, and VectorCrossProd().

Referenced by NumericalMaxwellFluxTE(), and NumericalMaxwellFluxTM().

◆ ComputeNtimesFz()

void Nektar::SolverUtils::MMFSystem::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 
)
protected

Definition at line 915 of file MMFSystem.cpp.

923 {
924  int nTraceNumPoints = GetTraceTotPoints();
925 
926  NekDouble tmpFwd, tmpBwd;
927  NekDouble Aver, ntimesz;
928 
929  for (int i = 0; i < nTraceNumPoints; ++i)
930  {
931  Aver = 0.5 * (imFwd[i] + imBwd[i]);
932 
933  tmpFwd = 0.0;
934  tmpBwd = 0.0;
935  for (int k = 0; k < m_spacedim; ++k)
936  {
937  ntimesz = 0.5 * (imFwd[i] * Fwd[2][i] + imBwd[i] * Bwd[2][i]);
938 
939  tmpFwd +=
940  m_MFtraceFwd[dir][k][i] * m_ntimesMFFwd[2][k][i] * ntimesz;
941  tmpBwd +=
942  m_MFtraceBwd[dir][k][i] * m_ntimesMFBwd[2][k][i] * ntimesz;
943  }
944 
945  outarrayFwd[i] = tmpFwd / Aver;
946  outarrayBwd[i] = tmpBwd / Aver;
947  }
948 }
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFBwd
Definition: MMFSystem.h:203
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFFwd
Definition: MMFSystem.h:202

References Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), m_MFtraceBwd, m_MFtraceFwd, m_ntimesMFBwd, m_ntimesMFFwd, and Nektar::SolverUtils::EquationSystem::m_spacedim.

Referenced by NumericalMaxwellFluxTE(), and NumericalMaxwellFluxTM().

◆ ComputeNtimesMF()

void Nektar::SolverUtils::MMFSystem::ComputeNtimesMF ( )
protected

Definition at line 618 of file MMFSystem.cpp.

619 {
620  int MFdim = 3;
621  int nTracePointsTot = GetTraceNpoints();
622 
623  m_ntimesMFFwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
624  m_ntimesMFBwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
626  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
628  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
629  for (int j = 0; j < MFdim; ++j)
630  {
631  m_ntimesMFFwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
632  m_ntimesMFBwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
634  Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
636  Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
637 
638  for (int k = 0; k < m_spacedim; ++k)
639  {
640  m_ntimesMFFwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
641  m_ntimesMFBwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
642  m_ntimes_ntimesMFFwd[j][k] =
643  Array<OneD, NekDouble>(nTracePointsTot);
644  m_ntimes_ntimesMFBwd[j][k] =
645  Array<OneD, NekDouble>(nTracePointsTot);
646  }
647 
654  }
655 
656  std::cout << "*** m_ntimesMFFwd = ( "
657  << VectorAvgMagnitude(m_ntimesMFFwd[0]) << " , "
658  << VectorAvgMagnitude(m_ntimesMFFwd[1]) << " , "
659  << VectorAvgMagnitude(m_ntimesMFFwd[2]) << " ) " << std::endl;
660  std::cout << "*** m_ntimesMFBwd = ( "
661  << VectorAvgMagnitude(m_ntimesMFBwd[0]) << " , "
662  << VectorAvgMagnitude(m_ntimesMFBwd[1]) << " , "
663  << VectorAvgMagnitude(m_ntimesMFBwd[2]) << " ) " << std::endl;
664  std::cout << "*** m_ntimes_ntimesMFFwd = ( "
668  << std::endl;
669  std::cout << "*** m_ntimes_ntimesMFBwd = ( "
673  << std::endl;
674 }
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFBwd
Definition: MMFSystem.h:205
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFFwd
Definition: MMFSystem.h:204

References Nektar::SolverUtils::EquationSystem::GetTraceNpoints(), m_MFtraceBwd, m_MFtraceFwd, m_ntimes_ntimesMFBwd, m_ntimes_ntimesMFFwd, m_ntimesMFBwd, m_ntimesMFFwd, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_traceNormals, VectorAvgMagnitude(), and VectorCrossProd().

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

◆ ComputeNtimestimesdF12()

void Nektar::SolverUtils::MMFSystem::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 
)
protected

Definition at line 1080 of file MMFSystem.cpp.

1088 {
1089  int nTraceNumPoints = GetTraceTotPoints();
1090 
1091  Array<OneD, NekDouble> directFwd(nTraceNumPoints);
1092  Array<OneD, NekDouble> directBwd(nTraceNumPoints);
1093 
1094  NekDouble Aver1, Aver2;
1095  for (int i = 0; i < nTraceNumPoints; ++i)
1096  {
1097  Aver1 = im1Fwd[i] + im1Bwd[i];
1098  Aver2 = im2Fwd[i] + im2Bwd[i];
1099 
1100  outarrayFwd[i] =
1101  -m_alpha * (Fwd[2][i] - Bwd[2][i]) * (1.0 / Aver1 + 1.0 / Aver2);
1102  outarrayBwd[i] =
1103  -m_alpha * (Fwd[2][i] - Bwd[2][i]) * (1.0 / Aver1 + 1.0 / Aver2);
1104  }
1105 }

References Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), and m_alpha.

Referenced by NumericalMaxwellFluxTE(), and NumericalMaxwellFluxTM().

◆ ComputeNtimestimesdFz()

void Nektar::SolverUtils::MMFSystem::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 
)
protected

Definition at line 1014 of file MMFSystem.cpp.

1020 {
1021  int nTraceNumPoints = GetTraceTotPoints();
1022 
1023  Array<OneD, NekDouble> dH(m_spacedim);
1024  Array<OneD, NekDouble> nFwd(m_spacedim);
1025 
1026  Array<OneD, NekDouble> eiFwd(m_spacedim);
1027  Array<OneD, NekDouble> eiBwd(m_spacedim);
1028 
1029  Array<OneD, NekDouble> eitimesdHFwd(m_spacedim);
1030  Array<OneD, NekDouble> eitimesdHBwd(m_spacedim);
1031 
1032  Array<OneD, NekDouble> ntimeseitimesdHFwd(m_spacedim);
1033  Array<OneD, NekDouble> ntimeseitimesdHBwd(m_spacedim);
1034 
1035  NekDouble Aver, HFwdk, HBwdk, tmpFwd, tmpBwd;
1036  for (int i = 0; i < nTraceNumPoints; ++i)
1037  {
1038  Aver = 0.5 * (imFwd[i] + imBwd[i]);
1039 
1040  // Get [H]
1041  for (int k = 0; k < m_spacedim; k++)
1042  {
1043  HFwdk = Fwd[0][i] * m_MFtraceFwd[0][k][i] +
1044  Fwd[1][i] * m_MFtraceFwd[1][k][i];
1045  HBwdk = Bwd[0][i] * m_MFtraceBwd[0][k][i] +
1046  Bwd[1][i] * m_MFtraceBwd[1][k][i];
1047  dH[k] = HFwdk - HBwdk;
1048 
1049  eiFwd[k] = m_MFtraceFwd[dir][k][i];
1050  eiBwd[k] = m_MFtraceBwd[dir][k][i];
1051 
1052  nFwd[k] = m_traceNormals[k][i];
1053  }
1054 
1055  // MFtraceFwd (MFtraceBwd) \times [H]
1056  // VectorCrossProd(eiFwd, dH, eitimesdHFwd);
1057  // VectorCrossProd(eiBwd, dH, eitimesdHBwd);
1058  VectorCrossProd(nFwd, dH, eitimesdHFwd);
1059  VectorCrossProd(nFwd, dH, eitimesdHBwd);
1060 
1061  // n times eitimesdH
1062  VectorCrossProd(nFwd, eitimesdHFwd, ntimeseitimesdHFwd);
1063  VectorCrossProd(nFwd, eitimesdHBwd, ntimeseitimesdHBwd);
1064 
1065  // MFtraceFwd \cdot ntimeseitimesdH
1066  tmpFwd = 0.0;
1067  tmpBwd = 0.0;
1068  for (int k = 0; k < m_spacedim; k++)
1069  {
1070  tmpFwd += eiFwd[k] * ntimeseitimesdHFwd[k];
1071  tmpBwd += eiBwd[k] * ntimeseitimesdHBwd[k];
1072  }
1073 
1074  outarrayFwd[i] = 0.5 * m_alpha * tmpFwd / Aver;
1075  outarrayBwd[i] = 0.5 * m_alpha * tmpBwd / Aver;
1076  }
1077 }

References Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), m_alpha, m_MFtraceBwd, m_MFtraceFwd, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_traceNormals, and VectorCrossProd().

Referenced by NumericalMaxwellFluxTE(), and NumericalMaxwellFluxTM().

◆ ComputeZimYim()

void Nektar::SolverUtils::MMFSystem::ComputeZimYim ( Array< OneD, Array< OneD, NekDouble >> &  epsvec,
Array< OneD, Array< OneD, NekDouble >> &  muvec 
)
protected

Definition at line 1108 of file MMFSystem.cpp.

1111 {
1112  int nTraceNumPoints = GetTraceNpoints();
1113 
1114  switch (m_TestMaxwellType)
1115  {
1116  case eMaxwell1D:
1117  case eScatField1D:
1118  {
1119  Array<OneD, NekDouble> Fwdeps(nTraceNumPoints, 1.0);
1120  Array<OneD, NekDouble> Bwdeps(nTraceNumPoints, 1.0);
1121  Array<OneD, NekDouble> Fwdmu(nTraceNumPoints, 1.0);
1122  Array<OneD, NekDouble> Bwdmu(nTraceNumPoints, 1.0);
1123  m_fields[0]->GetFwdBwdTracePhys(epsvec[0], Fwdeps, Bwdeps);
1124  m_fields[0]->GetFwdBwdTracePhys(muvec[0], Fwdeps, Bwdeps);
1125 
1126  CopyBoundaryTrace(Fwdeps, Bwdeps, eFwdEQBwd, 0);
1127  CopyBoundaryTrace(Fwdmu, Bwdmu, eFwdEQBwd, 1);
1128 
1129  m_ZimFwd = Array<OneD, Array<OneD, NekDouble>>(1);
1130  m_ZimBwd = Array<OneD, Array<OneD, NekDouble>>(1);
1131  m_YimFwd = Array<OneD, Array<OneD, NekDouble>>(1);
1132  m_YimBwd = Array<OneD, Array<OneD, NekDouble>>(1);
1133 
1134  m_ZimFwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1135  m_ZimBwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1136  m_YimFwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1137  m_YimBwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1138 
1139  // ZimFwd = sqrt( muFwd / epsFwd), ZimBwd = sqrt( muBwd / epsBwd)
1140  for (int i = 0; i < nTraceNumPoints; ++i)
1141  {
1142  m_ZimFwd[0][i] = sqrt(Fwdmu[i] / Fwdeps[i]);
1143  m_ZimBwd[0][i] = sqrt(Bwdmu[i] / Bwdeps[i]);
1144 
1145  m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1146  m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1147  }
1148 
1149  std::cout << "*** ZimFwd = " << RootMeanSquare(m_ZimFwd[0])
1150  << ", ZimBwd = " << RootMeanSquare(m_ZimBwd[0])
1151  << ", YimFwd = " << RootMeanSquare(m_YimFwd[0])
1152  << ", YimBwd = " << RootMeanSquare(m_YimBwd[0])
1153  << std::endl;
1154  }
1155  break;
1156 
1157  case eTestMaxwell2DPEC:
1159  case eTestMaxwell2DPMC:
1160  case eScatField2D:
1161  case eTotField2D:
1162  case eMaxwellSphere:
1163  case eELF2DSurface:
1164  {
1165  m_ZimFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1166  m_ZimBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1167  m_YimFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1168  m_YimBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1169 
1170  for (int j = 0; j < m_shapedim; ++j)
1171  {
1172  m_ZimFwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1173  m_ZimBwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1174  m_YimFwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1175  m_YimBwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1176  }
1177 
1178  switch (m_PolType)
1179  {
1180  case eTransMagnetic:
1181  {
1182  Array<OneD, NekDouble> Fwdmu1(nTraceNumPoints);
1183  Array<OneD, NekDouble> Bwdmu1(nTraceNumPoints);
1184  Array<OneD, NekDouble> Fwdmu2(nTraceNumPoints);
1185  Array<OneD, NekDouble> Bwdmu2(nTraceNumPoints);
1186  Array<OneD, NekDouble> Fwdeps3(nTraceNumPoints);
1187  Array<OneD, NekDouble> Bwdeps3(nTraceNumPoints);
1188 
1189  m_fields[0]->GetFwdBwdTracePhys(muvec[0], Fwdmu1, Bwdmu1);
1190  m_fields[0]->GetFwdBwdTracePhys(muvec[1], Fwdmu2, Bwdmu2);
1191  m_fields[0]->GetFwdBwdTracePhys(epsvec[2], Fwdeps3,
1192  Bwdeps3);
1193 
1194  CopyBoundaryTrace(Fwdmu1, Bwdmu1, eFwdEQBwd, 0);
1195  CopyBoundaryTrace(Fwdmu2, Bwdmu2, eFwdEQBwd, 1);
1196  CopyBoundaryTrace(Fwdeps3, Bwdeps3, eFwdEQBwd, 2);
1197 
1198  // ZimFwd = sqrt( muFwd / epsFwd), ZimBwd = sqrt( muBwd /
1199  // epsBwd)
1200  for (int i = 0; i < nTraceNumPoints; ++i)
1201  {
1202  m_ZimFwd[0][i] = sqrt(Fwdmu2[i] / Fwdeps3[i]);
1203  m_ZimBwd[0][i] = sqrt(Bwdmu2[i] / Bwdeps3[i]);
1204 
1205  m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1206  m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1207 
1208  m_ZimFwd[1][i] = sqrt(Fwdmu1[i] / Fwdeps3[i]);
1209  m_ZimBwd[1][i] = sqrt(Bwdmu1[i] / Bwdeps3[i]);
1210 
1211  m_YimFwd[1][i] = 1.0 / m_ZimFwd[1][i];
1212  m_YimBwd[1][i] = 1.0 / m_ZimBwd[1][i];
1213  }
1214  }
1215  break; // eTransMagnetic
1216 
1217  case eTransElectric:
1218  {
1219  Array<OneD, NekDouble> Fwdeps1(nTraceNumPoints);
1220  Array<OneD, NekDouble> Bwdeps1(nTraceNumPoints);
1221  Array<OneD, NekDouble> Fwdeps2(nTraceNumPoints);
1222  Array<OneD, NekDouble> Bwdeps2(nTraceNumPoints);
1223  Array<OneD, NekDouble> Fwdmu3(nTraceNumPoints);
1224  Array<OneD, NekDouble> Bwdmu3(nTraceNumPoints);
1225 
1226  m_fields[0]->GetFwdBwdTracePhys(epsvec[0], Fwdeps1,
1227  Bwdeps1);
1228  m_fields[0]->GetFwdBwdTracePhys(epsvec[1], Fwdeps2,
1229  Bwdeps2);
1230  m_fields[0]->GetFwdBwdTracePhys(muvec[2], Fwdmu3, Bwdmu3);
1231 
1232  CopyBoundaryTrace(Fwdeps1, Bwdeps1, eFwdEQBwd, 0);
1233  CopyBoundaryTrace(Fwdeps2, Bwdeps2, eFwdEQBwd, 1);
1234  CopyBoundaryTrace(Fwdmu3, Bwdmu3, eFwdEQBwd, 2);
1235 
1236  for (int i = 0; i < nTraceNumPoints; ++i)
1237  {
1238  m_ZimFwd[0][i] = sqrt(Fwdmu3[i] / Fwdeps2[i]);
1239  m_ZimBwd[0][i] = sqrt(Bwdmu3[i] / Bwdeps2[i]);
1240 
1241  m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1242  m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1243 
1244  m_ZimFwd[1][i] = sqrt(Fwdmu3[i] / Fwdeps1[i]);
1245  m_ZimBwd[1][i] = sqrt(Bwdmu3[i] / Bwdeps1[i]);
1246 
1247  m_YimFwd[1][i] = 1.0 / m_ZimFwd[1][i];
1248  m_YimBwd[1][i] = 1.0 / m_ZimBwd[1][i];
1249  }
1250  }
1251  break; // eTransELectric
1252 
1253  default:
1254  break;
1255  } // PolType
1256 
1257  std::cout << "*** ZimFwd0 = [ "
1258  << Vmath::Vmin(nTraceNumPoints, m_ZimFwd[0], 1) << " , "
1259  << Vmath::Vmax(nTraceNumPoints, m_ZimFwd[0], 1)
1260  << " ], ZimBwd0 = [ "
1261  << Vmath::Vmin(nTraceNumPoints, m_ZimBwd[0], 1) << " , "
1262  << Vmath::Vmax(nTraceNumPoints, m_ZimBwd[0], 1) << " ] "
1263  << std::endl;
1264  std::cout << "*** ZimFwd1 = [ "
1265  << Vmath::Vmin(nTraceNumPoints, m_ZimFwd[1], 1) << " , "
1266  << Vmath::Vmax(nTraceNumPoints, m_ZimFwd[1], 1)
1267  << " ], ZimBwd1 = [ "
1268  << Vmath::Vmin(nTraceNumPoints, m_ZimBwd[1], 1) << " , "
1269  << Vmath::Vmax(nTraceNumPoints, m_ZimBwd[1], 1) << " ] "
1270  << std::endl;
1271  }
1272  break; // eMaxwell2D
1273 
1274  default:
1275  break;
1276  } // TestMaxwellType
1277 }
Array< OneD, Array< OneD, NekDouble > > m_ZimFwd
Definition: MMFSystem.h:207
Array< OneD, Array< OneD, NekDouble > > m_YimBwd
Definition: MMFSystem.h:210
Array< OneD, Array< OneD, NekDouble > > m_ZimBwd
Definition: MMFSystem.h:208
Array< OneD, Array< OneD, NekDouble > > m_YimFwd
Definition: MMFSystem.h:209
TestMaxwellType m_TestMaxwellType
Definition: MMFSystem.h:156
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:1050
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:945

References CopyBoundaryTrace(), Nektar::SolverUtils::eELF2DSurface, Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::eMaxwell1D, Nektar::SolverUtils::eMaxwellSphere, Nektar::SolverUtils::eScatField1D, Nektar::SolverUtils::eScatField2D, Nektar::SolverUtils::eTestMaxwell2DPEC, Nektar::SolverUtils::eTestMaxwell2DPECAVGFLUX, Nektar::SolverUtils::eTestMaxwell2DPMC, Nektar::SolverUtils::eTotField2D, Nektar::SolverUtils::eTransElectric, Nektar::SolverUtils::eTransMagnetic, Nektar::SolverUtils::EquationSystem::GetTraceNpoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_PolType, m_shapedim, m_TestMaxwellType, m_YimBwd, m_YimFwd, m_ZimBwd, m_ZimFwd, RootMeanSquare(), tinysimd::sqrt(), Vmath::Vmax(), and Vmath::Vmin().

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

◆ CopyBoundaryTrace()

void Nektar::SolverUtils::MMFSystem::CopyBoundaryTrace ( const Array< OneD, const NekDouble > &  Fwd,
Array< OneD, NekDouble > &  Bwd,
const BoundaryCopyType  BDCopyType,
const int  var = 0,
const std::string  btype = "NoUserDefined" 
)

Definition at line 838 of file MMFSystem.cpp.

842 {
843  int id1, id2, npts, nptselem, cnt = 0, bdrycnt = 0;
844  Array<OneD, NekDouble> Dirichlet, x0, x1, x2;
845 
846  // loop over Boundary Regions
847  for (int n = 0; n < m_fields[var]->GetBndConditions().size(); ++n)
848  {
849  nptselem = m_fields[var]->GetBndCondExpansions()[n]->GetNpoints();
850 
851  Dirichlet = Array<OneD, NekDouble>(nptselem);
852  x0 = Array<OneD, NekDouble>(nptselem);
853  x1 = Array<OneD, NekDouble>(nptselem);
854  x2 = Array<OneD, NekDouble>(nptselem);
855 
856  if (BDCopyType == eDirichlet)
857  {
858  m_fields[var]->GetBndCondExpansions()[n]->GetCoords(x0, x1, x2);
860  m_session->GetFunction("BoundaryConditions", 0);
861  ifunc->Evaluate(x0, x1, x2, 0.0, Dirichlet);
862  }
863 
864  for (int e = 0;
865  e < m_fields[var]->GetBndCondExpansions()[n]->GetExpSize(); ++e)
866  {
867  npts = m_fields[var]
868  ->GetBndCondExpansions()[n]
869  ->GetExp(e)
870  ->GetNumPoints(0);
871  id1 = m_fields[var]->GetBndCondExpansions()[n]->GetPhys_Offset(e);
872  id2 = m_fields[var]->GetTrace()->GetPhys_Offset(
873  m_fields[var]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt +
874  e));
875 
876  if (m_fields[var]->GetBndConditions()[n]->GetUserDefined() ==
877  BDtype ||
878  BDtype == "NoUserDefined")
879  {
880  switch (BDCopyType)
881  {
882  case eDirichlet:
883  {
884  Vmath::Vcopy(npts, &Dirichlet[id1], 1, &Bwd[id2], 1);
885  bdrycnt++;
886  }
887  break;
888 
889  case eFwdEQBwd:
890  {
891  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
892  bdrycnt++;
893  }
894  break;
895 
896  case eFwdEQNegBwd:
897  {
898  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
899  Vmath::Neg(npts, &Bwd[id2], 1);
900  bdrycnt++;
901  }
902  break;
903 
904  default:
905  break;
906  }
907  }
908  }
909 
910  cnt += m_fields[var]->GetBndCondExpansions()[n]->GetExpSize();
911  }
912 }
SOLVER_UTILS_EXPORT int GetExpSize()
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:129

References Nektar::SolverUtils::eDirichlet, Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::eFwdEQNegBwd, Nektar::SolverUtils::EquationSystem::GetExpSize(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_session, Vmath::Neg(), and Vmath::Vcopy().

Referenced by AverageMaxwellFlux1D(), ComputeMFtrace(), ComputencdotMF(), ComputeZimYim(), LaxFriedrichMaxwellFlux1D(), NumericalMaxwellFluxTE(), NumericalMaxwellFluxTM(), Nektar::MMFSWE::NumericalSWEFlux(), and UpwindMaxwellFlux1D().

◆ DeriveCrossProductMF()

void Nektar::SolverUtils::MMFSystem::DeriveCrossProductMF ( Array< OneD, Array< OneD, NekDouble >> &  CrossProductMF)
protected

Definition at line 571 of file MMFSystem.cpp.

573 {
574  int MFdim = 3;
575  int nq = GetTotPoints();
576 
577  CrossProductMF = Array<OneD, Array<OneD, NekDouble>>(MFdim);
578 
579  Array<OneD, Array<OneD, NekDouble>> MF1tmp(m_spacedim);
580  Array<OneD, Array<OneD, NekDouble>> MF2tmp(m_spacedim);
581  Array<OneD, Array<OneD, NekDouble>> MF3tmp(m_spacedim);
582  Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtmpCurl(MFdim);
583  for (int j = 0; j < MFdim; ++j)
584  {
585  CrossProductMF[j] = Array<OneD, NekDouble>(nq * m_spacedim);
586  MFtmpCurl[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
587  for (int k = 0; k < m_spacedim; ++k)
588  {
589  MFtmpCurl[j][k] = Array<OneD, NekDouble>(nq);
590  }
591  }
592 
593  for (int k = 0; k < m_spacedim; ++k)
594  {
595  MF1tmp[k] = Array<OneD, NekDouble>(nq);
596  MF2tmp[k] = Array<OneD, NekDouble>(nq);
597  MF3tmp[k] = Array<OneD, NekDouble>(nq);
598 
599  Vmath::Vcopy(nq, &m_movingframes[0][k * nq], 1, &MF1tmp[k][0], 1);
600  Vmath::Vcopy(nq, &m_movingframes[1][k * nq], 1, &MF2tmp[k][0], 1);
601  Vmath::Vcopy(nq, &m_movingframes[2][k * nq], 1, &MF3tmp[k][0], 1);
602  }
603 
604  VectorCrossProd(MF3tmp, MF1tmp, MFtmpCurl[0]);
605  VectorCrossProd(MF2tmp, MF3tmp, MFtmpCurl[1]);
606  VectorCrossProd(MF1tmp, MF2tmp, MFtmpCurl[2]);
607 
608  for (int j = 0; j < MFdim; ++j)
609  {
610  for (int k = 0; k < m_spacedim; ++k)
611  {
612  Vmath::Vcopy(nq, &MFtmpCurl[j][k][0], 1, &CrossProductMF[j][k * nq],
613  1);
614  }
615  }
616 }

References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_movingframes, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vcopy(), and VectorCrossProd().

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

◆ GetIncidentField()

Array< OneD, NekDouble > Nektar::SolverUtils::MMFSystem::GetIncidentField ( const int  var,
const NekDouble  time 
)
protected

Definition at line 2011 of file MMFSystem.cpp.

2013 {
2014  int nq = m_fields[0]->GetNpoints();
2015 
2016  Array<OneD, NekDouble> x0(nq);
2017  Array<OneD, NekDouble> x1(nq);
2018  Array<OneD, NekDouble> x2(nq);
2019 
2020  m_fields[0]->GetCoords(x0, x1, x2);
2021 
2022  // GetSmoothFactor such that wave propages from the left to the object.
2023  // a = 0.1, ta = 1, f = 1.0./(1.0 + exp( -0.5.*(time-ta)/a ));
2024  Array<OneD, NekDouble> SmoothFactor(nq, 1.0);
2025  switch (m_SmoothFactor)
2026  {
2027  case 0:
2028  {
2029  for (int i = 0; i < nq; i++)
2030  {
2031  SmoothFactor[i] = 1.0 / (1.0 + exp(-1.0 * (time - 1.0) / 0.1));
2032  }
2033  }
2034  break;
2035 
2036  case 1:
2037  {
2038  NekDouble xmin = Vmath::Vmin(nq, x0, 1);
2039  NekDouble xp;
2040  for (int i = 0; i < nq; i++)
2041  {
2042  xp = x0[i] - xmin - time - m_SFinit;
2043  if (xp > 0.0)
2044  {
2045  SmoothFactor[i] =
2046  2.0 / (1.0 + exp(0.5 * (sqrt(xp * xp) - 0.1)));
2047  }
2048 
2049  else
2050  {
2051  SmoothFactor[i] = 1.0;
2052  }
2053  }
2054  }
2055  break;
2056 
2057  default:
2058  break;
2059  }
2060 
2061  // Generate a factor for smoothly increasing wave
2062  Array<OneD, NekDouble> F1(nq);
2063  Array<OneD, NekDouble> F2(nq);
2064  Array<OneD, NekDouble> F3(nq);
2065 
2066  Array<OneD, NekDouble> dF1dt(nq);
2067  Array<OneD, NekDouble> dF2dt(nq);
2068  Array<OneD, NekDouble> dF3dt(nq);
2069 
2070  Array<OneD, NekDouble> F1int(nq);
2071  Array<OneD, NekDouble> F2int(nq);
2072  Array<OneD, NekDouble> F3int(nq);
2073 
2074  Array<OneD, NekDouble> outarray(nq);
2075 
2076  switch (m_IncType)
2077  {
2078  case ePlaneWave:
2079  {
2080  NekDouble cs, sn;
2081  NekDouble e1y, e2y;
2082  switch (m_PolType)
2083  {
2084  case eTransMagnetic:
2085  {
2086  // H = 0 \hat{x} - e^{ikx} \hat{y}
2087  // E = e^{ikx} \hat{z}
2088  // outarray1 = Hr1inc
2089  // outarray2 = Hr2inc
2090  // outarray3 = Ezrinc
2091  for (int i = 0; i < nq; i++)
2092  {
2093  e1y = m_movingframes[0][nq + i];
2094  e2y = m_movingframes[1][nq + i];
2095 
2096  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2097  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2098 
2099  F1[i] = -1.0 * cs * e1y;
2100  F2[i] = -1.0 * cs * e2y;
2101  F3[i] = cs;
2102 
2103  dF1dt[i] = -1.0 * m_Incfreq * sn * e1y;
2104  dF2dt[i] = -1.0 * m_Incfreq * sn * e2y;
2105  dF3dt[i] = 1.0 * m_Incfreq * sn;
2106 
2107  F1int[i] = (1.0 / m_Incfreq) * sn * e1y;
2108  F2int[i] = (1.0 / m_Incfreq) * sn * e2y;
2109  F3int[i] = (-1.0 / m_Incfreq) * sn;
2110  }
2111  }
2112  break;
2113 
2114  case eTransElectric:
2115  {
2116  // E = 0 \hat{x} + e^{ikx} \hat{y}
2117  // H = e^{ikx} \hat{z}
2118  // outarray1 = Er1inc
2119  // outarray2 = Er2inc
2120  // outarray3 = Hzrinc
2121  // outarray1 = Ei1inc
2122  // outarray2 = Ei2inc
2123  // outarray3 = Hziinc
2124  for (int i = 0; i < nq; i++)
2125  {
2126  e1y = m_movingframes[0][nq + i];
2127  e2y = m_movingframes[1][nq + i];
2128 
2129  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2130  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2131 
2132  F1[i] = cs * e1y;
2133  F2[i] = cs * e2y;
2134  F3[i] = cs;
2135 
2136  dF1dt[i] = m_Incfreq * sn * e1y;
2137  dF2dt[i] = m_Incfreq * sn * e2y;
2138  dF3dt[i] = m_Incfreq * sn;
2139 
2140  F1int[i] = (-1.0 / m_Incfreq) * sn * e1y;
2141  F2int[i] = (-1.0 / m_Incfreq) * sn * e2y;
2142  F3int[i] = (-1.0 / m_Incfreq) * sn;
2143  }
2144  }
2145  break;
2146 
2147  default:
2148  break;
2149  }
2150  }
2151  break;
2152 
2153  case ePlaneWaveImag:
2154  {
2155  NekDouble cs, sn;
2156  NekDouble e1y, e2y;
2157  switch (m_PolType)
2158  {
2159  case eTransMagnetic:
2160  {
2161  // H = 0 \hat{x} - e^{ikx} \hat{y}
2162  // E = e^{ikx} \hat{z}
2163  // outarray1 = Hr1inc
2164  // outarray2 = Hr2inc
2165  // outarray3 = Ezrinc
2166  for (int i = 0; i < nq; i++)
2167  {
2168  e1y = m_movingframes[0][nq + i];
2169  e2y = m_movingframes[1][nq + i];
2170 
2171  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2172  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2173 
2174  F1[i] = -1.0 * sn * e1y;
2175  F2[i] = -1.0 * sn * e2y;
2176  F3[i] = sn;
2177 
2178  dF1dt[i] = m_Incfreq * cs * e1y;
2179  dF2dt[i] = m_Incfreq * cs * e2y;
2180  dF3dt[i] = -1.0 * m_Incfreq * cs;
2181 
2182  F1int[i] = (-1.0 / m_Incfreq) * cs * e1y;
2183  F2int[i] = (-1.0 / m_Incfreq) * cs * e2y;
2184  F3int[i] = (1.0 / m_Incfreq) * cs;
2185  }
2186  }
2187  break;
2188 
2189  case eTransElectric:
2190  {
2191  // E = 0 \hat{x} + e^{ikx} \hat{y}
2192  // H = e^{ikx} \hat{z}
2193  // outarray1 = Er1inc
2194  // outarray2 = Er2inc
2195  // outarray3 = Hzrinc
2196  // outarray1 = Ei1inc
2197  // outarray2 = Ei2inc
2198  // outarray3 = Hziinc
2199  for (int i = 0; i < nq; i++)
2200  {
2201  e1y = m_movingframes[0][nq + i];
2202  e2y = m_movingframes[1][nq + i];
2203 
2204  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2205  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2206 
2207  F1[i] = sn * e1y;
2208  F2[i] = sn * e2y;
2209  F3[i] = sn;
2210 
2211  dF1dt[i] = -1.0 * m_Incfreq * cs * e1y;
2212  dF2dt[i] = -1.0 * m_Incfreq * cs * e2y;
2213  dF3dt[i] = -1.0 * m_Incfreq * cs;
2214 
2215  F1int[i] = (1.0 / m_Incfreq) * cs * e1y;
2216  F2int[i] = (1.0 / m_Incfreq) * cs * e2y;
2217  F3int[i] = (1.0 / m_Incfreq) * cs;
2218  }
2219  }
2220  break;
2221 
2222  default:
2223  break;
2224  }
2225  }
2226  break;
2227 
2228  default:
2229  break;
2230  }
2231 
2232  switch (var)
2233  {
2234  case 0:
2235  {
2236  outarray = F1;
2237  }
2238  break;
2239 
2240  case 1:
2241  {
2242  outarray = F2;
2243  }
2244  break;
2245 
2246  case 2:
2247  {
2248  outarray = F3;
2249  }
2250  break;
2251 
2252  case 10:
2253  {
2254  outarray = dF1dt;
2255  }
2256  break;
2257 
2258  case 11:
2259  {
2260  outarray = dF2dt;
2261  }
2262  break;
2263 
2264  case 12:
2265  {
2266  outarray = dF3dt;
2267  }
2268  break;
2269 
2270  case 20:
2271  {
2272  outarray = F1int;
2273  }
2274  break;
2275 
2276  case 21:
2277  {
2278  outarray = F2int;
2279  }
2280  break;
2281 
2282  case 22:
2283  {
2284  outarray = F3int;
2285  }
2286  break;
2287 
2288  default:
2289  {
2290  Vmath::Zero(nq, outarray, 1);
2291  }
2292  break;
2293  }
2294 
2295  return outarray;
2296 }
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492

References Nektar::SolverUtils::ePlaneWave, Nektar::SolverUtils::ePlaneWaveImag, Nektar::SolverUtils::eTransElectric, Nektar::SolverUtils::eTransMagnetic, Nektar::SolverUtils::EquationSystem::m_fields, m_Incfreq, m_IncType, m_movingframes, m_PolType, m_SFinit, m_SmoothFactor, tinysimd::sqrt(), Vmath::Vmin(), and Vmath::Zero().

Referenced by Nektar::MMFMaxwell::Checkpoint_TotalFieldOutput(), Nektar::MMFMaxwell::Checkpoint_TotPlotOutput(), Nektar::MMFMaxwell::ComputeSurfaceCurrent(), Nektar::MMFMaxwell::DoOdeRhs(), NumericalMaxwellFluxTE(), NumericalMaxwellFluxTM(), and Nektar::MMFMaxwell::v_DoSolve().

◆ GetMaxwellFlux1D()

void Nektar::SolverUtils::MMFSystem::GetMaxwellFlux1D ( const int  var,
const Array< OneD, const Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  flux 
)
protected

Definition at line 1639 of file MMFSystem.cpp.

1642 {
1643  int nq = m_fields[0]->GetTotPoints();
1644 
1645  switch (var)
1646  {
1647  case 0:
1648  {
1649  // H in flux 0
1650  Vmath::Vcopy(nq, physfield[1], 1, flux[0], 1);
1651 
1652  // E in flux 1
1653  Vmath::Zero(nq, flux[1], 1);
1654  }
1655  break;
1656 
1657  case 1:
1658  {
1659  // E in flux 0
1660  Vmath::Vcopy(nq, physfield[0], 1, flux[0], 1);
1661 
1662  // H in flux 1
1663  Vmath::Zero(nq, flux[1], 1);
1664  }
1665  break;
1666  //----------------------------------------------------
1667 
1668  default:
1669  break;
1670  }
1671 }

References Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Vcopy(), and Vmath::Zero().

Referenced by GetMaxwellFluxVector().

◆ GetMaxwellFlux2D()

void Nektar::SolverUtils::MMFSystem::GetMaxwellFlux2D ( const int  var,
const Array< OneD, const Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  flux 
)
protected

Definition at line 1673 of file MMFSystem.cpp.

1676 {
1677  int nq = m_fields[0]->GetTotPoints();
1678 
1679  NekDouble sign = 1.0;
1680  switch (m_PolType)
1681  {
1682  // TransMagnetic
1683  case 0:
1684  {
1685  sign = -1.0;
1686  }
1687  break;
1688 
1689  // TransElectric
1690  case 1:
1691  {
1692  sign = 1.0;
1693  }
1694  break;
1695 
1696  default:
1697  break;
1698  }
1699 
1700  switch (var)
1701  {
1702  case 0:
1703  {
1704  // -Ez in flux 1
1705  Vmath::Smul(nq, sign, physfield[2], 1, flux[0], 1);
1706  Vmath::Zero(nq, flux[1], 1);
1707  }
1708  break;
1709 
1710  case 1:
1711  {
1712  // Ez in flux 0
1713  Vmath::Zero(nq, flux[0], 1);
1714  Vmath::Smul(nq, -sign, physfield[2], 1, flux[1], 1);
1715  }
1716  break;
1717 
1718  case 2:
1719  {
1720  Vmath::Smul(nq, sign, physfield[0], 1, flux[0], 1);
1721  Vmath::Smul(nq, -sign, physfield[1], 1, flux[1], 1);
1722  }
1723  break;
1724 
1725  default:
1726  ASSERTL0(false, "GetFluxVector2D: illegal vector index");
1727  }
1728 }
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:49
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248

References ASSERTL0, Nektar::SolverUtils::EquationSystem::m_fields, m_PolType, sign, Vmath::Smul(), and Vmath::Zero().

Referenced by GetMaxwellFluxVector().

◆ GetMaxwellFluxVector()

void Nektar::SolverUtils::MMFSystem::GetMaxwellFluxVector ( const int  var,
const Array< OneD, const Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  flux 
)
protected

Definition at line 1609 of file MMFSystem.cpp.

1612 {
1613  switch (m_TestMaxwellType)
1614  {
1615  case eMaxwell1D:
1616  case eScatField1D:
1617  {
1618  GetMaxwellFlux1D(var, physfield, flux);
1619  }
1620  break;
1621 
1622  case eTestMaxwell2DPEC:
1624  case eTestMaxwell2DPMC:
1625  case eScatField2D:
1626  case eTotField2D:
1627  case eMaxwellSphere:
1628  case eELF2DSurface:
1629  {
1630  GetMaxwellFlux2D(var, physfield, flux);
1631  }
1632  break;
1633 
1634  default:
1635  break;
1636  }
1637 }
SOLVER_UTILS_EXPORT void GetMaxwellFlux2D(const int var, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
Definition: MMFSystem.cpp:1673
SOLVER_UTILS_EXPORT void GetMaxwellFlux1D(const int var, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
Definition: MMFSystem.cpp:1639

References Nektar::SolverUtils::eELF2DSurface, Nektar::SolverUtils::eMaxwell1D, Nektar::SolverUtils::eMaxwellSphere, Nektar::SolverUtils::eScatField1D, Nektar::SolverUtils::eScatField2D, Nektar::SolverUtils::eTestMaxwell2DPEC, Nektar::SolverUtils::eTestMaxwell2DPECAVGFLUX, Nektar::SolverUtils::eTestMaxwell2DPMC, Nektar::SolverUtils::eTotField2D, GetMaxwellFlux1D(), GetMaxwellFlux2D(), and m_TestMaxwellType.

Referenced by Nektar::MMFMaxwell::AddGreenDerivCompensate(), and Nektar::MMFMaxwell::WeakDGMaxwellDirDeriv().

◆ GramSchumitz()

void Nektar::SolverUtils::MMFSystem::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 
)
protected

Definition at line 2405 of file MMFSystem.cpp.

2409 {
2410 
2411  int nq = v1[0].size();
2412  Array<OneD, NekDouble> tmp(nq, 0.0);
2413  Array<OneD, NekDouble> mag(nq, 0.0);
2414 
2415  for (int i = 0; i < m_spacedim; ++i)
2416  {
2417  // u2 = v2 - < u1 , v2 > ( u1 / < u1, u1 > )
2418  Vmath::Vvtvp(nq, &v1[i][0], 1, &v2[i][0], 1, &tmp[0], 1, &tmp[0], 1);
2419  Vmath::Vvtvp(nq, &v1[i][0], 1, &v1[i][0], 1, &mag[0], 1, &mag[0], 1);
2420  }
2421  Vmath::Vdiv(nq, &tmp[0], 1, &mag[0], 1, &tmp[0], 1);
2422  Vmath::Neg(nq, &tmp[0], 1);
2423 
2424  // outarray = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
2425 
2426  // u2 = v2 - < u1 , v2 > ( u1 / < u1, u1 > )
2427  for (int i = 0; i < m_spacedim; ++i)
2428  {
2429  outarray[i] = Array<OneD, NekDouble>(nq, 0.0);
2430  Vmath::Vvtvp(nq, &tmp[0], 1, &v1[i][0], 1, &v2[i][0], 1,
2431  &outarray[i][0], 1);
2432  }
2433 
2434  if (KeepTheMagnitude)
2435  {
2436  Array<OneD, NekDouble> magorig(nq, 0.0);
2437  Array<OneD, NekDouble> magnew(nq, 0.0);
2438 
2439  for (int i = 0; i < m_spacedim; ++i)
2440  {
2441  Vmath::Vmul(nq, &v2[0][0], 1, &v2[0][0], 1, &magorig[0], 1);
2442  Vmath::Vvtvp(nq, &v2[1][0], 1, &v2[1][0], 1, &magorig[0], 1,
2443  &magorig[0], 1);
2444  Vmath::Vvtvp(nq, &v2[2][0], 1, &v2[2][0], 1, &magorig[0], 1,
2445  &magorig[0], 1);
2446 
2447  Vmath::Vmul(nq, &outarray[0][0], 1, &outarray[0][0], 1, &magnew[0],
2448  1);
2449  Vmath::Vvtvp(nq, &outarray[1][0], 1, &outarray[1][0], 1, &magnew[0],
2450  1, &magnew[0], 1);
2451  Vmath::Vvtvp(nq, &outarray[2][0], 1, &outarray[2][0], 1, &magnew[0],
2452  1, &magnew[0], 1);
2453  }
2454 
2455  for (int i = 0; i < m_spacedim; ++i)
2456  {
2457  for (int j = 0; j < nq; ++j)
2458  {
2459  if (fabs(magnew[j]) > 0.000000001)
2460  {
2461  outarray[i][j] =
2462  outarray[i][j] * sqrt(magorig[j] / magnew[j]);
2463  }
2464  }
2465  }
2466  }
2467 }
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:284

References Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Neg(), tinysimd::sqrt(), Vmath::Vdiv(), Vmath::Vmul(), and Vmath::Vvtvp().

Referenced by Nektar::MMFAdvection::EvaluateAdvectionVelocity().

◆ LaxFriedrichMaxwellFlux1D()

void Nektar::SolverUtils::MMFSystem::LaxFriedrichMaxwellFlux1D ( Array< OneD, Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  numfluxFwd,
Array< OneD, Array< OneD, NekDouble >> &  numfluxBwd 
)
protected

Definition at line 1515 of file MMFSystem.cpp.

1519 {
1520  int i;
1521  int nTraceNumPoints = GetTraceTotPoints();
1522  int nvar = 2;
1523 
1524  // get temporary arrays
1525  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1526  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1527 
1528  for (i = 0; i < nvar; ++i)
1529  {
1530  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1531  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1532  }
1533 
1534  // get the physical values at the trace from the dependent variables
1535  for (i = 0; i < nvar; ++i)
1536  {
1537  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1538  }
1539 
1540  // E = 0 at the boundaries
1541  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1542 
1543  // d H / d n = 0 at the boundaries
1544  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1545 
1546  Array<OneD, NekDouble> dE(nTraceNumPoints);
1547  Array<OneD, NekDouble> dH(nTraceNumPoints);
1548 
1549  Vmath::Vsub(nTraceNumPoints, &Fwd[0][0], 1, &Bwd[0][0], 1, &dE[0], 1);
1550  Vmath::Vsub(nTraceNumPoints, &Fwd[1][0], 1, &Bwd[1][0], 1, &dH[0], 1);
1551 
1552  NekDouble nx;
1553  for (i = 0; i < nTraceNumPoints; ++i)
1554  {
1555  nx = m_traceNormals[0][i];
1556  numfluxFwd[0][i] =
1557  0.5 * nx * ((Fwd[1][i] + Bwd[1][i]) - nx * m_YimFwd[0][i] * dE[i]);
1558  numfluxFwd[1][i] =
1559  0.5 * nx * ((Fwd[0][i] + Bwd[0][i]) - nx * m_ZimFwd[0][i] * dH[i]);
1560 
1561  numfluxBwd[0][i] =
1562  0.5 * nx * ((Fwd[1][i] + Bwd[1][i]) - nx * m_YimFwd[0][i] * dE[i]);
1563  numfluxBwd[1][i] =
1564  0.5 * nx * ((Fwd[0][i] + Bwd[0][i]) - nx * m_ZimFwd[0][i] * dH[i]);
1565  }
1566 }

References CopyBoundaryTrace(), Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::eFwdEQNegBwd, Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_traceNormals, m_YimFwd, m_ZimFwd, and Vmath::Vsub().

Referenced by NumericalMaxwellFlux().

◆ MMFInitObject()

void Nektar::SolverUtils::MMFSystem::MMFInitObject ( const Array< OneD, const Array< OneD, NekDouble >> &  Anisotropy,
const int  TangentXelem = -1 
)

Definition at line 53 of file MMFSystem.cpp.

56 {
57  m_pi = 3.14159265358979323846;
58  m_shapedim = m_fields[0]->GetShapeDimension();
59 
60  ASSERTL0(m_spacedim == 3, "Space Dimension should be 3.");
61 
62  // Factor for Numerical Flux
63  m_session->LoadParameter("alpha", m_alpha, 1.0);
64 
65  // Factor for Numerical Flux
66  m_session->LoadParameter("Incfreq", m_Incfreq, 1.0);
67 
68  // SmoothFactor
69  m_session->LoadParameter("SmoothFactor", m_SmoothFactor, 1);
70  m_session->LoadParameter("SFinit", m_SFinit, 0.0);
71 
72  // Define SurfaceType
73  if (m_session->DefinesSolverInfo("SURFACETYPE"))
74  {
75  std::string SurfaceTypeStr = m_session->GetSolverInfo("SURFACETYPE");
76  int i;
77  for (i = 0; i < (int)SIZE_SurfaceType; ++i)
78  {
79  if (boost::iequals(SurfaceTypeMap[i], SurfaceTypeStr))
80  {
82  break;
83  }
84  }
85  }
86  else
87  {
89  }
90 
91  // if discontinuous Galerkin determine numerical flux to use
92  for (int i = 0; i < (int)SIZE_UpwindType; ++i)
93  {
94  bool match;
95  m_session->MatchSolverInfo("UPWINDTYPE", UpwindTypeMap[i], match,
96  false);
97  if (match)
98  {
100  break;
101  }
102  }
103 
104  // SetUpMovingFrames: To generate m_movingframes
105  SetUpMovingFrames(Anisotropy, TangentXelem);
106 
107  // Derive movingfrmaes at interfaces
108  // Get: m_MFtraceFwd and m_MFtraceBwd
109  ComputeMFtrace();
110 
111  // Check Movingframes and surfraceNormal
112  // Get: m_ncdotMFFwd,m_ncdotMFBwd,m_nperpcdotMFFwd,m_nperpcdotMFBwd
113  ComputencdotMF();
114 
115  // Compute nabla cdot m_movingframes
116  // Get: m_DivMF, m_CurlMF
118 }
void SetUpMovingFrames(const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem)
Definition: MMFSystem.cpp:120
SOLVER_UTILS_EXPORT void ComputeDivCurlMF()
Definition: MMFSystem.cpp:445
SOLVER_UTILS_EXPORT void ComputeMFtrace()
Definition: MMFSystem.cpp:522
SOLVER_UTILS_EXPORT void ComputencdotMF()
Definition: MMFSystem.cpp:320
const char *const UpwindTypeMap[]
Definition: MMFSystem.h:89
@ SIZE_UpwindType
Length of enum list.
Definition: MMFSystem.h:86
const char *const SurfaceTypeMap[]
Definition: MMFSystem.h:57

References ASSERTL0, ComputeDivCurlMF(), ComputeMFtrace(), ComputencdotMF(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_Incfreq, m_pi, Nektar::SolverUtils::EquationSystem::m_session, m_SFinit, m_shapedim, m_SmoothFactor, Nektar::SolverUtils::EquationSystem::m_spacedim, m_surfaceType, m_upwindType, SetUpMovingFrames(), Nektar::SolverUtils::SIZE_SurfaceType, Nektar::SolverUtils::SIZE_UpwindType, Nektar::SolverUtils::SurfaceTypeMap, and Nektar::SolverUtils::UpwindTypeMap.

Referenced by Nektar::MMFDiffusion::v_InitObject(), Nektar::MMFAdvection::v_InitObject(), Nektar::MMFMaxwell::v_InitObject(), and Nektar::MMFSWE::v_InitObject().

◆ NumericalMaxwellFlux()

void Nektar::SolverUtils::MMFSystem::NumericalMaxwellFlux ( Array< OneD, Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  numfluxFwd,
Array< OneD, Array< OneD, NekDouble >> &  numfluxBwd,
const NekDouble  time = 0.0 
)
protected

Definition at line 1730 of file MMFSystem.cpp.

1734 {
1735 
1736  switch (m_TestMaxwellType)
1737  {
1738  case eMaxwell1D:
1739  case eScatField1D:
1740  {
1741  switch (m_upwindType)
1742  {
1743  case eAverage:
1744  {
1745  AverageMaxwellFlux1D(physfield, numfluxFwd, numfluxBwd);
1746  }
1747  break;
1748 
1749  case eLaxFriedrich:
1750  {
1751  LaxFriedrichMaxwellFlux1D(physfield, numfluxFwd,
1752  numfluxBwd);
1753  }
1754  break;
1755 
1756  case eUpwind:
1757  {
1758  UpwindMaxwellFlux1D(physfield, numfluxFwd, numfluxBwd);
1759  }
1760  break;
1761 
1762  default:
1763  {
1764  ASSERTL0(false,
1765  "populate switch statement for upwind flux");
1766  }
1767  break; // upwindType
1768  }
1769  }
1770  break; // eMaxwell1D
1771 
1772  case eTestMaxwell2DPEC:
1774  case eTestMaxwell2DPMC:
1775  case eScatField2D:
1776  case eTotField2D:
1777  case eMaxwellSphere:
1778  case eELF2DSurface:
1779  {
1780  switch (m_PolType)
1781  {
1782  case eTransMagnetic:
1783  {
1784  NumericalMaxwellFluxTM(physfield, numfluxFwd, numfluxBwd,
1785  time);
1786  }
1787  break;
1788 
1789  case eTransElectric:
1790  {
1791  NumericalMaxwellFluxTE(physfield, numfluxFwd, numfluxBwd,
1792  time);
1793  }
1794  break;
1795 
1796  default:
1797  break;
1798  }
1799  }
1800  break; // eMaxwell2D
1801 
1802  default:
1803  break;
1804  } // m_TestMaxwellType
1805 }
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)
Definition: MMFSystem.cpp:1807
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)
Definition: MMFSystem.cpp:1908
SOLVER_UTILS_EXPORT void AverageMaxwellFlux1D(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
Definition: MMFSystem.cpp:1568
SOLVER_UTILS_EXPORT void LaxFriedrichMaxwellFlux1D(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
Definition: MMFSystem.cpp:1515
SOLVER_UTILS_EXPORT void UpwindMaxwellFlux1D(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
Definition: MMFSystem.cpp:1461
@ eAverage
averaged (or centred) flux
Definition: MMFSystem.h:80
@ eLaxFriedrich
Lax-Friedrich flux.
Definition: MMFSystem.h:81

References ASSERTL0, AverageMaxwellFlux1D(), Nektar::SolverUtils::eAverage, Nektar::SolverUtils::eELF2DSurface, Nektar::SolverUtils::eLaxFriedrich, Nektar::SolverUtils::eMaxwell1D, Nektar::SolverUtils::eMaxwellSphere, Nektar::SolverUtils::eScatField1D, Nektar::SolverUtils::eScatField2D, Nektar::SolverUtils::eTestMaxwell2DPEC, Nektar::SolverUtils::eTestMaxwell2DPECAVGFLUX, Nektar::SolverUtils::eTestMaxwell2DPMC, Nektar::SolverUtils::eTotField2D, Nektar::SolverUtils::eTransElectric, Nektar::SolverUtils::eTransMagnetic, Nektar::SolverUtils::eUpwind, LaxFriedrichMaxwellFlux1D(), m_PolType, m_TestMaxwellType, m_upwindType, NumericalMaxwellFluxTE(), NumericalMaxwellFluxTM(), and UpwindMaxwellFlux1D().

Referenced by Nektar::MMFMaxwell::WeakDGMaxwellDirDeriv().

◆ NumericalMaxwellFluxTE()

void Nektar::SolverUtils::MMFSystem::NumericalMaxwellFluxTE ( Array< OneD, Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  numfluxFwd,
Array< OneD, Array< OneD, NekDouble >> &  numfluxBwd,
const NekDouble  time 
)
protected

Definition at line 1908 of file MMFSystem.cpp.

1913 {
1914  int nq = m_fields[0]->GetNpoints();
1915  int nTraceNumPoints = GetTraceTotPoints();
1916  int nvar = physfield.size();
1917 
1918  // Get temporary arrays
1919  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1920  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1921  for (int i = 0; i < nvar; ++i)
1922  {
1923  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1924  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1925 
1926  // get the physical values at the trace
1927  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1928  }
1929 
1930  // E = 0 at the PEC boundaries:
1931  Array<OneD, NekDouble> IncField(nq, 0.0);
1932  Array<OneD, NekDouble> IncFieldBwd(nTraceNumPoints, 0.0);
1933  Array<OneD, Array<OneD, NekDouble>> IncFieldFwd(m_spacedim);
1934  for (int i = 0; i < m_spacedim; ++i)
1935  {
1936  IncFieldFwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1937 
1938  IncField = GetIncidentField(i, time);
1939  m_fields[0]->GetFwdBwdTracePhys(IncField, IncFieldFwd[i], IncFieldBwd);
1940 
1941  Vmath::Svtvp(nTraceNumPoints, 2.0, &IncFieldFwd[i][0], 1, &Fwd[i][0], 1,
1942  &IncFieldFwd[i][0], 1);
1943  }
1944 
1945  CopyBoundaryTrace(IncFieldFwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1946  "PEC_Forces");
1947  CopyBoundaryTrace(IncFieldFwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1,
1948  "PEC_Forces");
1949  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PEC_Forces");
1950 
1951  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PMC_Forces");
1952  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PMC_Forces");
1953  CopyBoundaryTrace(IncFieldFwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2,
1954  "PMC_Forces");
1955 
1956  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1957  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1, "PEC");
1958  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PEC");
1959 
1960  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PMC");
1961  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PMC");
1962  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2, "PMC");
1963 
1964  Array<OneD, NekDouble> e1Fwd_cdot_ncrossdE(nTraceNumPoints);
1965  Array<OneD, NekDouble> e1Bwd_cdot_ncrossdE(nTraceNumPoints);
1966  Array<OneD, NekDouble> e2Fwd_cdot_ncrossdE(nTraceNumPoints);
1967  Array<OneD, NekDouble> e2Bwd_cdot_ncrossdE(nTraceNumPoints);
1968  Array<OneD, NekDouble> e3Fwd_cdot_dHe3(nTraceNumPoints);
1969  Array<OneD, NekDouble> e3Bwd_cdot_dHe3(nTraceNumPoints);
1970 
1971  // Compute numfluxFwd[dir] = (eFwd^[dir] \cdot n \times e^3) * (ZimFwd *
1972  // HFwd^3 + ZimBwd * HBwd^3 )
1973  // Compute numfluxBwd[dir] = (eBwd^[dir] \cdot n \times e^3) * (ZimFwd *
1974  // HFwd^3 + ZimBwd * HBwd^3 )
1975  ComputeNtimesFz(0, Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], numfluxFwd[0],
1976  numfluxBwd[0]);
1977  ComputeNtimesFz(1, Fwd, Bwd, m_ZimFwd[1], m_ZimBwd[1], numfluxFwd[1],
1978  numfluxBwd[1]);
1979 
1980  // Compute numfluxFwd[2] = eFwd^3 \cdot ( n1e1 \times ( imFwd EFwd + imBwd
1981  // EBwd ) ) / 2 {{Y_i}}
1982  ComputeNtimesF12(Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], m_YimFwd[1],
1983  m_YimBwd[1], numfluxFwd[2], numfluxBwd[2]);
1984 
1985  // Compute e1Fwd_cdot_ncrossdE = eFwd[dir] \cdot \alpha n \times n \times
1986  // [E] / 2 {{ZimFwd}}
1987  ComputeNtimestimesdFz(0, Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0],
1988  e1Fwd_cdot_ncrossdE, e1Bwd_cdot_ncrossdE);
1989  ComputeNtimestimesdFz(1, Fwd, Bwd, m_ZimFwd[1], m_ZimBwd[1],
1990  e2Fwd_cdot_ncrossdE, e2Bwd_cdot_ncrossdE);
1991 
1992  // Compute - \alpha [H3] * ( 1/2{{Yim1}} + 1/2{{Yim2}} )
1993  ComputeNtimestimesdF12(Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], m_YimFwd[1],
1994  m_YimBwd[1], e3Fwd_cdot_dHe3, e3Bwd_cdot_dHe3);
1995 
1996  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[0], 1, e1Fwd_cdot_ncrossdE, 1,
1997  numfluxFwd[0], 1);
1998  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[1], 1, e2Fwd_cdot_ncrossdE, 1,
1999  numfluxFwd[1], 1);
2000  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[2], 1, e3Fwd_cdot_dHe3, 1,
2001  numfluxFwd[2], 1);
2002 
2003  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[0], 1, e1Bwd_cdot_ncrossdE, 1,
2004  numfluxBwd[0], 1);
2005  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[1], 1, e2Bwd_cdot_ncrossdE, 1,
2006  numfluxBwd[1], 1);
2007  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[2], 1, e3Bwd_cdot_dHe3, 1,
2008  numfluxBwd[2], 1);
2009 }
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)
Definition: MMFSystem.cpp:1014
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)
Definition: MMFSystem.cpp:1080
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetIncidentField(const int var, const NekDouble time)
Definition: MMFSystem.cpp:2011
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)
Definition: MMFSystem.cpp:951
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)
Definition: MMFSystem.cpp:915
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:622

References ComputeNtimesF12(), ComputeNtimesFz(), ComputeNtimestimesdF12(), ComputeNtimestimesdFz(), CopyBoundaryTrace(), Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::eFwdEQNegBwd, GetIncidentField(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, m_YimBwd, m_YimFwd, m_ZimBwd, m_ZimFwd, and Vmath::Svtvp().

Referenced by NumericalMaxwellFlux().

◆ NumericalMaxwellFluxTM()

void Nektar::SolverUtils::MMFSystem::NumericalMaxwellFluxTM ( Array< OneD, Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  numfluxFwd,
Array< OneD, Array< OneD, NekDouble >> &  numfluxBwd,
const NekDouble  time 
)
protected

Definition at line 1807 of file MMFSystem.cpp.

1811 {
1812  int nq = m_fields[0]->GetNpoints();
1813  int nTraceNumPoints = GetTraceTotPoints();
1814  int nvar = physfield.size();
1815 
1816  // get temporary arrays
1817  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1818  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1819  for (int i = 0; i < nvar; ++i)
1820  {
1821  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1822  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1823 
1824  // get the physical values at the trace
1825  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1826  }
1827 
1828  // E^|| = 0 at the PEC boundaries vs. H^|| = 0 at PMC boundaries
1829  Array<OneD, NekDouble> IncField(nq, 0.0);
1830  Array<OneD, NekDouble> IncFieldBwd(nTraceNumPoints, 0.0);
1831  Array<OneD, Array<OneD, NekDouble>> IncFieldFwd(m_spacedim);
1832  for (int i = 0; i < m_spacedim; ++i)
1833  {
1834  IncFieldFwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1835 
1836  IncField = GetIncidentField(i, time);
1837  m_fields[0]->GetFwdBwdTracePhys(IncField, IncFieldFwd[i], IncFieldBwd);
1838 
1839  Vmath::Svtvp(nTraceNumPoints, 2.0, &IncFieldFwd[i][0], 1, &Fwd[i][0], 1,
1840  &IncFieldFwd[i][0], 1);
1841  }
1842 
1843  // Total Field Formulation
1844  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PEC_Forces");
1845  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC_Forces");
1846  CopyBoundaryTrace(IncFieldFwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2,
1847  "PEC_Forces");
1848 
1849  CopyBoundaryTrace(IncFieldFwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1850  "PMC_Forces");
1851  CopyBoundaryTrace(IncFieldFwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1,
1852  "PMC_Forces");
1853  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PMC_Forces");
1854 
1855  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PEC");
1856  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1857  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2, "PEC");
1858 
1859  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PMC");
1860  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1, "PMC");
1861  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PMC");
1862 
1863  Array<OneD, NekDouble> e1Fwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1864  Array<OneD, NekDouble> e1Bwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1865  Array<OneD, NekDouble> e2Fwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1866  Array<OneD, NekDouble> e2Bwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1867  Array<OneD, NekDouble> e3Fwd_cdot_dEe3(nTraceNumPoints, 0.0);
1868  Array<OneD, NekDouble> e3Bwd_cdot_dEe3(nTraceNumPoints, 0.0);
1869 
1870  // Compute numfluxFwd[dir] = (eFwd^[dir] \cdot n \times e^3) * (YimFwd *
1871  // EFwd^3 + YimBwd * EBwd^3 )
1872  ComputeNtimesFz(0, Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], numfluxFwd[0],
1873  numfluxBwd[0]);
1874  ComputeNtimesFz(1, Fwd, Bwd, m_YimFwd[1], m_YimBwd[1], numfluxFwd[1],
1875  numfluxBwd[1]);
1876 
1877  // Compute numfluxFwd[2] = eFwd^3 \cdot ( n1e1 \times ( ZimFwd HFwd + ZimBwd
1878  // HBwd ) ) / 2 {{Z_i}}
1879  ComputeNtimesF12(Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], m_ZimFwd[1],
1880  m_ZimBwd[1], numfluxFwd[2], numfluxBwd[2]);
1881 
1882  // Compute e1Fwd_cdot_ncrossdE = eFwd[dir] \cdot \alpha n \times n \times
1883  // [H] / 2 {{YimFwd}}
1884  ComputeNtimestimesdFz(0, Fwd, Bwd, m_YimFwd[0], m_YimBwd[0],
1885  e1Fwd_cdot_ncrossdH, e1Bwd_cdot_ncrossdH);
1886  ComputeNtimestimesdFz(1, Fwd, Bwd, m_YimFwd[1], m_YimBwd[1],
1887  e2Fwd_cdot_ncrossdH, e2Bwd_cdot_ncrossdH);
1888 
1889  // Compute \alpha [E3] * ( 1/2{{Zim1}} + 1/2{{Zim2}} )
1890  ComputeNtimestimesdF12(Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], m_ZimFwd[1],
1891  m_ZimBwd[1], e3Fwd_cdot_dEe3, e3Bwd_cdot_dEe3);
1892 
1893  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[0], 1, e1Fwd_cdot_ncrossdH,
1894  1, numfluxFwd[0], 1);
1895  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[1], 1, e2Fwd_cdot_ncrossdH,
1896  1, numfluxFwd[1], 1);
1897  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[2], 1, e3Fwd_cdot_dEe3, 1,
1898  numfluxFwd[2], 1);
1899 
1900  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[0], 1, e1Bwd_cdot_ncrossdH,
1901  1, numfluxBwd[0], 1);
1902  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[1], 1, e2Bwd_cdot_ncrossdH,
1903  1, numfluxBwd[1], 1);
1904  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[2], 1, e3Bwd_cdot_dEe3, 1,
1905  numfluxBwd[2], 1);
1906 }

References ComputeNtimesF12(), ComputeNtimesFz(), ComputeNtimestimesdF12(), ComputeNtimestimesdFz(), CopyBoundaryTrace(), Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::eFwdEQNegBwd, GetIncidentField(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, m_YimBwd, m_YimFwd, m_ZimBwd, m_ZimFwd, and Vmath::Svtvp().

Referenced by NumericalMaxwellFlux().

◆ RootMeanSquare()

NekDouble Nektar::SolverUtils::MMFSystem::RootMeanSquare ( const Array< OneD, const NekDouble > &  inarray)
protected

Definition at line 2344 of file MMFSystem.cpp.

2345 {
2346  int Ntot = inarray.size();
2347 
2348  NekDouble reval = 0.0;
2349  for (int i = 0; i < Ntot; ++i)
2350  {
2351  reval = reval + inarray[i] * inarray[i];
2352  }
2353  reval = sqrt(reval / Ntot);
2354 
2355  return reval;
2356 }

References tinysimd::sqrt().

Referenced by Computedemdxicdote(), ComputencdotMF(), ComputeZimYim(), Nektar::MMFAdvection::EvaluateAdvectionVelocity(), Nektar::MMFMaxwell::v_DoSolve(), Nektar::MMFAdvection::v_InitObject(), Nektar::MMFMaxwell::v_InitObject(), and VectorAvgMagnitude().

◆ SetUpMovingFrames()

void Nektar::SolverUtils::MMFSystem::SetUpMovingFrames ( const Array< OneD, const Array< OneD, NekDouble >> &  Anisotropy,
const int  TangentXelem 
)
protected

Definition at line 120 of file MMFSystem.cpp.

123 {
124 
125  int nq = m_fields[0]->GetNpoints();
126  int MFdim = 3;
127 
128  // Construct The Moving Frames
129  m_movingframes = Array<OneD, Array<OneD, NekDouble>>(MFdim);
130 
131  for (int j = 0; j < MFdim; ++j)
132  {
133  m_movingframes[j] = Array<OneD, NekDouble>(m_spacedim * nq, 0.0);
134  }
135 
136  // Read MMF Geom Info
137  std::string conn = "TangentX";
138  m_MMFfactors = Array<OneD, NekDouble>(4);
139 
140  m_session->LoadSolverInfo("MMFDir", conn, "LOCAL");
141 
142  // (x-x_0)^2/a^2 + (y-y_0)^2/b^2 = 1
143  // factors[0] = a
144  // factors[1] = b
145  // factors[2] = x_0
146  // factors[3] = y_0
147  m_session->LoadParameter("MMFCircAxisX", m_MMFfactors[0], 1.0);
148  m_session->LoadParameter("MMFCircAxisY", m_MMFfactors[1], 1.0);
149  m_session->LoadParameter("MMFCircCentreX", m_MMFfactors[2], 0.0);
150  m_session->LoadParameter("MMFCircCentreY", m_MMFfactors[3], 0.0);
151 
152  if (conn == "TangentX")
154  if (conn == "TangentY")
156  if (conn == "TangentXY")
158  if (conn == "TangentZ")
160  if (conn == "TangentCircular")
162  if (conn == "TangentIrregular")
164  if (conn == "TangentNonconvex")
166  if (conn == "LOCAL")
168 
169  /*for (int j=0; j<m_shapedim; j++)
170  {
171  for (int k=0; k<m_spacedim; k++)
172  {
173  std::cout << "before ehsan " << nq << "\t"<< m_shapedim << "\t" <<
174  m_spacedim << "\t" <<"m_movingframes " << m_movingframes[j][k*nq] <<
175  std::endl;
176  }
177  }*/
178  // Get Tangetn vectors from GeomFactors2D, Orthonormalized = true
179  m_fields[0]->GetMovingFrames(m_MMFdir, m_MMFfactors, m_movingframes);
180  /* for (int j=0; j<m_shapedim; j++)
181  {
182  for (int k=0; k<m_spacedim; k++)
183  {
184  std::cout << "ehsan " << nq << "\t"<< m_shapedim << "\t" <<
185  m_spacedim << "\t" <<"m_movingframes " << m_movingframes[j][k*nq] <<
186  std::endl;
187  }
188  }*/
189  // Align the tangentX direction after TangentXelem
190  if (TangentXelem > 0)
191  {
192  Array<OneD, NekDouble> tmp(nq);
193  m_fields[0]->GenerateElementVector(TangentXelem, 1.0, 0.0, tmp);
194  for (int j = 0; j < m_shapedim; j++)
195  {
196  for (int k = 0; k < m_spacedim; k++)
197  {
198  Vmath::Vmul(nq, &tmp[0], 1, &m_movingframes[j][k * nq], 1,
199  &m_movingframes[j][k * nq], 1);
200  }
201  }
202 
203  m_fields[0]->GenerateElementVector(TangentXelem, 0.0, 1.0, tmp);
204  Vmath::Vadd(nq, &tmp[0], 1, &m_movingframes[0][0 * nq], 1,
205  &m_movingframes[0][0 * nq], 1);
206  Vmath::Vadd(nq, &tmp[0], 1, &m_movingframes[1][1 * nq], 1,
207  &m_movingframes[1][1 * nq], 1);
208 
209  int indxtmp = Vmath::Imax(nq, tmp, 1);
210  std::cout << "*** MF in PML Region is aligned as MF1 = ( "
211  << m_movingframes[0][indxtmp] << " , "
212  << m_movingframes[0][nq + indxtmp] << " , "
213  << m_movingframes[0][2 * nq + indxtmp] << " ) "
214  << ", MF2 = ( " << m_movingframes[1][indxtmp] << " , "
215  << m_movingframes[1][nq + indxtmp] << " , "
216  << m_movingframes[1][2 * nq + indxtmp] << " ) " << std::endl;
217  }
218 
219  // Multiply Anisotropy to movingframes
220  for (int j = 0; j < m_shapedim; ++j)
221  {
222  for (int k = 0; k < m_spacedim; ++k)
223  {
224  Vmath::Vmul(nq, &Anisotropy[j][0], 1, &m_movingframes[j][k * nq], 1,
225  &m_movingframes[j][k * nq], 1);
226  }
227  }
228 
229  // Test the moving frames
231 }
void CheckMovingFrames(const Array< OneD, const Array< OneD, NekDouble >> &movingframes)
Definition: MMFSystem.cpp:233
Array< OneD, NekDouble > m_MMFfactors
Definition: MMFSystem.h:160
@ eLOCAL
No Principal direction.
@ eTangentIrregular
Circular around the centre of domain.
@ eTangentX
X coordinate direction.
@ eTangentCircular
Circular around the centre of domain.
@ eTangentNonconvex
Circular around the centre of domain.
@ eTangentXY
XY direction.
@ eTangentZ
Z coordinate direction.
@ eTangentY
Y coordinate direction.
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:918

References CheckMovingFrames(), Nektar::SpatialDomains::eLOCAL, Nektar::SpatialDomains::eTangentCircular, Nektar::SpatialDomains::eTangentIrregular, Nektar::SpatialDomains::eTangentNonconvex, Nektar::SpatialDomains::eTangentX, Nektar::SpatialDomains::eTangentXY, Nektar::SpatialDomains::eTangentY, Nektar::SpatialDomains::eTangentZ, Vmath::Imax(), Nektar::SolverUtils::EquationSystem::m_fields, m_MMFdir, m_MMFfactors, m_movingframes, Nektar::SolverUtils::EquationSystem::m_session, m_shapedim, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vadd(), and Vmath::Vmul().

Referenced by MMFInitObject().

◆ UpwindMaxwellFlux1D()

void Nektar::SolverUtils::MMFSystem::UpwindMaxwellFlux1D ( Array< OneD, Array< OneD, NekDouble >> &  physfield,
Array< OneD, Array< OneD, NekDouble >> &  numfluxFwd,
Array< OneD, Array< OneD, NekDouble >> &  numfluxBwd 
)
protected

Definition at line 1461 of file MMFSystem.cpp.

1465 {
1466  int i;
1467  int nTraceNumPoints = GetTraceTotPoints();
1468  int nvar = 2;
1469 
1470  // get temporary arrays
1471  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1472  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1473 
1474  for (i = 0; i < nvar; ++i)
1475  {
1476  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1477  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1478  }
1479 
1480  // get the physical values at the trace from the dependent variables
1481  for (i = 0; i < nvar; ++i)
1482  {
1483  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1484  }
1485 
1486  // E = 0 at the boundaries
1487  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1488 
1489  // d H / d n = 0 at the boundaries
1490  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1491 
1492  Array<OneD, NekDouble> dE(nTraceNumPoints);
1493  Array<OneD, NekDouble> dH(nTraceNumPoints);
1494 
1495  Vmath::Vsub(nTraceNumPoints, &Fwd[0][0], 1, &Bwd[0][0], 1, &dE[0], 1);
1496  Vmath::Vsub(nTraceNumPoints, &Fwd[1][0], 1, &Bwd[1][0], 1, &dH[0], 1);
1497 
1498  NekDouble nx, AverZ, AverY, AverZH, AverYE;
1499  for (i = 0; i < nTraceNumPoints; ++i)
1500  {
1501  nx = m_traceNormals[0][i];
1502  AverZ = m_ZimFwd[0][i] + m_ZimBwd[0][i];
1503  AverY = m_YimFwd[0][i] + m_YimBwd[0][i];
1504  AverZH = m_ZimFwd[0][i] * Fwd[1][i] + m_ZimBwd[0][i] * Bwd[1][i];
1505  AverYE = m_YimFwd[0][i] * Fwd[0][i] + m_YimBwd[0][i] * Bwd[0][i];
1506 
1507  numfluxFwd[0][i] = nx / AverZ * (AverZH - nx * dE[i]);
1508  numfluxFwd[1][i] = nx / AverY * (AverYE - nx * dH[i]);
1509 
1510  numfluxBwd[0][i] = nx / AverZ * (AverZH - nx * dE[i]);
1511  numfluxBwd[1][i] = nx / AverY * (AverYE - nx * dH[i]);
1512  }
1513 }

References CopyBoundaryTrace(), Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::eFwdEQNegBwd, Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_traceNormals, m_YimBwd, m_YimFwd, m_ZimBwd, m_ZimFwd, and Vmath::Vsub().

Referenced by NumericalMaxwellFlux().

◆ v_GenerateSummary()

void Nektar::SolverUtils::MMFSystem::v_GenerateSummary ( SummaryList s)
overridevirtual

Print a summary of time stepping parameters.

Prints a summary with some information regards the time-stepping.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Reimplemented in Nektar::MMFSWE, Nektar::MMFMaxwell, Nektar::MMFDiffusion, and Nektar::MMFAdvection.

Definition at line 2469 of file MMFSystem.cpp.

2470 {
2471  int nq = m_fields[0]->GetNpoints();
2473 
2474  AddSummaryItem(s, "Total grids", nq);
2475  AddSummaryItem(s, "Shape Dimension", m_shapedim);
2477  if (m_surfaceType == eSphere)
2478  {
2479  NekDouble MeshError;
2480 
2481  Array<OneD, NekDouble> x(nq);
2482  Array<OneD, NekDouble> y(nq);
2483  Array<OneD, NekDouble> z(nq);
2484  Array<OneD, NekDouble> rad(nq, 0.0);
2485 
2486  m_fields[0]->GetCoords(x, y, z);
2487 
2488  Vmath::Vvtvp(nq, x, 1, x, 1, rad, 1, rad, 1);
2489  Vmath::Vvtvp(nq, y, 1, y, 1, rad, 1, rad, 1);
2490  Vmath::Vvtvp(nq, z, 1, z, 1, rad, 1, rad, 1);
2491  Vmath::Vsqrt(nq, rad, 1, rad, 1);
2492 
2493  Vmath::Sadd(nq, -1.0, rad, 1, rad, 1);
2494  Vmath::Vabs(nq, rad, 1, rad, 1);
2495 
2496  MeshError = m_fields[0]->Integral(rad);
2497  SolverUtils::AddSummaryItem(s, "Mesh Error", MeshError);
2498  }
2499 
2501  AddSummaryItem(s, "SmoothFactor", m_SmoothFactor);
2502  AddSummaryItem(s, "SFinit", m_SFinit);
2503 
2504  if (fabs(m_alpha - 1.0) > 0.001)
2505  {
2506  AddSummaryItem(s, "Alpha", m_alpha);
2507  }
2508 
2509  if (m_Incfreq > 0.0)
2510  {
2511  AddSummaryItem(s, "Incfreq", m_Incfreq);
2512  }
2513 
2514  if (m_MMFdir == 4)
2515  {
2516  AddSummaryItem(s, "MMFCircAxisX", m_MMFfactors[0]);
2517  AddSummaryItem(s, "MMFCircAxisY", m_MMFfactors[1]);
2518  AddSummaryItem(s, "MMFCircCentreX", m_MMFfactors[2]);
2519  AddSummaryItem(s, "MMFCircCentreY", m_MMFfactors[3]);
2520  }
2521 }
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
static NekDouble rad(NekDouble x, NekDouble y)
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
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add scalar y = alpha + x.
Definition: Vmath.cpp:384

References Nektar::SolverUtils::AddSummaryItem(), Nektar::SolverUtils::eSphere, Nektar::SpatialDomains::GeomMMFMap, m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_Incfreq, m_MMFdir, m_MMFfactors, m_SFinit, m_shapedim, m_SmoothFactor, m_surfaceType, Nektar::LibUtilities::rad(), Vmath::Sadd(), Nektar::SolverUtils::SurfaceTypeMap, Nektar::SolverUtils::UnsteadySystem::v_GenerateSummary(), Vmath::Vabs(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by Nektar::MMFAdvection::v_GenerateSummary(), Nektar::MMFDiffusion::v_GenerateSummary(), Nektar::MMFMaxwell::v_GenerateSummary(), and Nektar::MMFSWE::v_GenerateSummary().

◆ VectorAvgMagnitude()

NekDouble Nektar::SolverUtils::MMFSystem::VectorAvgMagnitude ( const Array< OneD, const Array< OneD, NekDouble >> &  inarray)
protected

Definition at line 2358 of file MMFSystem.cpp.

2360 {
2361  int nq = inarray[0].size();
2362 
2363  Array<OneD, NekDouble> tmp(nq, 0.0);
2364  for (int k = 0; k < m_spacedim; k++)
2365  {
2366  Vmath::Vvtvp(nq, &inarray[k][0], 1, &inarray[k][0], 1, &tmp[0], 1,
2367  &tmp[0], 1);
2368  }
2369  Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
2370 
2371  return RootMeanSquare(tmp);
2372 }

References Nektar::SolverUtils::EquationSystem::m_spacedim, RootMeanSquare(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by ComputeMFtrace(), and ComputeNtimesMF().

◆ VectorCrossProd() [1/2]

void Nektar::SolverUtils::MMFSystem::VectorCrossProd ( const Array< OneD, const Array< OneD, NekDouble >> &  v1,
const Array< OneD, const Array< OneD, NekDouble >> &  v2,
Array< OneD, Array< OneD, NekDouble >> &  v3 
)
protected

Computes the vector cross-product in 3D of v1 and v2, storing the result in v3.

Parameters
v1First input vector.
v2Second input vector.
v3Output vector computed to be orthogonal to both v1 and v2.

Definition at line 699 of file MMFSystem.cpp.

703 {
704  ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
705  ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
706  ASSERTL0(v3.size() == 3, "Output vector has dimension not equal to 3.");
707 
708  int nq = v1[0].size();
709  Array<OneD, NekDouble> temp(nq);
710 
711  Vmath::Vmul(nq, v1[2], 1, v2[1], 1, temp, 1);
712  Vmath::Vvtvm(nq, v1[1], 1, v2[2], 1, temp, 1, v3[0], 1);
713 
714  Vmath::Vmul(nq, v1[0], 1, v2[2], 1, temp, 1);
715  Vmath::Vvtvm(nq, v1[2], 1, v2[0], 1, temp, 1, v3[1], 1);
716 
717  Vmath::Vmul(nq, v1[1], 1, v2[0], 1, temp, 1);
718  Vmath::Vvtvm(nq, v1[0], 1, v2[1], 1, temp, 1, v3[2], 1);
719 }
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 minus vector): z = w*x - y
Definition: Vmath.cpp:598

References ASSERTL0, Vmath::Vmul(), and Vmath::Vvtvm().

Referenced by ComputencdotMF(), ComputeNtimesF12(), ComputeNtimesMF(), ComputeNtimestimesdFz(), DeriveCrossProductMF(), and Nektar::MMFAdvection::EvaluateAdvectionVelocity().

◆ VectorCrossProd() [2/2]

void Nektar::SolverUtils::MMFSystem::VectorCrossProd ( const Array< OneD, NekDouble > &  v1,
const Array< OneD, NekDouble > &  v2,
Array< OneD, NekDouble > &  v3 
)
protected

Definition at line 721 of file MMFSystem.cpp.

724 {
725  ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
726  ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
727  ASSERTL0(v3.size() == 3, "Output vector has dimension not equal to 3.");
728 
729  v3[0] = v1[1] * v2[2] - v1[2] * v2[1];
730  v3[1] = v1[2] * v2[0] - v1[0] * v2[2];
731  v3[2] = v1[0] * v2[1] - v1[1] * v2[0];
732 }

References ASSERTL0.

◆ VectorDotProd()

void Nektar::SolverUtils::MMFSystem::VectorDotProd ( const Array< OneD, const Array< OneD, NekDouble >> &  v1,
const Array< OneD, const Array< OneD, NekDouble >> &  v2,
Array< OneD, NekDouble > &  v3 
)
protected

Definition at line 676 of file MMFSystem.cpp.

680 {
681  int coordim = v1.size();
682  int nq = v1[0].size();
683 
684  v3 = Array<OneD, NekDouble>(nq, 0.0);
685  for (int i = 0; i < coordim; ++i)
686  {
687  Vmath::Vvtvp(nq, &v1[i][0], 1, &v2[i][0], 1, &v3[0], 1, &v3[0], 1);
688  }
689 }

References Vmath::Vvtvp().

Referenced by ComputencdotMF().

Member Data Documentation

◆ m_alpha

NekDouble Nektar::SolverUtils::MMFSystem::m_alpha
protected

◆ m_CurlMF

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::MMFSystem::m_CurlMF
protected

◆ m_dedxi_cdot_e

Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble> > > > Nektar::SolverUtils::MMFSystem::m_dedxi_cdot_e
protected

Definition at line 220 of file MMFSystem.h.

Referenced by AdddedtMaxwell(), and Computedemdxicdote().

◆ m_DivMF

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_DivMF
protected

◆ m_epsvec

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_epsvec
protected

◆ m_Incfreq

NekDouble Nektar::SolverUtils::MMFSystem::m_Incfreq
protected

◆ m_IncType

IncType Nektar::SolverUtils::MMFSystem::m_IncType

◆ m_MFlength

Array<OneD, NekDouble> Nektar::SolverUtils::MMFSystem::m_MFlength
protected

Definition at line 224 of file MMFSystem.h.

◆ m_MFtraceBwd

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::MMFSystem::m_MFtraceBwd
protected

◆ m_MFtraceFwd

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::MMFSystem::m_MFtraceFwd
protected

◆ m_MMFdir

SpatialDomains::GeomMMF Nektar::SolverUtils::MMFSystem::m_MMFdir
protected

◆ m_MMFfactors

Array<OneD, NekDouble> Nektar::SolverUtils::MMFSystem::m_MMFfactors

◆ m_movingframes

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_movingframes
protected

◆ m_muvec

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_muvec
protected

◆ m_ncdotMFBwd

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd
protected

◆ m_ncdotMFFwd

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd
protected

◆ m_negepsvecminus1

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_negepsvecminus1
protected

Definition at line 215 of file MMFSystem.h.

Referenced by Nektar::MMFMaxwell::DoOdeRhs(), and Nektar::MMFMaxwell::v_InitObject().

◆ m_negmuvecminus1

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_negmuvecminus1
protected

Definition at line 216 of file MMFSystem.h.

Referenced by Nektar::MMFMaxwell::DoOdeRhs(), and Nektar::MMFMaxwell::v_InitObject().

◆ m_nperpcdotMFBwd

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_nperpcdotMFBwd
protected

Definition at line 193 of file MMFSystem.h.

Referenced by ComputencdotMF(), and Nektar::MMFSWE::NumericalSWEFlux().

◆ m_nperpcdotMFFwd

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_nperpcdotMFFwd
protected

◆ m_ntimes_ntimesMFBwd

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::MMFSystem::m_ntimes_ntimesMFBwd
protected

Definition at line 205 of file MMFSystem.h.

Referenced by ComputeNtimesMF().

◆ m_ntimes_ntimesMFFwd

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::MMFSystem::m_ntimes_ntimesMFFwd
protected

Definition at line 204 of file MMFSystem.h.

Referenced by ComputeNtimesMF().

◆ m_ntimesMFBwd

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::MMFSystem::m_ntimesMFBwd
protected

Definition at line 203 of file MMFSystem.h.

Referenced by ComputeNtimesFz(), and ComputeNtimesMF().

◆ m_ntimesMFFwd

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::MMFSystem::m_ntimesMFFwd
protected

◆ m_pi

NekDouble Nektar::SolverUtils::MMFSystem::m_pi

◆ m_PolType

PolType Nektar::SolverUtils::MMFSystem::m_PolType

◆ m_SFinit

NekDouble Nektar::SolverUtils::MMFSystem::m_SFinit
protected

Definition at line 184 of file MMFSystem.h.

Referenced by GetIncidentField(), MMFInitObject(), and v_GenerateSummary().

◆ m_shapedim

int Nektar::SolverUtils::MMFSystem::m_shapedim

◆ m_SmoothFactor

int Nektar::SolverUtils::MMFSystem::m_SmoothFactor
protected

Definition at line 183 of file MMFSystem.h.

Referenced by GetIncidentField(), MMFInitObject(), and v_GenerateSummary().

◆ m_surfaceNormal

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_surfaceNormal
protected

Definition at line 187 of file MMFSystem.h.

◆ m_surfaceType

SurfaceType Nektar::SolverUtils::MMFSystem::m_surfaceType

◆ m_TestMaxwellType

TestMaxwellType Nektar::SolverUtils::MMFSystem::m_TestMaxwellType

◆ m_upwindType

UpwindType Nektar::SolverUtils::MMFSystem::m_upwindType

◆ m_YimBwd

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_YimBwd
protected

◆ m_YimFwd

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_YimFwd
protected

◆ m_ZimBwd

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_ZimBwd
protected

◆ m_ZimFwd

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::MMFSystem::m_ZimFwd
protected