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

#include <MMFMaxwell.h>

Inheritance diagram for Nektar::MMFMaxwell:
[legend]

Public Member Functions

virtual void v_InitObject (bool DeclareFields=true)
 Initialise the object. More...
 
virtual void v_DoSolve ()
 Solves an unsteady problem. More...
 
virtual ~MMFMaxwell ()
 Destructor. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::MMFSystem
SOLVER_UTILS_EXPORT MMFSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 
virtual SOLVER_UTILS_EXPORT ~MMFSystem ()
 
SOLVER_UTILS_EXPORT void MMFInitObject (const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem=-1)
 
SOLVER_UTILS_EXPORT void CopyBoundaryTrace (const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const BoundaryCopyType BDCopyType, const int var=0, const std::string btype="NoUserDefined")
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
virtual SOLVER_UTILS_EXPORT ~UnsteadySystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble >> &inarray)
 Calculate the larger time-step mantaining the problem stable. More...
 
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, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT Array< OneD, const Array< OneD, NekDouble > > GetTraceNormals ()
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp ()
 Virtual function to identify if operator is negated in DoSolve. More...
 

Static Public Member Functions

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

Public Attributes

CloakType m_CloakType
 
SourceType m_SourceType
 
bool m_DispersiveCloak
 
- Public Attributes inherited from Nektar::SolverUtils::MMFSystem
NekDouble m_pi
 
int m_shapedim
 
SurfaceType m_surfaceType
 
UpwindType m_upwindType
 
TestMaxwellType m_TestMaxwellType
 
PolType m_PolType
 
IncType m_IncType
 
Array< OneD, 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...
 

Static Public Attributes

static std::string className
 Name of class. More...
 
- Static Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
static std::string cmdSetStartTime
 
static std::string cmdSetStartChkNum
 

Protected Member Functions

 MMFMaxwell (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Session reader. More...
 
void DoOdeRhs (const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
 Compute the RHS. More...
 
void DoOdeProjection (const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
 Compute the projection. More...
 
void AddGreenDerivCompensate (const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
 
void WeakDGMaxwellDirDeriv (const Array< OneD, const Array< OneD, NekDouble >> &InField, Array< OneD, Array< OneD, NekDouble >> &OutField, const NekDouble time=0.0)
 Calculate weak DG advection in the form \( \langle\phi, \hat{F}\cdot n\rangle - (\nabla \phi \cdot F) \). More...
 
NekDouble ComputeEnergyDensity (Array< OneD, Array< OneD, NekDouble >> &fields)
 
Array< OneD, NekDoubleTestMaxwell1D (const NekDouble time, unsigned int field)
 
Array< OneD, NekDoubleTestMaxwell2DPEC (const NekDouble time, unsigned int field, const SolverUtils::PolType Polarization)
 
Array< OneD, NekDoubleTestMaxwell2DPMC (const NekDouble time, unsigned int field, const SolverUtils::PolType Polarization)
 
Array< OneD, NekDoubleTestMaxwellSphere (const NekDouble time, const NekDouble omega, unsigned int field)
 
void Printout_SurfaceCurrent (Array< OneD, Array< OneD, NekDouble >> &fields, const int time)
 
Array< OneD, NekDoubleComputeSurfaceCurrent (const int time, const Array< OneD, const Array< OneD, NekDouble >> &fields)
 
void GenerateSigmaPML (const NekDouble PMLthickness, const NekDouble PMLstart, const NekDouble PMLmaxsigma, Array< OneD, Array< OneD, NekDouble >> &SigmaPML)
 
void ComputeMaterialVector (Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
 
void ComputeMaterialOpticalCloak (const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec, const bool Dispersion=false)
 
void ComputeMaterialMicroWaveCloak (const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
 
void Checkpoint_TotalFieldOutput (const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
 
void Checkpoint_PlotOutput (const int n, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
 
void Checkpoint_TotPlotOutput (const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
 
void Checkpoint_EDFluxOutput (const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
 
void Checkpoint_EnergyOutput (const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
 
Array< OneD, NekDoubleGaussianPulse (const NekDouble time, const NekDouble Psx, const NekDouble Psy, const NekDouble Psz, const NekDouble Gaussianradius)
 
void AddPML (const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
 
Array< OneD, NekDoubleEvaluateCoriolis ()
 
void AddCoriolis (Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
 
Array< OneD, NekDoubleComputeRadCloak (const int Nlayer=5)
 
virtual void v_GenerateSummary (SolverUtils::SummaryList &s)
 Print Summary. More...
 
virtual void v_SetInitialConditions (const NekDouble initialtime, bool dumpInitialConditions, const int domain)
 
virtual void v_EvaluateExactSolution (unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 
void print_MMF (Array< OneD, Array< OneD, NekDouble >> &inarray)
 
- Protected Member Functions inherited from Nektar::SolverUtils::MMFSystem
void SetUpMovingFrames (const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem)
 
void CheckMovingFrames (const Array< OneD, const Array< OneD, NekDouble >> &movingframes)
 
SOLVER_UTILS_EXPORT void ComputencdotMF ()
 
SOLVER_UTILS_EXPORT void ComputeDivCurlMF ()
 
SOLVER_UTILS_EXPORT void ComputeMFtrace ()
 
SOLVER_UTILS_EXPORT void VectorDotProd (const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, NekDouble > &v3)
 
SOLVER_UTILS_EXPORT void VectorCrossProd (const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, Array< OneD, NekDouble >> &v3)
 
SOLVER_UTILS_EXPORT void VectorCrossProd (const Array< OneD, NekDouble > &v1, const Array< OneD, NekDouble > &v2, Array< OneD, NekDouble > &v3)
 
SOLVER_UTILS_EXPORT void ComputeCurl (const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
 
SOLVER_UTILS_EXPORT Array< OneD, 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...
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise ()
 Sets up initial conditions. More...
 
virtual SOLVER_UTILS_EXPORT void v_AppendOutput1D (Array< OneD, Array< OneD, NekDouble >> &solution1D)
 Print the solution at each solution point in a txt file. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_GetTimeStep (const Array< OneD, const Array< OneD, NekDouble >> &inarray)
 Return the timestep to be used for the next step in the time-marching loop. More...
 
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_RequireFwdTrans ()
 
virtual SOLVER_UTILS_EXPORT void v_SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time, int &nchk)
 
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff (const Array< OneD, Array< OneD, NekDouble >> vel, StdRegions::VarCoeffMap &varCoeffMap)
 Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h time velocity. More...
 
virtual SOLVER_UTILS_EXPORT bool UpdateTimeStepCheck ()
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members. More...
 
virtual SOLVER_UTILS_EXPORT 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_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

int m_ElemtGroup0
 
int m_ElemtGroup1
 
int m_boundaryforSF
 
int m_PrintoutSurfaceCurrent
 
int m_AddPML
 
int m_PMLorder
 
int m_AddRotation
 
bool m_Cloaking
 
NekDouble m_CloakNlayer
 
NekDouble m_Cloakraddelta
 
NekDouble m_wp2Tol
 
Array< OneD, NekDoublem_wp2
 
Array< OneD, NekDoublem_SourceVector
 
NekDouble m_Psx
 
NekDouble m_Psy
 
NekDouble m_Psz
 
NekDouble m_PSduration
 
NekDouble m_Gaussianradius
 
Array< OneD, Array< OneD, NekDouble > > m_CrossProductMF
 
NekDouble m_freq
 
NekDouble m_n1
 
NekDouble m_n2
 
NekDouble m_n3
 
Array< OneD, NekDoublem_varepsilon
 
Array< OneD, NekDoublem_mu
 
int m_TestPML
 
int m_PMLelement
 
int m_RecPML
 
NekDouble m_PMLthickness
 
NekDouble m_PMLstart
 
NekDouble m_PMLmaxsigma
 
Array< OneD, Array< OneD, NekDouble > > m_SigmaPML
 
int m_NoInc
 
Array< OneD, NekDoublem_coriolis
 
- Protected Attributes inherited from Nektar::SolverUtils::MMFSystem
NekDouble m_alpha
 
NekDouble m_Incfreq
 
int m_SmoothFactor
 
NekDouble m_SFinit
 
Array< OneD, Array< OneD, NekDouble > > m_movingframes
 
Array< OneD, Array< OneD, NekDouble > > m_surfaceNormal
 
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
 
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
 
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
 
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFBwd
 
Array< OneD, Array< OneD, NekDouble > > m_DivMF
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFFwd
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFBwd
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFFwd
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFBwd
 
Array< OneD, Array< OneD, NekDouble > > m_ZimFwd
 
Array< OneD, Array< OneD, NekDouble > > m_ZimBwd
 
Array< OneD, Array< OneD, NekDouble > > m_YimFwd
 
Array< OneD, Array< OneD, NekDouble > > m_YimBwd
 
Array< OneD, Array< OneD, NekDouble > > m_epsvec
 
Array< OneD, Array< OneD, NekDouble > > m_muvec
 
Array< OneD, Array< OneD, NekDouble > > m_negepsvecminus1
 
Array< OneD, Array< OneD, NekDouble > > m_negmuvecminus1
 
Array< OneD, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > > m_dedxi_cdot_e
 
SpatialDomains::GeomMMF m_MMFdir
 
Array< OneD, NekDoublem_MFlength
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
int m_infosteps
 Number of time steps between outputting status information. More...
 
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
int m_nanSteps
 
LibUtilities::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
 
bool m_root
 
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_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
Array< OneD, 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...
 

Friends

class MemoryManager< MMFMaxwell >
 

Additional Inherited Members

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

Detailed Description

Definition at line 73 of file MMFMaxwell.h.

Constructor & Destructor Documentation

◆ ~MMFMaxwell()

Nektar::MMFMaxwell::~MMFMaxwell ( )
virtual

Destructor.

Unsteady linear advection equation destructor.

Definition at line 445 of file MMFMaxwell.cpp.

446 {
447 }

◆ MMFMaxwell()

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

Session reader.

Definition at line 56 of file MMFMaxwell.cpp.

58  : UnsteadySystem(pSession, pGraph), MMFSystem(pSession, pGraph)
59 {
60 }
SOLVER_UTILS_EXPORT MMFSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition: MMFSystem.cpp:43
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.

Member Function Documentation

◆ AddCoriolis()

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

Definition at line 2866 of file MMFMaxwell.cpp.

2868 {
2869  int nq = physarray[0].size();
2870 
2871  Array<OneD, NekDouble> tmp(nq);
2872 
2873  int indx;
2874  for (int j = 0; j < m_shapedim; ++j)
2875  {
2876  if (j == 0)
2877  {
2878  indx = 2;
2879  }
2880 
2881  else if (j == 1)
2882  {
2883  indx = 1;
2884  }
2885 
2886  Vmath::Vmul(nq, m_coriolis, 1, physarray[indx], 1, tmp, 1);
2887 
2888  switch (m_PolType)
2889  {
2891  {
2892  Vmath::Vmul(nq, m_muvec[indx], 1, tmp, 1, tmp, 1);
2893  }
2894  break;
2895 
2897  {
2898  Vmath::Vmul(nq, m_epsvec[indx], 1, tmp, 1, tmp, 1);
2899  }
2900  break;
2901 
2902  default:
2903  break;
2904  }
2905 
2906  if (j == 1)
2907  {
2908  Vmath::Neg(nq, tmp, 1);
2909  }
2910 
2911  Vmath::Vadd(nq, tmp, 1, outarray[j], 1, outarray[j], 1);
2912  }
2913 }
Array< OneD, NekDouble > m_coriolis
Definition: MMFMaxwell.h:144
Array< OneD, Array< OneD, NekDouble > > m_muvec
Definition: MMFSystem.h:213
Array< OneD, Array< OneD, NekDouble > > m_epsvec
Definition: MMFSystem.h:212
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 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, m_coriolis, Nektar::SolverUtils::MMFSystem::m_epsvec, Nektar::SolverUtils::MMFSystem::m_muvec, Nektar::SolverUtils::MMFSystem::m_PolType, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Neg(), Vmath::Vadd(), and Vmath::Vmul().

Referenced by DoOdeRhs().

◆ AddGreenDerivCompensate()

void Nektar::MMFMaxwell::AddGreenDerivCompensate ( const Array< OneD, const Array< OneD, NekDouble >> &  physarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray 
)
protected

Definition at line 1169 of file MMFMaxwell.cpp.

1172 {
1173  // routine works for both primitive and conservative formulations
1174  int ncoeffs = outarray[0].size();
1175  int nq = physarray[0].size();
1176 
1177  Array<OneD, NekDouble> tmp(nq);
1178  Array<OneD, NekDouble> tmpc(ncoeffs);
1179 
1180  Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
1181  for (int j = 0; j < m_shapedim; ++j)
1182  {
1183  fluxvector[j] = Array<OneD, NekDouble>(nq);
1184  }
1185 
1186  // m_CurlMF[0][0] = e^3 \cdot (\nabla \times e^1) [ NEW m_CurlMF[0][2] ]
1187  // m_CurlMF[0][1] = 0.0
1188  // m_CurlMF[1][0] = 0.0,
1189  // m_CurlMF[1][1] = e^3 \cdot (\nabla \times e^2) [ NEW m_CurlMF[1][2] ]
1190  // m_CurlMF[2][0] = e^1 \cdot (\nabla \times e^3) [ NEW m_CurlMF[2][0] ]
1191  // m_CurlMF[2][1] = e^2 \cdot (\nabla \times e^3) [ NEW m_CurlMF[2][1] ]
1192 
1193  int var;
1194 
1195  switch (m_TestMaxwellType)
1196  {
1204  {
1205  var = 0;
1206  GetMaxwellFluxVector(var, physarray, fluxvector);
1207  Vmath::Vmul(nq, &fluxvector[0][0], 1, &m_CurlMF[0][2][0], 1,
1208  &tmp[0], 1);
1209  m_fields[var]->IProductWRTBase(tmp, tmpc);
1210  Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1211 
1212  var = 1;
1213  GetMaxwellFluxVector(var, physarray, fluxvector);
1214  Vmath::Vmul(nq, &fluxvector[1][0], 1, &m_CurlMF[1][2][0], 1,
1215  &tmp[0], 1);
1216  Vmath::Neg(nq, tmp, 1);
1217  m_fields[var]->IProductWRTBase(tmp, tmpc);
1218  Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1219 
1220  var = 2;
1221  GetMaxwellFluxVector(var, physarray, fluxvector);
1222  Vmath::Vmul(nq, &fluxvector[0][0], 1, &m_CurlMF[2][0][0], 1,
1223  &tmp[0], 1);
1224  Vmath::Vvtvm(nq, &fluxvector[1][0], 1, &m_CurlMF[2][1][0], 1,
1225  &tmp[0], 1, &tmp[0], 1);
1226  m_fields[var]->IProductWRTBase(tmp, tmpc);
1227  Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1228  }
1229  break;
1230 
1231  default:
1232  break;
1233  }
1234 }
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT void GetMaxwellFluxVector(const int var, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
Definition: MMFSystem.cpp:1609
TestMaxwellType m_TestMaxwellType
Definition: MMFSystem.h:156
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Definition: MMFSystem.h:196
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector plus vector): z = w*x - y
Definition: Vmath.cpp:598

References Nektar::SolverUtils::eELF2DSurface, Nektar::SolverUtils::eMaxwellSphere, Nektar::SolverUtils::eScatField2D, Nektar::SolverUtils::eTestMaxwell2DPEC, Nektar::SolverUtils::eTestMaxwell2DPECAVGFLUX, Nektar::SolverUtils::eTestMaxwell2DPMC, Nektar::SolverUtils::eTotField2D, Nektar::SolverUtils::MMFSystem::GetMaxwellFluxVector(), Nektar::SolverUtils::MMFSystem::m_CurlMF, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_shapedim, Nektar::SolverUtils::MMFSystem::m_TestMaxwellType, Vmath::Neg(), Vmath::Vadd(), Vmath::Vmul(), and Vmath::Vvtvm().

Referenced by DoOdeRhs().

◆ AddPML()

void Nektar::MMFMaxwell::AddPML ( const Array< OneD, const Array< OneD, NekDouble >> &  physarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray 
)
protected

Definition at line 2422 of file MMFMaxwell.cpp.

2425 {
2426  int nq = m_fields[0]->GetTotPoints();
2427  Array<OneD, NekDouble> tmp(nq);
2428 
2429  Array<OneD, NekDouble> Sigma0plus1Neg(nq);
2430  Array<OneD, NekDouble> Sigma0minus1(nq);
2431  Array<OneD, NekDouble> Sigma1minus0(nq);
2432 
2433  Vmath::Vsub(nq, &m_SigmaPML[1][0], 1, &m_SigmaPML[0][0], 1,
2434  &Sigma1minus0[0], 1);
2435  Vmath::Vsub(nq, &m_SigmaPML[0][0], 1, &m_SigmaPML[1][0], 1,
2436  &Sigma0minus1[0], 1);
2437  Vmath::Vadd(nq, &m_SigmaPML[0][0], 1, &m_SigmaPML[1][0], 1,
2438  &Sigma0plus1Neg[0], 1);
2439  Vmath::Neg(nq, Sigma0plus1Neg, 1);
2440 
2441  switch (m_PolType)
2442  {
2444  {
2445  int indxH0 = 0;
2446  int indxH1 = 1;
2447  int indxE2 = 2;
2448  int indxQ0 = 3;
2449  int indxQ1 = 4;
2450  int indxP2 = 5;
2451 
2452  // dH0/dt: Add (sigma_0 - \sigma_1) H0 + Q0
2453  Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxH0][0], 1,
2454  &outarray[indxH0][0], 1, &outarray[indxH0][0], 1);
2455  Vmath::Vadd(nq, &physarray[indxQ0][0], 1, &outarray[indxH0][0], 1,
2456  &outarray[indxH0][0], 1);
2457 
2458  // dH1/dt: Add (sigma_1 - \sigma_0) H1 + Q1
2459  Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxH1][0], 1,
2460  &outarray[indxH1][0], 1, &outarray[indxH1][0], 1);
2461  Vmath::Vadd(nq, &physarray[indxQ1][0], 1, &outarray[indxH1][0], 1,
2462  &outarray[indxH1][0], 1);
2463 
2464  // dHz/dt: Add -(\sigma_0 + \sigma_1) Ez + Pz
2465  Vmath::Vvtvp(nq, &Sigma0plus1Neg[0], 1, &physarray[indxE2][0], 1,
2466  &outarray[indxE2][0], 1, &outarray[indxE2][0], 1);
2467  Vmath::Vadd(nq, &physarray[indxP2][0], 1, &outarray[indxE2][0], 1,
2468  &outarray[indxE2][0], 1);
2469 
2470  // dQ0/dt: Assign -\sigma_0 * ( Q0 + (\sigma_0 - \sigma_1) * H0 )
2471  Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxH0][0], 1,
2472  &physarray[indxQ0][0], 1, &outarray[indxQ0][0], 1);
2473  Vmath::Vmul(nq, &m_SigmaPML[0][0], 1, &outarray[indxQ0][0], 1,
2474  &outarray[indxQ0][0], 1);
2475  Vmath::Neg(nq, &outarray[indxQ0][0], 1);
2476 
2477  // dQ1/dt: Assign -\sigma_1 * ( Q1 + (\sigma_1 - \sigma_0) * H1 )
2478  Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxH1][0], 1,
2479  &physarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2480  Vmath::Vmul(nq, &m_SigmaPML[1][0], 1, &outarray[indxQ1][0], 1,
2481  &outarray[indxQ1][0], 1);
2482  Vmath::Neg(nq, &outarray[indxQ1][0], 1);
2483 
2484  if (m_DispersiveCloak)
2485  {
2486  Vmath::Vvtvp(nq, &m_wp2[0], 1, &physarray[indxH1][0], 1,
2487  &outarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2488  }
2489 
2490  // dP3/dt: Assign - \sigma_1 * \sigma_2 * E_z
2491  Vmath::Vmul(nq, &m_SigmaPML[0][0], 1, &m_SigmaPML[1][0], 1,
2492  &outarray[indxP2][0], 1);
2493  Vmath::Vmul(nq, &physarray[indxE2][0], 1, &outarray[indxP2][0], 1,
2494  &outarray[indxP2][0], 1);
2495  Vmath::Neg(nq, &outarray[indxP2][0], 1);
2496  }
2497  break;
2498 
2500  {
2501  int indxE0 = 0;
2502  int indxE1 = 1;
2503  int indxH2 = 2;
2504  int indxQ0 = 3;
2505  int indxQ1 = 4;
2506  int indxP2 = 5;
2507 
2508  // dE0/dt: Add (sigma_0 - \sigma_1) E0 - Q0
2509  Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxE0][0], 1,
2510  &outarray[indxE0][0], 1, &outarray[indxE0][0], 1);
2511  Vmath::Vsub(nq, &outarray[indxE0][0], 1, &physarray[indxQ0][0], 1,
2512  &outarray[indxE0][0], 1);
2513 
2514  // dE1/dt: Add (sigma_1 - \sigma_0) E1 - Q1
2515  Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxE1][0], 1,
2516  &outarray[indxE1][0], 1, &outarray[indxE1][0], 1);
2517  Vmath::Vsub(nq, &outarray[indxE1][0], 1, &physarray[indxQ1][0], 1,
2518  &outarray[indxE1][0], 1);
2519 
2520  // dHz/dt: Add -(\sigma_0 + \sigma_1) Hz - Pz
2521  Vmath::Vvtvp(nq, &Sigma0plus1Neg[0], 1, &physarray[indxH2][0], 1,
2522  &outarray[indxH2][0], 1, &outarray[indxH2][0], 1);
2523  Vmath::Vsub(nq, &outarray[indxH2][0], 1, &physarray[indxP2][0], 1,
2524  &outarray[indxH2][0], 1);
2525 
2526  // dQ0/dt: Assign -\sigma_0 * ( Q0 + (\sigma_1 - \sigma_0) * E0 )
2527  Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxE0][0], 1,
2528  &physarray[indxQ0][0], 1, &outarray[indxQ0][0], 1);
2529  Vmath::Vmul(nq, &m_SigmaPML[0][0], 1, &outarray[indxQ0][0], 1,
2530  &outarray[indxQ0][0], 1);
2531  Vmath::Neg(nq, &outarray[indxQ0][0], 1);
2532 
2533  // dQ1/dt: Assign -\sigma_1 * ( Q1 + (\sigma_0 - \sigma_1) * E1 )
2534  Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxE1][0], 1,
2535  &physarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2536  Vmath::Vmul(nq, &m_SigmaPML[1][0], 1, &outarray[indxQ1][0], 1,
2537  &outarray[indxQ1][0], 1);
2538  Vmath::Neg(nq, &outarray[indxQ1][0], 1);
2539 
2540  if (m_DispersiveCloak)
2541  {
2542  Vmath::Vvtvp(nq, &m_wp2[0], 1, &physarray[indxE1][0], 1,
2543  &outarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2544  }
2545 
2546  // dP3/dt: Assign \sigma_1 * \sigma_2 * H_z
2547  Vmath::Vmul(nq, &m_SigmaPML[0][0], 1, &m_SigmaPML[1][0], 1,
2548  &outarray[indxP2][0], 1);
2549  Vmath::Vmul(nq, &physarray[indxH2][0], 1, &outarray[indxP2][0], 1,
2550  &outarray[indxP2][0], 1);
2551  }
2552  break;
2553 
2554  default:
2555  break;
2556  }
2557 }
Array< OneD, Array< OneD, NekDouble > > m_SigmaPML
Definition: MMFMaxwell.h:140
Array< OneD, NekDouble > m_wp2
Definition: MMFMaxwell.h:119
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 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::eTransElectric, Nektar::SolverUtils::eTransMagnetic, m_DispersiveCloak, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_PolType, m_SigmaPML, m_wp2, Vmath::Neg(), Vmath::Vadd(), Vmath::Vmul(), Vmath::Vsub(), and Vmath::Vvtvp().

Referenced by DoOdeRhs().

◆ Checkpoint_EDFluxOutput()

void Nektar::MMFMaxwell::Checkpoint_EDFluxOutput ( const int  n,
const NekDouble  time,
const Array< OneD, const Array< OneD, NekDouble >> &  fieldphys 
)
protected

Definition at line 2681 of file MMFMaxwell.cpp.

2684 {
2685  boost::ignore_unused(time);
2686 
2687  int nvar = m_fields.size();
2688  int nq = m_fields[0]->GetTotPoints();
2689  int ncoeffs = m_fields[0]->GetNcoeffs();
2690 
2691  std::string outname = m_sessionName + "EDFlux_" +
2692  boost::lexical_cast<std::string>(n) + ".chk";
2693 
2694  std::vector<std::string> variables(nvar);
2695  variables[0] = "EDFx";
2696  variables[1] = "EDFy";
2697  variables[2] = "EDFz";
2698  for (int i = 3; i < nvar; ++i)
2699  {
2700  variables[i] = m_boundaryConditions->GetVariable(i);
2701  }
2702 
2703  std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2704  for (int i = 0; i < nvar; ++i)
2705  {
2706  fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2707  }
2708 
2709  Array<OneD, NekDouble> tmp(nq);
2710 
2711  // TE: H^3 (E^2 e^1 - E^1 e^2 )
2712  // TM: -E^3 (H^2 e^1 - H^1 e^2 )
2713  for (int k = 0; k < m_spacedim; ++k)
2714  {
2715  Vmath::Vmul(nq, &fieldphys[0][0], 1, &m_movingframes[1][k * nq], 1,
2716  &tmp[0], 1);
2717  Vmath::Vvtvm(nq, &fieldphys[1][0], 1, &m_movingframes[0][k * nq], 1,
2718  &tmp[0], 1, &tmp[0], 1);
2719 
2720  Vmath::Vmul(nq, &fieldphys[2][0], 1, &tmp[0], 1, &tmp[0], 1);
2721 
2723  {
2724  Vmath::Neg(nq, tmp, 1);
2725  }
2726 
2727  m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2728  }
2729 
2730  WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2731 }
int m_spacedim
Spatial dimension (>= expansion dim).
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
std::string m_sessionName
Name of the session.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:186

References Nektar::SolverUtils::eTransMagnetic, Nektar::SolverUtils::EquationSystem::m_boundaryConditions, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_PolType, Nektar::SolverUtils::EquationSystem::m_sessionName, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Neg(), Vmath::Vmul(), Vmath::Vvtvm(), and Nektar::SolverUtils::EquationSystem::WriteFld().

Referenced by v_DoSolve().

◆ Checkpoint_EnergyOutput()

void Nektar::MMFMaxwell::Checkpoint_EnergyOutput ( const int  n,
const NekDouble  time,
const Array< OneD, const Array< OneD, NekDouble >> &  fieldphys 
)
protected

Definition at line 2733 of file MMFMaxwell.cpp.

2736 {
2737  boost::ignore_unused(time);
2738 
2739  int nvar = m_fields.size();
2740  int nq = m_fields[0]->GetTotPoints();
2741  int ncoeffs = m_fields[0]->GetNcoeffs();
2742 
2743  std::string outname = m_sessionName + "Energy_" +
2744  boost::lexical_cast<std::string>(n) + ".chk";
2745 
2746  std::vector<std::string> variables(nvar);
2747  variables[0] = "Energy";
2748  variables[1] = "EnergyFlux";
2749  variables[2] = "Zero";
2750  for (int i = 3; i < nvar; ++i)
2751  {
2752  variables[i] = m_boundaryConditions->GetVariable(i);
2753  }
2754 
2755  std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2756  for (int i = 0; i < nvar; ++i)
2757  {
2758  fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2759  }
2760 
2761  // Energy = 0.5 * ( E^2 + H^2 )
2762  Array<OneD, NekDouble> energy(nq, 0.0);
2763  Array<OneD, NekDouble> totfield(nq);
2764  for (int k = 0; k < m_spacedim; ++k)
2765  {
2766  // totfield = GetIncidentField(k,time);
2767  // Vmath::Vadd(nq, &fieldphys[k][0], 1, &totfield[0], 1, &totfield[0],
2768  // 1);
2769 
2770  Vmath::Vvtvp(nq, &fieldphys[k][0], 1, &fieldphys[k][0], 1, &energy[0],
2771  1, &energy[0], 1);
2772  }
2773  Vmath::Smul(nq, 0.5, energy, 1, energy, 1);
2774  m_fields[0]->FwdTrans(energy, fieldcoeffs[0]);
2775 
2776  // EnergyFlux = F3* sqrt( F1^2 + F2^2 )
2777  Array<OneD, NekDouble> energyflux(nq, 0.0);
2778  Array<OneD, NekDouble> Zero(nq, 0.0);
2779  for (int k = 0; k < 2; ++k)
2780  {
2781  // totfield = GetIncidentField(k,time);
2782  // Vmath::Vadd(nq, &fieldphys[k][0], 1, &totfield[0], 1, &totfield[0],
2783  // 1);
2784 
2785  Vmath::Vvtvp(nq, &fieldphys[k][0], 1, &fieldphys[k][0], 1,
2786  &energyflux[0], 1, &energyflux[0], 1);
2787  }
2788 
2789  Vmath::Vsqrt(nq, energyflux, 1, energyflux, 1);
2790  Vmath::Vmul(nq, &fieldphys[2][0], 1, &energyflux[0], 1, &energyflux[0], 1);
2791 
2792  m_fields[1]->FwdTrans(energyflux, fieldcoeffs[1]);
2793  m_fields[2]->FwdTrans(Zero, fieldcoeffs[2]);
2794 
2795  WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2796 }
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
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
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492

References Nektar::SolverUtils::EquationSystem::m_boundaryConditions, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_sessionName, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Smul(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vvtvp(), Nektar::SolverUtils::EquationSystem::WriteFld(), and Vmath::Zero().

Referenced by v_DoSolve().

◆ Checkpoint_PlotOutput()

void Nektar::MMFMaxwell::Checkpoint_PlotOutput ( const int  n,
const Array< OneD, const Array< OneD, NekDouble >> &  fieldphys 
)
protected

Definition at line 2592 of file MMFMaxwell.cpp.

2594 {
2595  int nvar = m_fields.size();
2596  int nq = m_fields[0]->GetTotPoints();
2597  int ncoeffs = m_fields[0]->GetNcoeffs();
2598 
2599  std::string outname =
2600  m_sessionName + "Plot_" + boost::lexical_cast<std::string>(n) + ".chk";
2601 
2602  std::vector<std::string> variables(nvar);
2603  variables[0] = "Fx";
2604  variables[1] = "Fy";
2605  variables[2] = "Fz";
2606 
2607  std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2608  for (int i = 0; i < nvar; ++i)
2609  {
2610  fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2611  }
2612 
2613  Array<OneD, NekDouble> tmp(nq);
2614  for (int k = 0; k < m_spacedim; ++k)
2615  {
2616  Vmath::Vmul(nq, &fieldphys[0][0], 1, &m_movingframes[0][k * nq], 1,
2617  &tmp[0], 1);
2618  Vmath::Vvtvp(nq, &fieldphys[1][0], 1, &m_movingframes[1][k * nq], 1,
2619  &tmp[0], 1, &tmp[0], 1);
2620 
2621  m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2622  }
2623 
2624  WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2625 }

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

Referenced by v_DoSolve(), and v_SetInitialConditions().

◆ Checkpoint_TotalFieldOutput()

void Nektar::MMFMaxwell::Checkpoint_TotalFieldOutput ( const int  n,
const NekDouble  time,
const Array< OneD, const Array< OneD, NekDouble >> &  fieldphys 
)
protected

Definition at line 2559 of file MMFMaxwell.cpp.

2562 {
2563  int nvar = m_fields.size();
2564  int nq = m_fields[0]->GetTotPoints();
2565  int ncoeffs = m_fields[0]->GetNcoeffs();
2566 
2567  std::string outname =
2568  m_sessionName + "Tot_" + boost::lexical_cast<std::string>(n) + ".chk";
2569 
2570  std::vector<std::string> variables(nvar);
2571  std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2572 
2573  for (int i = 0; i < nvar; ++i)
2574  {
2575  fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2576  }
2577 
2578  Array<OneD, NekDouble> totfield(nq);
2579  for (int i = 0; i < nvar; ++i)
2580  {
2581  totfield = GetIncidentField(i, time);
2582  Vmath::Vadd(nq, fieldphys[i], 1, totfield, 1, totfield, 1);
2583 
2584  m_fields[i]->FwdTrans(totfield, fieldcoeffs[i]);
2585  variables[i] = m_boundaryConditions->GetVariable(i);
2586  }
2587 
2588  WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2589 }
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetIncidentField(const int var, const NekDouble time)
Definition: MMFSystem.cpp:2011

References Nektar::SolverUtils::MMFSystem::GetIncidentField(), Nektar::SolverUtils::EquationSystem::m_boundaryConditions, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_sessionName, Vmath::Vadd(), and Nektar::SolverUtils::EquationSystem::WriteFld().

Referenced by v_DoSolve().

◆ Checkpoint_TotPlotOutput()

void Nektar::MMFMaxwell::Checkpoint_TotPlotOutput ( const int  n,
const NekDouble  time,
const Array< OneD, const Array< OneD, NekDouble >> &  fieldphys 
)
protected

Definition at line 2627 of file MMFMaxwell.cpp.

2630 {
2631  int nvar = m_fields.size();
2632  int nq = m_fields[0]->GetTotPoints();
2633  int ncoeffs = m_fields[0]->GetNcoeffs();
2634 
2635  std::string outname = m_sessionName + "TotPlot_" +
2636  boost::lexical_cast<std::string>(n) + ".chk";
2637 
2638  std::vector<std::string> variables(nvar);
2639  variables[0] = "Frx";
2640  variables[1] = "Fry";
2641  variables[2] = "Frz";
2642  for (int i = 3; i < nvar; ++i)
2643  {
2644  variables[i] = m_boundaryConditions->GetVariable(i);
2645  }
2646 
2647  std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2648  for (int i = 0; i < nvar; ++i)
2649  {
2650  fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2651  }
2652 
2653  Array<OneD, NekDouble> tmp(nq);
2654  Array<OneD, NekDouble> totfield0(nq);
2655  Array<OneD, NekDouble> totfield1(nq);
2656 
2657  totfield0 = GetIncidentField(0, time);
2658  Vmath::Vadd(nq, fieldphys[0], 1, totfield0, 1, totfield0, 1);
2659 
2660  totfield1 = GetIncidentField(1, time);
2661  Vmath::Vadd(nq, fieldphys[1], 1, totfield1, 1, totfield1, 1);
2662 
2663  for (int k = 0; k < m_spacedim; ++k)
2664  {
2665  Vmath::Vmul(nq, &totfield0[0], 1, &m_movingframes[0][k * nq], 1,
2666  &tmp[0], 1);
2667  Vmath::Vvtvp(nq, &totfield1[0], 1, &m_movingframes[1][k * nq], 1,
2668  &tmp[0], 1, &tmp[0], 1);
2669 
2670  m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2671  }
2672 
2673  for (int j = 3; j < nvar; ++j)
2674  {
2675  m_fields[j]->FwdTrans(fieldphys[j], fieldcoeffs[j]);
2676  }
2677 
2678  WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2679 }

References Nektar::SolverUtils::MMFSystem::GetIncidentField(), Nektar::SolverUtils::EquationSystem::m_boundaryConditions, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::EquationSystem::m_sessionName, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vadd(), Vmath::Vmul(), Vmath::Vvtvp(), and Nektar::SolverUtils::EquationSystem::WriteFld().

Referenced by v_DoSolve().

◆ ComputeEnergyDensity()

NekDouble Nektar::MMFMaxwell::ComputeEnergyDensity ( Array< OneD, Array< OneD, NekDouble >> &  fields)
protected

Definition at line 2182 of file MMFMaxwell.cpp.

2184 {
2185  int nq = GetTotPoints();
2186  NekDouble energy;
2187 
2188  Array<OneD, NekDouble> tmp(nq, 0.0);
2189 
2190  for (int i = 0; i < 3; ++i)
2191  {
2192  Vmath::Vvtvp(nq, &fields[i][0], 1, &fields[i][0], 1, &tmp[0], 1,
2193  &tmp[0], 1);
2194  }
2195 
2196  energy = 0.5 * (m_fields[0]->Integral(tmp));
2197  return energy;
2198 }
SOLVER_UTILS_EXPORT int GetTotPoints()
double NekDouble

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

Referenced by v_DoSolve().

◆ ComputeMaterialMicroWaveCloak()

void Nektar::MMFMaxwell::ComputeMaterialMicroWaveCloak ( const Array< OneD, const NekDouble > &  radvec,
Array< OneD, Array< OneD, NekDouble >> &  epsvec,
Array< OneD, Array< OneD, NekDouble >> &  muvec 
)
protected

Definition at line 2368 of file MMFMaxwell.cpp.

2372 {
2373  int nq = GetNpoints();
2374 
2375  NekDouble m_b = 2.67;
2376  NekDouble m_a = 1.33;
2377  NekDouble m_adel;
2378 
2379  m_adel = m_a - m_Cloakraddelta;
2380 
2381  Array<OneD, NekDouble> Cloakregion(nq, 0.0);
2382  NekDouble ExactCloakArea = m_pi * (m_b * m_b - m_a * m_a);
2383  m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0, 0.0, Cloakregion);
2384 
2385  if (m_ElemtGroup0 > 0)
2386  {
2387  Array<OneD, NekDouble> Vacregion(nq, 0.0);
2388  m_fields[0]->GenerateElementVector(m_ElemtGroup0, 1.0, 0.0, Vacregion);
2389 
2390  Vmath::Vsub(nq, Cloakregion, 1, Vacregion, 1, Cloakregion, 1);
2391  }
2392 
2393  ExactCloakArea = ExactCloakArea - (m_fields[0]->Integral(Cloakregion));
2394  std::cout << "*** Error of Cloakregion area = " << ExactCloakArea
2395  << std::endl;
2396 
2397  epsvec[0] = Array<OneD, NekDouble>(nq, 1.0);
2398  epsvec[1] = Array<OneD, NekDouble>(nq, 1.0);
2399 
2400  muvec[0] = Array<OneD, NekDouble>(nq, 1.0);
2401  muvec[1] = Array<OneD, NekDouble>(nq, 1.0);
2402  for (int i = 0; i < nq; ++i)
2403  {
2404  if (Cloakregion[i] > 0)
2405  {
2406  // relrad = m_a +
2407  // (m_b-m_a)*(radvec[i]-Cloakradmin)/(Cloakradmax-Cloakradmin);
2408  // ratio = (relrad - m_a + m_Cloakraddelta)/relrad;
2409 
2410  epsvec[0][i] = radvec[i] / (radvec[i] - m_adel);
2411  epsvec[1][i] = (radvec[i] - m_adel) / radvec[i];
2412  muvec[2][i] = (m_b / (m_b - m_adel)) * (m_b / (m_b - m_adel)) *
2413  (radvec[i] - m_adel) / radvec[i];
2414 
2415  muvec[0][i] = epsvec[0][i];
2416  muvec[1][i] = epsvec[1][i];
2417  epsvec[2][i] = muvec[2][i];
2418  }
2419  }
2420 }
NekDouble m_Cloakraddelta
Definition: MMFMaxwell.h:116
SOLVER_UTILS_EXPORT int GetNpoints()

References Nektar::SolverUtils::EquationSystem::GetNpoints(), m_Cloakraddelta, m_ElemtGroup0, m_ElemtGroup1, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_pi, and Vmath::Vsub().

Referenced by v_InitObject().

◆ ComputeMaterialOpticalCloak()

void Nektar::MMFMaxwell::ComputeMaterialOpticalCloak ( const Array< OneD, const NekDouble > &  radvec,
Array< OneD, Array< OneD, NekDouble >> &  epsvec,
Array< OneD, Array< OneD, NekDouble >> &  muvec,
const bool  Dispersion = false 
)
protected

Definition at line 2295 of file MMFMaxwell.cpp.

2299 {
2300  boost::ignore_unused(muvec, Dispersion);
2301 
2302  int nq = GetNpoints();
2303 
2304  // Cloaking metamaterial
2305  // \varepsilon_\theta = (b/(b-a))^2
2306  // \varepsilon_r = (b/(b-a))^2 ((r-a)/r)^2
2307  // \mu_z = 1.0m_CloakingOn
2308 
2309  NekDouble m_b = 1.0 / 0.314;
2310  NekDouble m_a = 1.0;
2311 
2312  NekDouble m_adel = m_a - m_Cloakraddelta;
2313 
2314  NekDouble boveradel = m_b / (m_b - m_adel);
2315 
2316  Array<OneD, NekDouble> Cloakregion(nq, 0.0);
2317 
2318  NekDouble la = m_MMFfactors[0];
2319  NekDouble lb = m_MMFfactors[1];
2320 
2321  NekDouble ExactCloakArea;
2322  if (fabs(la * lb - 1.0) < 0.00001)
2323  {
2324  ExactCloakArea = m_pi * (m_b * m_b - m_a * m_a);
2325  }
2326 
2327  else
2328  {
2329  ExactCloakArea = m_pi * (3.0 * lb * 3.0 * la - lb * la);
2330  }
2331 
2332  m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0, 0.0, Cloakregion);
2333 
2334  ExactCloakArea = ExactCloakArea - (m_fields[0]->Integral(Cloakregion));
2335  std::cout << "*** Error of Cloakregion area = " << ExactCloakArea
2336  << std::endl;
2337 
2338  NekDouble ratio;
2339 
2340  epsvec[0] = Array<OneD, NekDouble>(nq, 1.0);
2341  epsvec[1] = Array<OneD, NekDouble>(nq, 1.0);
2342  m_wp2 = Array<OneD, NekDouble>(nq, 0.0);
2343  for (int i = 0; i < nq; ++i)
2344  {
2345  if (Cloakregion[i] > 0)
2346  {
2347  // relrad = m_a +
2348  // (m_b-m_adel)*(radvec[i]-m_adel)/(Cloakradmax-m_adel);
2349  ratio = (radvec[i] - m_adel) / radvec[i];
2350 
2351  epsvec[0][i] = boveradel * boveradel;
2352  if (m_DispersiveCloak)
2353  {
2354  epsvec[1][i] = 1.0;
2355  m_wp2[i] = m_Incfreq * m_Incfreq *
2356  (1.0 - boveradel * boveradel * ratio * ratio) +
2357  m_wp2Tol;
2358  }
2359 
2360  else
2361  {
2362  epsvec[1][i] = boveradel * boveradel * (ratio * ratio);
2363  }
2364  }
2365  }
2366 }
NekDouble m_wp2Tol
Definition: MMFMaxwell.h:118
Array< OneD, NekDouble > m_MMFfactors
Definition: MMFSystem.h:160

References Nektar::SolverUtils::EquationSystem::GetNpoints(), m_Cloakraddelta, m_DispersiveCloak, m_ElemtGroup1, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_Incfreq, Nektar::SolverUtils::MMFSystem::m_MMFfactors, Nektar::SolverUtils::MMFSystem::m_pi, m_wp2, and m_wp2Tol.

Referenced by v_InitObject().

◆ ComputeMaterialVector()

void Nektar::MMFMaxwell::ComputeMaterialVector ( Array< OneD, Array< OneD, NekDouble >> &  epsvec,
Array< OneD, Array< OneD, NekDouble >> &  muvec 
)
protected

Definition at line 2200 of file MMFMaxwell.cpp.

2203 {
2204  switch (m_TestMaxwellType)
2205  {
2207  {
2208  m_fields[0]->GenerateElementVector(m_ElemtGroup1, m_varepsilon[0],
2209  m_varepsilon[1], epsvec[0]);
2210  m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0, 1.0,
2211  muvec[0]);
2212  }
2213  break;
2214 
2220  {
2221  switch (m_PolType)
2222  {
2224  {
2225  m_fields[0]->GenerateElementVector(m_ElemtGroup1, m_mu[0],
2226  1.0, muvec[0]);
2227  m_fields[0]->GenerateElementVector(m_ElemtGroup1, m_mu[1],
2228  1.0, muvec[1]);
2229  m_fields[0]->GenerateElementVector(
2230  m_ElemtGroup1, m_varepsilon[2], 1.0, epsvec[2]);
2231 
2232  // // ONLY FOR VARIABLE ANISOTROPY TEST
2233  // int nq = GetTotPoints();
2234 
2235  // Array<OneD, NekDouble> tmpIN(nq);
2236 
2237  // m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0,
2238  // 0.0, tmpIN);
2239 
2240  // Array<OneD, NekDouble> x0(nq);
2241  // Array<OneD, NekDouble> x1(nq);
2242  // Array<OneD, NekDouble> x2(nq);
2243 
2244  // m_fields[0]->GetCoords(x0,x1,x2);
2245 
2246  // for (int i=0; i<nq; i++)
2247  // {
2248  // muvec[1][i] = tmpIN[i]*(1.0 - 2*x0[i]*x0[i]) +
2249  // (1.0-tmpIN[i]);
2250  // }
2251  }
2252  break;
2253 
2255  {
2256  m_fields[0]->GenerateElementVector(
2257  m_ElemtGroup1, m_varepsilon[0], 1.0, epsvec[0]);
2258  m_fields[0]->GenerateElementVector(
2259  m_ElemtGroup1, m_varepsilon[1], 1.0, epsvec[1]);
2260  m_fields[0]->GenerateElementVector(m_ElemtGroup1, m_mu[2],
2261  1.0, muvec[2]);
2262 
2263  // // ONLY FOR VARIABLE ANISOTROPY TEST
2264  // int nq = GetTotPoints();
2265 
2266  // Array<OneD, NekDouble> tmpIN(nq);
2267 
2268  // m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0,
2269  // 0.0, tmpIN);
2270 
2271  // Array<OneD, NekDouble> x0(nq);
2272  // Array<OneD, NekDouble> x1(nq);
2273  // Array<OneD, NekDouble> x2(nq);
2274 
2275  // m_fields[0]->GetCoords(x0,x1,x2);
2276 
2277  // for (int i=0; i<nq; i++)
2278  // {
2279  // epsvec[1][i] = tmpIN[i]*(1.0 - 2*x0[i]*x0[i]) +
2280  // (1.0-tmpIN[i]);
2281  // }
2282  }
2283  break;
2284 
2285  default:
2286  break; // Pol
2287  }
2288  }
2289 
2290  default:
2291  break; // TestType
2292  }
2293 }
Array< OneD, NekDouble > m_mu
Definition: MMFMaxwell.h:135
Array< OneD, NekDouble > m_varepsilon
Definition: MMFMaxwell.h:134

References Nektar::SolverUtils::eMaxwell1D, Nektar::SolverUtils::eScatField2D, Nektar::SolverUtils::eTestMaxwell2DPEC, Nektar::SolverUtils::eTestMaxwell2DPECAVGFLUX, Nektar::SolverUtils::eTestMaxwell2DPMC, Nektar::SolverUtils::eTotField2D, Nektar::SolverUtils::eTransElectric, Nektar::SolverUtils::eTransMagnetic, m_ElemtGroup1, Nektar::SolverUtils::EquationSystem::m_fields, m_mu, Nektar::SolverUtils::MMFSystem::m_PolType, Nektar::SolverUtils::MMFSystem::m_TestMaxwellType, and m_varepsilon.

Referenced by v_InitObject().

◆ ComputeRadCloak()

Array< OneD, NekDouble > Nektar::MMFMaxwell::ComputeRadCloak ( const int  Nlayer = 5)
protected

Definition at line 2915 of file MMFMaxwell.cpp.

2916 {
2917  int nq = GetNpoints();
2918 
2919  NekDouble m_b = 1.0 / 0.314;
2920  NekDouble m_a = 1.0;
2921 
2922  Array<OneD, int> Layer(CloakNlayer);
2923 
2924  NekDouble la = m_MMFfactors[0];
2925  NekDouble lb = m_MMFfactors[1];
2926  if (CloakNlayer == 8)
2927  {
2928  // Circular 8 layer
2929  if (fabs(la * lb - 1.0) < 0.0001)
2930  {
2931  Layer[0] = 119;
2932  Layer[1] = 239;
2933  Layer[2] = 323;
2934  Layer[3] = 459;
2935  Layer[4] = 575;
2936  Layer[5] = 727;
2937  Layer[6] = 891;
2938  Layer[7] = 1059;
2939  }
2940 
2941  // Ellipsoid 8 layer
2942  else
2943  {
2944  Layer[0] = 115;
2945  Layer[1] = 219;
2946  Layer[2] = 335;
2947  Layer[3] = 459;
2948  Layer[4] = 607;
2949  Layer[5] = 767;
2950  Layer[6] = 951;
2951  Layer[7] = 1155;
2952  }
2953  }
2954 
2955  if (CloakNlayer == 16)
2956  {
2957  // Circular 16 layer
2958  if (fabs(la * lb - 1.0) < 0.0001)
2959  {
2960  Layer[0] = 43;
2961  Layer[1] = 87;
2962  Layer[2] = 135;
2963  Layer[3] = 187;
2964  Layer[4] = 239;
2965  Layer[5] = 295;
2966  Layer[6] = 355;
2967  Layer[7] = 415;
2968  Layer[8] = 479;
2969  Layer[9] = 555;
2970  Layer[10] = 639;
2971  Layer[11] = 727;
2972  Layer[12] = 823;
2973  Layer[13] = 927;
2974  Layer[14] = 1039;
2975  Layer[15] = 1159;
2976  }
2977 
2978  // Ellipsoid 8 layer
2979  else
2980  {
2981  Layer[0] = 135;
2982  Layer[1] = 259;
2983  Layer[2] = 387;
2984  Layer[3] = 523;
2985  Layer[4] = 667;
2986  Layer[5] = 803;
2987  Layer[6] = 939;
2988  Layer[7] = 1083;
2989  Layer[8] = 1235;
2990  Layer[9] = 1387;
2991  Layer[10] = 1539;
2992  Layer[11] = 1699;
2993  Layer[12] = 1867;
2994  Layer[13] = 2035;
2995  Layer[14] = 2203;
2996  Layer[15] = 2379;
2997  }
2998  }
2999 
3000  Array<OneD, NekDouble> x0(nq);
3001  Array<OneD, NekDouble> x1(nq);
3002  Array<OneD, NekDouble> x2(nq);
3003 
3004  m_fields[0]->GetCoords(x0, x1, x2);
3005 
3006  Array<OneD, NekDouble> Formertmp(nq, 0.0);
3007  Array<OneD, NekDouble> Currenttmp(nq, 0.0);
3008  Array<OneD, NekDouble> Cloakregion(nq, 0.0);
3009 
3010  m_fields[0]->GenerateElementVector(Layer[0], 1.0, 0.0, Currenttmp);
3011  Vmath::Vcopy(nq, Currenttmp, 1, Cloakregion, 1);
3012  for (int i = 1; i < CloakNlayer; ++i)
3013  {
3014  m_fields[0]->GenerateElementVector(Layer[i], 1.0, 0.0, Currenttmp);
3015  m_fields[0]->GenerateElementVector(Layer[i - 1], 1.0, 0.0, Formertmp);
3016 
3017  Vmath::Vsub(nq, Currenttmp, 1, Formertmp, 1, Currenttmp, 1);
3018 
3019  Vmath::Svtvp(nq, 1.0 * (i + 1), Currenttmp, 1, Cloakregion, 1,
3020  Cloakregion, 1);
3021  }
3022 
3023  Array<OneD, NekDouble> radvec(nq);
3024 
3025  for (int i = 0; i < nq; ++i)
3026  {
3027  switch (m_MMFdir)
3028  {
3030  {
3031  if ((Cloakregion[i] > 0) && (CloakNlayer > 0))
3032  {
3033  radvec[i] = 1.0 + (m_b - m_a) / CloakNlayer *
3034  (Cloakregion[i] - 0.5);
3035  }
3036 
3037  else
3038  {
3039  radvec[i] =
3040  sqrt(x0[i] * x0[i] / la / la + x1[i] * x1[i] / lb / lb);
3041  }
3042  }
3043  break;
3044 
3046  {
3047  radvec[i] = sqrt(2.0 * x0[i] * x0[i] +
3048  x1[i] * x1[i] * x1[i] * x1[i] + x1[i] * x1[i]);
3049  }
3050  break;
3051 
3053  {
3054  radvec[i] = sqrt(3.0 * x0[i] * x0[i] +
3055  x1[i] * x1[i] * x1[i] * x1[i] - x1[i] * x1[i]);
3056  }
3057  break;
3058 
3059  default:
3060  break;
3061  }
3062  }
3063 
3064  return radvec;
3065 }
SpatialDomains::GeomMMF m_MMFdir
Definition: MMFSystem.h:222
@ eTangentIrregular
Circular around the centre of domain.
@ eTangentCircular
Circular around the centre of domain.
@ eTangentNonconvex
Circular around the centre of domain.
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
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291

References Nektar::SpatialDomains::eTangentCircular, Nektar::SpatialDomains::eTangentIrregular, Nektar::SpatialDomains::eTangentNonconvex, Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_MMFdir, Nektar::SolverUtils::MMFSystem::m_MMFfactors, tinysimd::sqrt(), Vmath::Svtvp(), Vmath::Vcopy(), and Vmath::Vsub().

Referenced by v_InitObject().

◆ ComputeSurfaceCurrent()

Array< OneD, NekDouble > Nektar::MMFMaxwell::ComputeSurfaceCurrent ( const int  time,
const Array< OneD, const Array< OneD, NekDouble >> &  fields 
)
protected

Definition at line 1981 of file MMFMaxwell.cpp.

1983 {
1984  int nq = m_fields[0]->GetTotPoints();
1985  int nTraceNumPoints = GetTraceTotPoints();
1986 
1987  Array<OneD, NekDouble> outfield(nTraceNumPoints, 0.0);
1988 
1989  switch (m_PolType)
1990  {
1991  // (n \times H^r)_z = H1r (n \times e^1)_z + H2r (n \times e^2)_z,
1993  {
1994  Array<OneD, NekDouble> tmp(nq);
1995  Array<OneD, NekDouble> tmpFwd(nTraceNumPoints);
1996 
1997  for (int i = 0; i < 2; i++)
1998  {
1999  tmp = GetIncidentField(i, time);
2000  Vmath::Vadd(nq, fields[i], 1, tmp, 1, tmp, 1);
2001 
2002  m_fields[0]->ExtractTracePhys(tmp, tmpFwd);
2003  Vmath::Vvtvp(nTraceNumPoints, &m_ntimesMFFwd[i][2][0], 1,
2004  &tmpFwd[0], 1, &outfield[0], 1, &outfield[0], 1);
2005  }
2006  }
2007  break;
2008 
2010  {
2011  Array<OneD, NekDouble> tmp(nq);
2012 
2013  tmp = GetIncidentField(2, time);
2014  Vmath::Vadd(nq, fields[2], 1, tmp, 1, tmp, 1);
2015  m_fields[0]->ExtractTracePhys(tmp, outfield);
2016  }
2017  break;
2018 
2019  default:
2020  break;
2021  }
2022 
2023  return outfield;
2024 }
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFFwd
Definition: MMFSystem.h:202

References Nektar::SolverUtils::eTransElectric, Nektar::SolverUtils::eTransMagnetic, Nektar::SolverUtils::MMFSystem::GetIncidentField(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_ntimesMFFwd, Nektar::SolverUtils::MMFSystem::m_PolType, Vmath::Vadd(), and Vmath::Vvtvp().

Referenced by Printout_SurfaceCurrent().

◆ create()

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

Creates an instance of this class.

Definition at line 83 of file MMFMaxwell.h.

86  {
89  p->InitObject();
90  return p;
91  }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.

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

◆ DoOdeProjection()

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

Compute the projection.

Compute the projection for the linear advection equation.

Parameters
inarrayGiven fields.
outarrayCalculated solution.
timeTime.

Definition at line 1318 of file MMFMaxwell.cpp.

1321 {
1322  boost::ignore_unused(time);
1323 
1324  int var = inarray.size();
1325 
1326  // SetBoundaryConditions(time);
1327 
1328  int nq = GetNpoints();
1329  for (int i = 0; i < var; ++i)
1330  {
1331  Vmath::Vcopy(nq, inarray[i], 1, outarray[i], 1);
1332  }
1333 }

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

Referenced by v_InitObject().

◆ DoOdeRhs()

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

Compute the RHS.

Compute the right-hand side for the linear advection equation.

Parameters
inarrayGiven fields.
outarrayCalculated solution.
timeTime.

Definition at line 887 of file MMFMaxwell.cpp.

890 {
891  int i;
892  int nvar = inarray.size();
893  int ncoeffs = GetNcoeffs();
894  int nq = GetTotPoints();
895 
896  Array<OneD, Array<OneD, NekDouble>> physarray(nvar);
897  Array<OneD, Array<OneD, NekDouble>> modarray(nvar);
898  for (i = 0; i < nvar; ++i)
899  {
900  physarray[i] = Array<OneD, NekDouble>(nq);
901  modarray[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
902 
903  Vmath::Vcopy(nq, &inarray[i][0], 1, &physarray[i][0], 1);
904  }
905 
906  for (i = 0; i < nvar; i++)
907  {
908  m_fields[i]->SetPhysState(true);
909  }
910 
911  // Compute Curl
912  switch (m_TestMaxwellType)
913  {
916 
917  {
918 
919  // Imaginary part is computed the same as Real part
920  Array<OneD, Array<OneD, NekDouble>> tmpin(3);
921  Array<OneD, Array<OneD, NekDouble>> tmpout(3);
922 
923  for (int i = 0; i < 3; ++i)
924  {
925  tmpin[i] = Array<OneD, NekDouble>(nq);
926  tmpout[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
927 
928  Vmath::Vcopy(nq, &physarray[i][0], 1, &tmpin[i][0], 1);
929  }
930 
931  WeakDGMaxwellDirDeriv(tmpin, tmpout, time);
932  AddGreenDerivCompensate(tmpin, tmpout);
933 
934  for (int i = 0; i < 3; ++i)
935  {
936  // For E and H
937  Vmath::Vcopy(ncoeffs, &tmpout[i][0], 1, &modarray[i][0], 1);
938  }
939  }
940  break;
941 
942  default:
943  {
944 
945  WeakDGMaxwellDirDeriv(physarray, modarray, time);
946  AddGreenDerivCompensate(physarray, modarray);
947  }
948  break;
949  }
950 
951  for (i = 0; i < nvar; ++i)
952  {
953  m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
954  m_fields[i]->BwdTrans(modarray[i], outarray[i]);
955  }
956 
958  {
959  Array<OneD, NekDouble> F(nq);
960  for (int j = 0; j < 2; ++j)
961  {
962  F = TestMaxwellSphere(time, m_freq, 3 + j);
963  Vmath::Vadd(nq, &F[0], 1, &outarray[j][0], 1, &outarray[j][0], 1);
964  }
965  }
966 
967  // Add Absorbing Boundary Conditions
968  if (m_AddPML > 0)
969  {
970  AddPML(physarray, outarray);
971  }
972 
973  // Add dedt component
974  if (m_AddRotation)
975  {
976  AddCoriolis(physarray, outarray);
977  AdddedtMaxwell(physarray, outarray);
978  }
979 
980  // Divide it by varepsilon or mu
981  Array<OneD, NekDouble> dFdt(nq, 0.0);
982  switch (m_TestMaxwellType)
983  {
985  {
986  Vmath::Vdiv(nq, outarray[0], 1, m_epsvec[0], 1, outarray[0], 1);
987  Vmath::Vdiv(nq, outarray[1], 1, m_muvec[0], 1, outarray[1], 1);
988  }
989  break;
990 
991  // TO BE CHANGED
993  {
994  Array<OneD, NekDouble> Hxdt(nq, 0.0);
995  Array<OneD, NekDouble> Hydt(nq, 0.0);
996 
997  Hxdt = TestMaxwell2DPEC(time, 10, m_PolType);
998  Hydt = TestMaxwell2DPEC(time, 11, m_PolType);
999 
1000  Array<OneD, NekDouble> x0(nq);
1001  Array<OneD, NekDouble> x1(nq);
1002  Array<OneD, NekDouble> x2(nq);
1003 
1004  m_fields[0]->GetCoords(x0, x1, x2);
1005 
1006  NekDouble theta, tmpx, tmpy;
1007  NekDouble uxx, uxy, uyy, detu, uti, uri;
1008  Array<OneD, NekDouble> utvec(nq, 1.0);
1009  Array<OneD, NekDouble> urvec(nq, 1.0);
1010  Array<OneD, NekDouble> tmpIN(nq);
1011 
1012  // Case I: ut = 4.0, ur = 0.5
1013  // NekDouble ut=4.0;
1014  // NekDouble ur=0.5;
1015 
1016  // m_fields[0]->GenerateElementVector(m_ElemtGroup1, ut, 1.0,
1017  // utvec);
1018  // m_fields[0]->GenerateElementVector(m_ElemtGroup1, ur, 1.0,
1019  // urvec);
1020 
1021  // Case II: ut = 0.5, ur = 1 - 2x^2
1022  NekDouble ut = 0.5;
1023  m_fields[0]->GenerateElementVector(m_ElemtGroup1, ut, 1.0, utvec);
1024 
1025  m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0, 0.0, tmpIN);
1026 
1027  for (int i = 0; i < nq; i++)
1028  {
1029  urvec[i] =
1030  tmpIN[i] * (1.0 - 2 * x0[i] * x0[i]) + (1.0 - tmpIN[i]);
1031  }
1032 
1033  for (int i = 0; i < nq; ++i)
1034  {
1035  theta = atan2((x1[i] + 2.0), (x0[i] + 2.0));
1036 
1037  uti = utvec[i];
1038  uri = urvec[i];
1039 
1040  uxx = uti * cos(theta) * cos(theta) +
1041  uri * sin(theta) * sin(theta);
1042  uyy = uti * sin(theta) * sin(theta) +
1043  uri * cos(theta) * cos(theta);
1044  uxy = (uti - uri) * cos(theta) * sin(theta);
1045 
1046  detu = uxx * uyy - uxy * uxy;
1047 
1048  tmpx = outarray[0][i] + (1.0 - uxx) * Hxdt[i] - uxy * Hydt[i];
1049  tmpy = outarray[1][i] - uxy * Hxdt[i] + (1.0 - uyy) * Hydt[i];
1050 
1051  outarray[0][i] = (1 / detu) * (uyy * tmpx - uxy * tmpy);
1052  outarray[1][i] = (1 / detu) * (-uxy * tmpx + uxx * tmpy);
1053  }
1054  }
1055  break;
1056 
1060  {
1061  switch (m_PolType)
1062  {
1064  {
1066  {
1067  dFdt = TestMaxwell2DPEC(time, 10, m_PolType);
1068  Vmath::Vvtvp(nq, m_negmuvecminus1[0], 1, dFdt, 1,
1069  outarray[0], 1, outarray[0], 1);
1070 
1071  dFdt = TestMaxwell2DPEC(time, 11, m_PolType);
1072  Vmath::Vvtvp(nq, m_negmuvecminus1[1], 1, dFdt, 1,
1073  outarray[1], 1, outarray[1], 1);
1074 
1075  dFdt = TestMaxwell2DPEC(time, 12, m_PolType);
1076  Vmath::Vvtvp(nq, m_negepsvecminus1[2], 1, dFdt, 1,
1077  outarray[2], 1, outarray[2], 1);
1078  }
1079 
1081  {
1082  dFdt = GetIncidentField(10, time);
1083  Vmath::Vvtvp(nq, m_negmuvecminus1[0], 1, dFdt, 1,
1084  outarray[0], 1, outarray[0], 1);
1085 
1086  dFdt = GetIncidentField(11, time);
1087  Vmath::Vvtvp(nq, m_negmuvecminus1[1], 1, dFdt, 1,
1088  outarray[1], 1, outarray[1], 1);
1089 
1090  dFdt = GetIncidentField(12, time);
1091  Vmath::Vvtvp(nq, m_negepsvecminus1[2], 1, dFdt, 1,
1092  outarray[2], 1, outarray[2], 1);
1093  }
1094 
1095  Vmath::Vdiv(nq, outarray[0], 1, m_muvec[0], 1, outarray[0],
1096  1);
1097  Vmath::Vdiv(nq, outarray[1], 1, m_muvec[1], 1, outarray[1],
1098  1);
1099  Vmath::Vdiv(nq, outarray[2], 1, m_epsvec[2], 1, outarray[2],
1100  1);
1101  }
1102  break;
1103 
1105  {
1107  {
1108  // (I - \mu^i) d F^{inc} / dt
1109  dFdt = TestMaxwell2DPEC(time, 10, m_PolType);
1110  Vmath::Vvtvp(nq, m_negepsvecminus1[0], 1, dFdt, 1,
1111  outarray[0], 1, outarray[0], 1);
1112 
1113  dFdt = TestMaxwell2DPEC(time, 11, m_PolType);
1114  Vmath::Vvtvp(nq, m_negepsvecminus1[1], 1, dFdt, 1,
1115  outarray[1], 1, outarray[1], 1);
1116 
1117  dFdt = TestMaxwell2DPEC(time, 12, m_PolType);
1118  Vmath::Vvtvp(nq, m_negmuvecminus1[2], 1, dFdt, 1,
1119  outarray[2], 1, outarray[2], 1);
1120  }
1121 
1123  {
1124  dFdt = GetIncidentField(10, time);
1125  Vmath::Vvtvp(nq, m_negepsvecminus1[0], 1, dFdt, 1,
1126  outarray[0], 1, outarray[0], 1);
1127 
1128  // Add - wp^2 \int E_2^{inc}
1129  dFdt = GetIncidentField(21, time);
1130  if (m_DispersiveCloak)
1131  {
1132  Vmath::Vmul(nq, m_wp2, 1, dFdt, 1, dFdt, 1);
1133  Vmath::Vsub(nq, outarray[1], 1, dFdt, 1,
1134  outarray[1], 1);
1135  }
1136 
1137  else
1138  {
1139  Vmath::Vvtvp(nq, m_negepsvecminus1[1], 1, dFdt, 1,
1140  outarray[1], 1, outarray[1], 1);
1141  }
1142 
1143  dFdt = GetIncidentField(12, time);
1144  Vmath::Vvtvp(nq, m_negmuvecminus1[2], 1, dFdt, 1,
1145  outarray[2], 1, outarray[2], 1);
1146  }
1147 
1148  Vmath::Vdiv(nq, outarray[0], 1, m_epsvec[0], 1, outarray[0],
1149  1);
1150  Vmath::Vdiv(nq, outarray[1], 1, m_epsvec[1], 1, outarray[1],
1151  1);
1152  Vmath::Vdiv(nq, outarray[2], 1, m_muvec[2], 1, outarray[2],
1153  1);
1154  }
1155  break;
1156 
1157  default:
1158  break;
1159  }
1160  }
1161  break; // case SolverUtils::eTestMaxwell2DPEC:
1162 
1163  default:
1164  break;
1165 
1166  } // switch(m_TestMaxwellType)
1167 }
Array< OneD, NekDouble > TestMaxwellSphere(const NekDouble time, const NekDouble omega, unsigned int field)
Array< OneD, NekDouble > TestMaxwell2DPEC(const NekDouble time, unsigned int field, const SolverUtils::PolType Polarization)
void AddPML(const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
void WeakDGMaxwellDirDeriv(const Array< OneD, const Array< OneD, NekDouble >> &InField, Array< OneD, Array< OneD, NekDouble >> &OutField, const NekDouble time=0.0)
Calculate weak DG advection in the form .
void AddGreenDerivCompensate(const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
NekDouble m_freq
Definition: MMFMaxwell.h:131
void AddCoriolis(Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
SOLVER_UTILS_EXPORT int GetNcoeffs()
Array< OneD, Array< OneD, NekDouble > > m_negepsvecminus1
Definition: MMFSystem.h:215
SOLVER_UTILS_EXPORT void AdddedtMaxwell(const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: MMFSystem.cpp:1370
Array< OneD, Array< OneD, NekDouble > > m_negmuvecminus1
Definition: MMFSystem.h:216
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 AddCoriolis(), Nektar::SolverUtils::MMFSystem::AdddedtMaxwell(), AddGreenDerivCompensate(), AddPML(), Nektar::SolverUtils::eMaxwell1D, Nektar::SolverUtils::eMaxwellSphere, Nektar::SolverUtils::eScatField2D, Nektar::SolverUtils::eTestMaxwell2DPEC, Nektar::SolverUtils::eTestMaxwell2DPECAVGFLUX, Nektar::SolverUtils::eTotField2D, Nektar::SolverUtils::eTransElectric, Nektar::SolverUtils::eTransMagnetic, Nektar::SolverUtils::MMFSystem::GetIncidentField(), Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_AddPML, m_AddRotation, m_DispersiveCloak, m_ElemtGroup1, Nektar::SolverUtils::MMFSystem::m_epsvec, Nektar::SolverUtils::EquationSystem::m_fields, m_freq, Nektar::SolverUtils::MMFSystem::m_muvec, Nektar::SolverUtils::MMFSystem::m_negepsvecminus1, Nektar::SolverUtils::MMFSystem::m_negmuvecminus1, Nektar::SolverUtils::MMFSystem::m_PolType, Nektar::SolverUtils::MMFSystem::m_TestMaxwellType, m_wp2, TestMaxwell2DPEC(), TestMaxwellSphere(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsub(), Vmath::Vvtvp(), and WeakDGMaxwellDirDeriv().

Referenced by v_InitObject().

◆ EvaluateCoriolis()

Array< OneD, NekDouble > Nektar::MMFMaxwell::EvaluateCoriolis ( )
protected

Definition at line 2835 of file MMFMaxwell.cpp.

2836 {
2837  int nq = GetTotPoints();
2838 
2839  NekDouble m_Omega = 1.5486 * 0.000001;
2840 
2841  NekDouble x0j, x1j, x2j;
2842  NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2843 
2844  Array<OneD, NekDouble> x(nq);
2845  Array<OneD, NekDouble> y(nq);
2846  Array<OneD, NekDouble> z(nq);
2847 
2848  m_fields[0]->GetCoords(x, y, z);
2849 
2850  Array<OneD, NekDouble> outarray(nq, 0.0);
2851  for (int j = 0; j < nq; ++j)
2852  {
2853  x0j = x[j];
2854  x1j = y[j];
2855  x2j = z[j];
2856 
2857  CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2858  cos_theta);
2859 
2860  outarray[j] = 2.0 * m_Omega * sin_theta;
2861  }
2862 
2863  return outarray;
2864 }
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795

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

Referenced by v_InitObject().

◆ GaussianPulse()

Array< OneD, NekDouble > Nektar::MMFMaxwell::GaussianPulse ( const NekDouble  time,
const NekDouble  Psx,
const NekDouble  Psy,
const NekDouble  Psz,
const NekDouble  Gaussianradius 
)
protected

Definition at line 2798 of file MMFMaxwell.cpp.

2803 {
2804  int nq = m_fields[0]->GetTotPoints();
2805  int ncoeffs = m_fields[0]->GetNcoeffs();
2806 
2807  Array<OneD, NekDouble> x(nq);
2808  Array<OneD, NekDouble> y(nq);
2809  Array<OneD, NekDouble> z(nq);
2810 
2811  m_fields[0]->GetCoords(x, y, z);
2812 
2813  Array<OneD, NekDouble> outarray(nq, 0.0);
2814  Array<OneD, NekDouble> tmpc(ncoeffs);
2815  NekDouble rad;
2816 
2817  NekDouble SFradius = m_PSduration * 0.1;
2818  NekDouble SmoothFactor =
2819  1.0 / (1.0 + exp(-0.5 * (time - m_PSduration) / SFradius));
2820 
2821  for (int j = 0; j < nq; ++j)
2822  {
2823  rad = sqrt((x[j] - Psx) * (x[j] - Psx) + (y[j] - Psy) * (y[j] - Psy) +
2824  (z[j] - Psz) * (z[j] - Psz));
2825  outarray[j] = SmoothFactor * exp(-1.0 * (rad / Gaussianradius) *
2826  (rad / Gaussianradius));
2827  }
2828 
2829  m_fields[0]->FwdTransLocalElmt(outarray, tmpc);
2830  m_fields[0]->BwdTrans(tmpc, outarray);
2831 
2832  return outarray;
2833 }
NekDouble m_PSduration
Definition: MMFMaxwell.h:123
static NekDouble rad(NekDouble x, NekDouble y)
Definition: Interpreter.cpp:86

References Nektar::SolverUtils::EquationSystem::m_fields, m_PSduration, Nektar::LibUtilities::rad(), and tinysimd::sqrt().

Referenced by v_DoSolve(), and v_SetInitialConditions().

◆ GenerateSigmaPML()

void Nektar::MMFMaxwell::GenerateSigmaPML ( const NekDouble  PMLthickness,
const NekDouble  PMLstart,
const NekDouble  PMLmaxsigma,
Array< OneD, Array< OneD, NekDouble >> &  SigmaPML 
)
protected

Definition at line 2026 of file MMFMaxwell.cpp.

2030 {
2031  int nq = m_fields[0]->GetNpoints();
2032 
2033  // Construct sigmaX and sigmaY for UPML
2034  Array<OneD, NekDouble> x(nq);
2035  Array<OneD, NekDouble> y(nq);
2036  Array<OneD, NekDouble> z(nq);
2037 
2038  m_fields[0]->GetCoords(x, y, z);
2039 
2040  // Construction of SigmaPML
2041  SigmaPML = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
2042  for (int j = 0; j < m_shapedim; j++)
2043  {
2044  SigmaPML[j] = Array<OneD, NekDouble>(nq, 0.0);
2045  }
2046 
2047  // PML region indicator: [ 0 : curvedPML : RecPML]
2048  // RecPML = [0 : 0 : 1]
2049  // PMLRegion = [0 : 1 : 1], CurvedPML = PMLRegion - RecPML = [0: 1: 0]
2050  Array<OneD, NekDouble> RecPML(nq);
2051  Array<OneD, NekDouble> CurvedPML(nq);
2052  Array<OneD, NekDouble> PMLRegion(nq);
2053  m_fields[0]->GenerateElementVector(m_RecPML, 0.0, 1.0, RecPML);
2054  m_fields[0]->GenerateElementVector(m_PMLelement, 0.0, 1.0, PMLRegion);
2055  Vmath::Vsub(nq, PMLRegion, 1, RecPML, 1, CurvedPML, 1);
2056 
2057  switch (m_AddPML)
2058  {
2059  // RecPML only
2060  case 1:
2061  {
2062  // Rectangular PML
2063  NekDouble xRlayer, xLlayer, yRlayer, yLlayer;
2064 
2065  xRlayer = Vmath::Vmax(nq, x, 1) - PMLthickness;
2066  xLlayer = Vmath::Vmin(nq, x, 1) + PMLthickness;
2067  yRlayer = Vmath::Vmax(nq, y, 1) - PMLthickness;
2068  yLlayer = Vmath::Vmin(nq, y, 1) + PMLthickness;
2069 
2070  NekDouble xd, yd;
2071  for (int i = 0; i < nq; i++)
2072  {
2073  // SimgaPML along e^1
2074  if (x[i] >= xRlayer)
2075  {
2076  xd = (x[i] - xRlayer) / PMLthickness;
2077  }
2078 
2079  else if (x[i] <= xLlayer)
2080  {
2081  xd = (xLlayer - x[i]) / PMLthickness;
2082  }
2083 
2084  else
2085  {
2086  xd = 0.0;
2087  }
2088 
2089  SigmaPML[0][i] = RecPML[i] * PMLmaxsigma * (xd * xd * xd);
2090 
2091  // SigmaPML along e^2
2092  if (y[i] >= yRlayer)
2093  {
2094  yd = (y[i] - yRlayer) / PMLthickness;
2095  }
2096 
2097  else if (y[i] <= yLlayer)
2098  {
2099  yd = (yLlayer - y[i]) / PMLthickness;
2100  }
2101 
2102  else
2103  {
2104  yd = 0.0;
2105  }
2106 
2107  SigmaPML[1][i] = PMLRegion[i] * PMLmaxsigma * (yd * yd * yd);
2108  }
2109  }
2110  break;
2111 
2112  // CurvedPML only
2113  case 2:
2114  {
2115  // Curved PML
2116  NekDouble relrad, rad;
2117  for (int i = 0; i < nq; i++)
2118  {
2119  rad = sqrt(x[i] * x[i] / m_MMFfactors[0] / m_MMFfactors[0] +
2120  y[i] * y[i] / m_MMFfactors[1] / m_MMFfactors[1]);
2121 
2122  if (rad >= PMLstart)
2123  {
2124  relrad = (rad - PMLstart) / PMLthickness;
2125  SigmaPML[1][i] =
2126  PMLRegion[i] * PMLmaxsigma * pow(relrad, m_PMLorder);
2127  }
2128  }
2129  }
2130  break;
2131 
2132  // Slanted PML
2133  case 3:
2134  {
2135  NekDouble relrad, radon, radtw, radth, radfo;
2136  for (int i = 0; i < nq; i++)
2137  {
2138  radon = -1.0 * x[i] + y[i] - 7;
2139  radtw = x[i] + y[i] - 7;
2140  radth = -x[i] - y[i] - 7;
2141  radfo = x[i] - y[i] - 7;
2142 
2143  if (radon >= 0.0)
2144  {
2145  relrad = radon / PMLthickness;
2146  SigmaPML[1][i] =
2147  PMLRegion[i] * PMLmaxsigma * pow(relrad, m_PMLorder);
2148  }
2149 
2150  if (radtw >= 0.0)
2151  {
2152  relrad = radtw / PMLthickness;
2153  SigmaPML[0][i] =
2154  PMLRegion[i] * PMLmaxsigma * pow(relrad, m_PMLorder);
2155  }
2156 
2157  if (radth >= 0.0)
2158  {
2159  relrad = radth / PMLthickness;
2160  SigmaPML[0][i] =
2161  PMLRegion[i] * PMLmaxsigma * pow(relrad, m_PMLorder);
2162  }
2163 
2164  if (radfo >= 0.0)
2165  {
2166  relrad = radfo / PMLthickness;
2167  SigmaPML[1][i] =
2168  PMLRegion[i] * PMLmaxsigma * pow(relrad, m_PMLorder);
2169  }
2170  }
2171  }
2172  break;
2173  }
2174 
2175  std::cout << "*** sigma1 = [ " << Vmath::Vmin(nq, &SigmaPML[0][0], 1)
2176  << " , " << Vmath::Vmax(nq, &SigmaPML[0][0], 1)
2177  << " ] , sigma2 = [ " << Vmath::Vmin(nq, &SigmaPML[1][0], 1)
2178  << " , " << Vmath::Vmax(nq, &SigmaPML[1][0], 1) << " ] "
2179  << std::endl;
2180 }
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 m_AddPML, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_MMFfactors, m_PMLelement, m_PMLorder, m_RecPML, Nektar::SolverUtils::MMFSystem::m_shapedim, Nektar::LibUtilities::rad(), tinysimd::sqrt(), Vmath::Vmax(), Vmath::Vmin(), and Vmath::Vsub().

Referenced by v_InitObject().

◆ print_MMF()

void Nektar::MMFMaxwell::print_MMF ( Array< OneD, Array< OneD, NekDouble >> &  inarray)
protected

Definition at line 3135 of file MMFMaxwell.cpp.

3136 {
3137  int Ntot = inarray.size();
3138 
3139  NekDouble reval = 0.0;
3140  for (int i = 0; i < Ntot; ++i)
3141  {
3142  std::cout << "[" << i << "] = " << inarray[2][i] << std::endl;
3143  // reval = reval + inarray[i]*inarray[i];
3144  }
3145  reval = sqrt(reval / Ntot);
3146 }

References tinysimd::sqrt().

◆ Printout_SurfaceCurrent()

void Nektar::MMFMaxwell::Printout_SurfaceCurrent ( Array< OneD, Array< OneD, NekDouble >> &  fields,
const int  time 
)
protected

Definition at line 1903 of file MMFMaxwell.cpp.

1905 {
1906  int nq = m_fields[0]->GetTotPoints();
1907  int nTraceNumPoints = GetTraceTotPoints();
1908 
1909  int totbdryexp =
1910  m_fields[0]->GetBndCondExpansions()[m_boundaryforSF]->GetExpSize();
1911  int npts = m_fields[0]
1912  ->GetBndCondExpansions()[m_boundaryforSF]
1913  ->GetExp(0)
1914  ->GetNumPoints(0);
1915  int totnpts = totbdryexp * npts;
1916 
1917  Array<OneD, NekDouble> Jphi(totnpts);
1918  Array<OneD, NekDouble> Jrad(totnpts);
1919  Array<OneD, NekDouble> Jcurrent(totnpts);
1920 
1921  Array<OneD, NekDouble> phiFwd(nTraceNumPoints);
1922  Array<OneD, NekDouble> radFwd(nTraceNumPoints);
1923 
1924  // Compute phiFwd = acos(-x/r) along the trace
1925  Array<OneD, NekDouble> x(nq);
1926  Array<OneD, NekDouble> y(nq);
1927  Array<OneD, NekDouble> z(nq);
1928 
1929  m_fields[0]->GetCoords(x, y, z);
1930 
1931  Array<OneD, NekDouble> xFwd(nTraceNumPoints);
1932  Array<OneD, NekDouble> yFwd(nTraceNumPoints);
1933 
1934  m_fields[0]->ExtractTracePhys(x, xFwd);
1935  m_fields[0]->ExtractTracePhys(y, yFwd);
1936 
1937  for (int i = 0; i < nTraceNumPoints; ++i)
1938  {
1939  radFwd[i] = sqrt(xFwd[i] * xFwd[i] / m_MMFfactors[0] / m_MMFfactors[0] +
1940  yFwd[i] * yFwd[i] / m_MMFfactors[1] / m_MMFfactors[1]);
1941  phiFwd[i] = atan2(yFwd[i] / radFwd[i], -1.0 * xFwd[i] / radFwd[i]);
1942  }
1943 
1944  Array<OneD, NekDouble> ntimesHFwd(nTraceNumPoints);
1945  ntimesHFwd = ComputeSurfaceCurrent(time, fields);
1946 
1947  // The surface for current should be the first boundary
1948  int id2, cnt = 0;
1949  for (int e = 0; e < totbdryexp; ++e)
1950  {
1951  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
1952  m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1953 
1954  Vmath::Vcopy(npts, &phiFwd[id2], 1, &Jphi[e * npts], 1);
1955  Vmath::Vcopy(npts, &radFwd[id2], 1, &Jrad[e * npts], 1);
1956  Vmath::Vcopy(npts, &ntimesHFwd[id2], 1, &Jcurrent[e * npts], 1);
1957  }
1958 
1959  // Vmath::Vmul(totnpts, tmpr, 1, tmpr, 1, Jcurrent, 1);
1960  // Vmath::Vvtvp(totnpts, tmpi, 1, tmpi, 1, Jcurrent, 1, Jcurrent, 1);
1961  // Vmath::Vsqrt(totnpts, Jcurrent, 1, Jcurrent, 1);
1962 
1963  std::cout << "========================================================"
1964  << std::endl;
1965 
1966  std::cout << "phi = " << std::endl;
1967  for (int i = 0; i < totnpts; ++i)
1968  {
1969  std::cout << Jphi[i] << ", ";
1970  }
1971  std::cout << std::endl << std::endl;
1972 
1973  std::cout << "J = " << std::endl;
1974  for (int i = 0; i < totnpts; ++i)
1975  {
1976  std::cout << Jcurrent[i] << ", ";
1977  }
1978  std::cout << std::endl << std::endl;
1979 }
Array< OneD, NekDouble > ComputeSurfaceCurrent(const int time, const Array< OneD, const Array< OneD, NekDouble >> &fields)

References ComputeSurfaceCurrent(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), m_boundaryforSF, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_MMFfactors, tinysimd::sqrt(), and Vmath::Vcopy().

Referenced by v_DoSolve().

◆ TestMaxwell1D()

Array< OneD, NekDouble > Nektar::MMFMaxwell::TestMaxwell1D ( const NekDouble  time,
unsigned int  field 
)
protected

Definition at line 1466 of file MMFMaxwell.cpp.

1468 {
1469  int nq = m_fields[0]->GetNpoints();
1470 
1471  Array<OneD, NekDouble> x0(nq);
1472  Array<OneD, NekDouble> x1(nq);
1473  Array<OneD, NekDouble> x2(nq);
1474 
1475  m_fields[0]->GetCoords(x0, x1, x2);
1476 
1477  Array<OneD, NekDouble> E(nq);
1478  Array<OneD, NekDouble> H(nq);
1479 
1480  // Derive the frequency \omega
1481  NekDouble omega;
1482  NekDouble Tol = 0.000000001;
1483  if (fabs(m_n1 - m_n2) < Tol)
1484  {
1485  omega = m_pi / m_n1;
1486  }
1487 
1488  else
1489  {
1490  omega = 2.0 * m_pi / m_n2;
1491 
1492  NekDouble newomega, F, Fprime;
1493  for (int i = 0; i < 10000; ++i)
1494  {
1495  F = m_n1 * tan(m_n2 * omega) + m_n2 * tan(m_n1 * omega);
1496  Fprime = m_n1 * m_n2 *
1497  (1.0 / cos(m_n2 * omega) / cos(m_n2 * omega) +
1498  1.0 / cos(m_n1 * omega) / cos(m_n1 * omega));
1499 
1500  newomega = omega - F / Fprime;
1501 
1502  if (fabs(newomega - omega) > Tol)
1503  {
1504  omega = newomega;
1505  }
1506 
1507  else
1508  {
1509  break;
1510  }
1511  }
1512  }
1513 
1514  // Generate A^k and B^k
1515  std::complex<double> im = sqrt(std::complex<double>(-1));
1516  std::complex<double> A1, A2, B1, B2;
1517  std::complex<double> Ak, Bk, nk;
1518  std::complex<double> Ec, Hc;
1519 
1520  A1 = m_n2 * cos(m_n2 * omega) / (m_n1 * cos(m_n1 * omega));
1521  A2 = exp(-1.0 * im * omega * (m_n1 + m_n2));
1522  B1 = A1 * exp(-2.0 * im * m_n1 * omega);
1523  B2 = A2 * exp(2.0 * im * m_n2 * omega);
1524 
1525  for (int i = 0; i < nq; ++i)
1526  {
1527  if (x0[i] > 0)
1528  {
1529  Ak = A2;
1530  Bk = B2;
1531  nk = m_n2;
1532  }
1533 
1534  else
1535  {
1536  Ak = A1;
1537  Bk = B1;
1538  nk = m_n1;
1539  }
1540 
1541  Ec = (Ak * exp(im * nk * omega * x0[i]) -
1542  Bk * exp(-im * nk * omega * x0[i])) *
1543  exp(im * omega * time);
1544  Hc = nk *
1545  (Ak * exp(im * nk * omega * x0[i]) +
1546  Bk * exp(-im * nk * omega * x0[i])) *
1547  exp(im * omega * time);
1548 
1549  E[i] = Ec.real();
1550  H[i] = Hc.real();
1551  }
1552 
1553  Array<OneD, NekDouble> outfield;
1554  switch (field)
1555  {
1556  case (0):
1557  {
1558  outfield = E;
1559  }
1560  break;
1561 
1562  case (1):
1563  {
1564  outfield = H;
1565  }
1566  break;
1567  }
1568 
1569  return outfield;
1570 }

References Nektar::SolverUtils::EquationSystem::m_fields, m_n1, m_n2, Nektar::SolverUtils::MMFSystem::m_pi, and tinysimd::sqrt().

Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().

◆ TestMaxwell2DPEC()

Array< OneD, NekDouble > Nektar::MMFMaxwell::TestMaxwell2DPEC ( const NekDouble  time,
unsigned int  field,
const SolverUtils::PolType  Polarization 
)
protected

Definition at line 1572 of file MMFMaxwell.cpp.

1575 {
1576  int nq = m_fields[0]->GetNpoints();
1577 
1578  Array<OneD, NekDouble> x0(nq);
1579  Array<OneD, NekDouble> x1(nq);
1580  Array<OneD, NekDouble> x2(nq);
1581 
1582  m_fields[0]->GetCoords(x0, x1, x2);
1583 
1584  NekDouble freqm = 1.0, freqn = 1.0;
1585  NekDouble omega = m_pi * sqrt(freqm * freqm + freqn * freqn);
1586  NekDouble mpi = freqm * m_pi;
1587  NekDouble npi = freqn * m_pi;
1588 
1589  Array<OneD, NekDouble> F1(nq);
1590  Array<OneD, NekDouble> F2(nq);
1591  Array<OneD, NekDouble> Fz(nq);
1592  Array<OneD, NekDouble> dF1dt(nq);
1593  Array<OneD, NekDouble> dF2dt(nq);
1594  Array<OneD, NekDouble> dFzdt(nq);
1595  NekDouble Fx, Fy, dFxdt, dFydt;
1596 
1597  for (int i = 0; i < nq; ++i)
1598  {
1599  switch (Polarization)
1600  {
1602  {
1603  Fx = -1.0 * (npi / omega) * sin(mpi * x0[i]) *
1604  cos(npi * x1[i]) * sin(omega * time);
1605  Fy = (mpi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1606  sin(omega * time);
1607 
1608  F1[i] =
1609  Fx * m_movingframes[0][i] + Fy * m_movingframes[0][i + nq];
1610  F2[i] =
1611  Fx * m_movingframes[1][i] + Fy * m_movingframes[1][i + nq];
1612  Fz[i] = sin(mpi * x0[i]) * sin(npi * x1[i]) * cos(omega * time);
1613 
1614  dFxdt = (npi)*sin(mpi * x0[i]) * cos(npi * x1[i]) *
1615  cos(omega * time);
1616  dFydt = -1.0 * (mpi)*cos(mpi * x0[i]) * sin(npi * x1[i]) *
1617  cos(omega * time);
1618 
1619  dF1dt[i] = dFxdt * m_movingframes[0][i] +
1620  dFydt * m_movingframes[0][i + nq];
1621  dF2dt[i] = dFxdt * m_movingframes[1][i] +
1622  dFydt * m_movingframes[1][i + nq];
1623  dFzdt[i] = omega * sin(mpi * x0[i]) * sin(npi * x1[i]) *
1624  sin(omega * time);
1625  }
1626  break;
1627 
1629  {
1630  Fx = -1.0 * (npi / omega) * cos(mpi * x0[i]) *
1631  sin(npi * x1[i]) * sin(omega * time);
1632  Fy = (mpi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1633  sin(omega * time);
1634 
1635  F1[i] =
1636  Fx * m_movingframes[0][i] + Fy * m_movingframes[0][i + nq];
1637  F2[i] =
1638  Fx * m_movingframes[1][i] + Fy * m_movingframes[1][i + nq];
1639  Fz[i] = cos(mpi * x0[i]) * cos(npi * x1[i]) * cos(omega * time);
1640 
1641  dFxdt = (npi)*cos(mpi * x0[i]) * sin(npi * x1[i]) *
1642  cos(omega * time);
1643  dFydt = -1.0 * (mpi)*sin(mpi * x0[i]) * cos(npi * x1[i]) *
1644  cos(omega * time);
1645 
1646  dF1dt[i] = dFxdt * m_movingframes[0][i] +
1647  dFydt * m_movingframes[0][i + nq];
1648  dF2dt[i] = dFxdt * m_movingframes[1][i] +
1649  dFydt * m_movingframes[1][i + nq];
1650  dFzdt[i] = omega * cos(mpi * x0[i]) * cos(npi * x1[i]) *
1651  sin(omega * time);
1652  }
1653  break;
1654 
1655  default:
1656  break;
1657  }
1658  }
1659 
1660  Array<OneD, NekDouble> outfield;
1661  switch (field)
1662  {
1663  case (0):
1664  {
1665  outfield = F1;
1666  }
1667  break;
1668 
1669  case (1):
1670  {
1671  outfield = F2;
1672  }
1673  break;
1674 
1675  case (2):
1676  {
1677  outfield = Fz;
1678  }
1679  break;
1680 
1681  case (10):
1682  {
1683  outfield = dF1dt;
1684  }
1685  break;
1686 
1687  case (11):
1688  {
1689  outfield = dF2dt;
1690  }
1691  break;
1692 
1693  case (12):
1694  {
1695  outfield = dFzdt;
1696  }
1697  break;
1698  }
1699 
1700  return outfield;
1701 }

References Nektar::SolverUtils::eTransElectric, Nektar::SolverUtils::eTransMagnetic, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_pi, and tinysimd::sqrt().

Referenced by DoOdeRhs(), v_EvaluateExactSolution(), and v_SetInitialConditions().

◆ TestMaxwell2DPMC()

Array< OneD, NekDouble > Nektar::MMFMaxwell::TestMaxwell2DPMC ( const NekDouble  time,
unsigned int  field,
const SolverUtils::PolType  Polarization 
)
protected

Definition at line 1703 of file MMFMaxwell.cpp.

1706 {
1707  int nq = m_fields[0]->GetNpoints();
1708 
1709  Array<OneD, NekDouble> x0(nq);
1710  Array<OneD, NekDouble> x1(nq);
1711  Array<OneD, NekDouble> x2(nq);
1712 
1713  m_fields[0]->GetCoords(x0, x1, x2);
1714 
1715  NekDouble freqm = 1.0, freqn = 1.0;
1716  NekDouble omega = m_pi * sqrt(freqm * freqm + freqn * freqn);
1717  NekDouble mpi = freqm * m_pi;
1718  NekDouble npi = freqn * m_pi;
1719 
1720  Array<OneD, NekDouble> F1(nq);
1721  Array<OneD, NekDouble> F2(nq);
1722  Array<OneD, NekDouble> Fz(nq);
1723  NekDouble Fx, Fy;
1724 
1725  for (int i = 0; i < nq; ++i)
1726  {
1727  switch (Polarization)
1728  {
1730  {
1731  Fx = (npi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1732  sin(omega * time);
1733  Fy = -(mpi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1734  sin(omega * time);
1735 
1736  F1[i] =
1737  Fx * m_movingframes[0][i] + Fy * m_movingframes[0][i + nq];
1738  F2[i] =
1739  Fx * m_movingframes[1][i] + Fy * m_movingframes[1][i + nq];
1740  Fz[i] = cos(mpi * x0[i]) * cos(npi * x1[i]) * cos(omega * time);
1741  }
1742  break;
1743 
1745  {
1746  Fx = (npi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1747  sin(omega * time);
1748  Fy = -(mpi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1749  sin(omega * time);
1750 
1751  F1[i] =
1752  Fx * m_movingframes[0][i] + Fy * m_movingframes[0][i + nq];
1753  F2[i] =
1754  Fx * m_movingframes[1][i] + Fy * m_movingframes[1][i + nq];
1755  Fz[i] = sin(mpi * x0[i]) * sin(npi * x1[i]) * cos(omega * time);
1756  }
1757  break;
1758 
1759  default:
1760  break;
1761  }
1762  }
1763 
1764  Array<OneD, NekDouble> outfield;
1765  switch (field)
1766  {
1767  case (0):
1768  {
1769  outfield = F1;
1770  }
1771  break;
1772 
1773  case (1):
1774  {
1775  outfield = F2;
1776  }
1777  break;
1778 
1779  case (2):
1780  {
1781  outfield = Fz;
1782  }
1783  break;
1784  }
1785 
1786  return outfield;
1787 }

References Nektar::SolverUtils::eTransElectric, Nektar::SolverUtils::eTransMagnetic, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_pi, and tinysimd::sqrt().

Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().

◆ TestMaxwellSphere()

Array< OneD, NekDouble > Nektar::MMFMaxwell::TestMaxwellSphere ( const NekDouble  time,
const NekDouble  omega,
unsigned int  field 
)
protected

Definition at line 1789 of file MMFMaxwell.cpp.

1792 {
1793  int nq = m_fields[0]->GetTotPoints();
1794 
1795  Array<OneD, NekDouble> outfield(nq);
1796 
1797  Array<OneD, NekDouble> x(nq);
1798  Array<OneD, NekDouble> y(nq);
1799  Array<OneD, NekDouble> z(nq);
1800 
1801  m_fields[0]->GetCoords(x, y, z);
1802 
1803  Array<OneD, NekDouble> H1(nq);
1804  Array<OneD, NekDouble> H2(nq);
1805  Array<OneD, NekDouble> E3(nq);
1806  Array<OneD, NekDouble> F1(nq);
1807  Array<OneD, NekDouble> F2(nq);
1808 
1809  Array<OneD, NekDouble> curlv(nq);
1810  Array<OneD, Array<OneD, NekDouble>> velvec(m_spacedim);
1811  Array<OneD, Array<OneD, NekDouble>> Fvec(m_spacedim);
1812  for (int i = 0; i < m_spacedim; ++i)
1813  {
1814  velvec[i] = Array<OneD, NekDouble>(nq, 0.0);
1815  Fvec[i] = Array<OneD, NekDouble>(nq, 0.0);
1816  }
1817 
1818  NekDouble xj, yj, zj, sin_varphi, cos_varphi, sin_theta, cos_theta;
1819  NekDouble vth, vphi, Fth, Fphi;
1820  for (int i = 0; i < nq; i++)
1821  {
1822  xj = x[i];
1823  yj = y[i];
1824  zj = z[i];
1825 
1826  CartesianToSpherical(xj, yj, zj, sin_varphi, cos_varphi, sin_theta,
1827  cos_theta);
1828 
1829  vth = -4.0 * sin_varphi * cos_varphi * cos_theta * cos_theta *
1830  cos_theta * sin_theta;
1831  vphi =
1832  -1.0 * sin_varphi * sin_varphi * cos_theta * cos_theta * cos_theta;
1833  velvec[0][i] = -vth * sin_theta * cos_varphi - vphi * sin_varphi;
1834  velvec[1][i] = -vth * sin_theta * sin_varphi + vphi * cos_varphi;
1835  velvec[2][i] = vth * cos_theta;
1836 
1837  E3[i] = (-4.0 * cos_theta * cos_theta * sin_theta * cos_varphi *
1838  cos_varphi) *
1839  (1.0 / omega * sin(omega * time));
1840 
1841  Fth = -omega * vth -
1842  (8.0 / omega) * cos_theta * sin_theta * cos_varphi * sin_varphi;
1843  Fphi = -omega * vphi +
1844  (4.0 / omega) * cos_varphi * cos_varphi * cos_theta *
1845  (2.0 * sin_theta * sin_theta - cos_theta * cos_theta);
1846  Fvec[0][i] = -Fth * sin_theta * cos_varphi - Fphi * sin_varphi;
1847  Fvec[1][i] = -Fth * sin_theta * sin_varphi + Fphi * cos_varphi;
1848  Fvec[2][i] = Fth * cos_theta;
1849  }
1850 
1851  H1 = CartesianToMovingframes(velvec, 0);
1852  H2 = CartesianToMovingframes(velvec, 1);
1853 
1854  Vmath::Smul(nq, cos(omega * time), H1, 1, H1, 1);
1855  Vmath::Smul(nq, cos(omega * time), H2, 1, H2, 1);
1856 
1857  F1 = CartesianToMovingframes(Fvec, 0);
1858  F2 = CartesianToMovingframes(Fvec, 1);
1859 
1860  Vmath::Smul(nq, sin(omega * time), F1, 1, F1, 1);
1861  Vmath::Smul(nq, sin(omega * time), F2, 1, F2, 1);
1862 
1863  switch (field)
1864  {
1865  // return H1
1866  case 0:
1867  {
1868  outfield = H1;
1869  }
1870  break;
1871 
1872  case 1:
1873  {
1874  outfield = H2;
1875  }
1876  break;
1877 
1878  case 2:
1879  {
1880  outfield = E3;
1881  }
1882  break;
1883 
1884  case 3:
1885  {
1886  outfield = F1;
1887  }
1888  break;
1889 
1890  case 4:
1891  {
1892  outfield = F2;
1893  }
1894  break;
1895 
1896  default:
1897  break;
1898  }
1899 
1900  return outfield;
1901 }
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
Definition: MMFSystem.cpp:774

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

Referenced by DoOdeRhs(), v_EvaluateExactSolution(), and v_SetInitialConditions().

◆ v_DoSolve()

void Nektar::MMFMaxwell::v_DoSolve ( void  )
virtual

Solves an unsteady problem.

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

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 449 of file MMFMaxwell.cpp.

450 {
451  ASSERTL0(m_intScheme != 0, "No time integration scheme.");
452 
453  int i, nchk = 1;
454  int nq = GetTotPoints();
455  int nvariables = 0;
456  int nfields = m_fields.size();
457 
458  if (m_intVariables.empty())
459  {
460  for (i = 0; i < nfields; ++i)
461  {
462  m_intVariables.push_back(i);
463  }
464  nvariables = nfields;
465  }
466  else
467  {
468  nvariables = m_intVariables.size();
469  }
470 
471  // Set up wrapper to fields data storage.
472  Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
473  Array<OneD, Array<OneD, NekDouble>> tmp(nvariables);
474 
475  // Order storage to list time-integrated fields first.
476  for (i = 0; i < nvariables; ++i)
477  {
478  fields[i] = m_fields[m_intVariables[i]]->GetPhys();
479  m_fields[m_intVariables[i]]->SetPhysState(false);
480  }
481 
482  // Initialise time integration scheme
483  m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
484 
485  // Check uniqueness of checkpoint output
486  ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
487  (m_checktime > 0.0 && m_checksteps == 0) ||
488  (m_checktime == 0.0 && m_checksteps > 0),
489  "Only one of IO_CheckTime and IO_CheckSteps "
490  "should be set!");
491 
492  int Ntot = m_steps / m_checksteps + 1;
493 
494  Array<OneD, NekDouble> TimeSeries(Ntot);
495  Array<OneD, NekDouble> Energy(Ntot);
496 
497  LibUtilities::Timer timer;
498  bool doCheckTime = false;
499  int step = 0;
500  NekDouble intTime = 0.0;
501  NekDouble cpuTime = 0.0;
502  NekDouble elapsed = 0.0;
503 
504  int cntap = 0;
505  Array<OneD, NekDouble> Ezantipod;
506  int indxantipod = 0;
507 
508  switch (m_SourceType)
509  {
510  case ePointSource:
511  {
512  Ezantipod = Array<OneD, NekDouble>(m_steps / m_checksteps);
513 
514  Array<OneD, NekDouble> x(nq);
515  Array<OneD, NekDouble> y(nq);
516  Array<OneD, NekDouble> z(nq);
517 
518  m_fields[0]->GetCoords(x, y, z);
519 
520  NekDouble Tol = 0.000001;
521  NekDouble rad;
522  for (i = 0; i < nq; ++i)
523  {
524  rad = sqrt((x[i] + m_Psx) * (x[i] + m_Psx) +
525  (y[i] + m_Psy) * (y[i] + m_Psy) +
526  (z[i] + m_Psz) * (z[i] + m_Psz));
527  std::cout << "rad" << rad << std::endl;
528  if (rad < Tol)
529  {
530  indxantipod = i;
531  break;
532  }
533  }
534  }
535  break;
536 
537  case ePlanarSource:
538  {
539  m_SourceVector = Array<OneD, NekDouble>(nq, 0.0);
540 
541  Array<OneD, NekDouble> x(nq);
542  Array<OneD, NekDouble> y(nq);
543  Array<OneD, NekDouble> z(nq);
544 
545  m_fields[0]->GetCoords(x, y, z);
546 
547  NekDouble Tol = 0.000001;
548  NekDouble rad;
549  for (i = 0; i < nq; ++i)
550  {
551  rad = sqrt((x[i] - m_Psx) * (x[i] - m_Psx));
552  if (rad < Tol)
553  {
554  m_SourceVector[i] = 1.0;
555  }
556  }
557 
558  std::cout << "*** Area of Planar Source = "
559  << m_fields[0]->Integral(m_SourceVector) << std::endl;
560  }
561  break;
562 
563  default:
564  break;
565  }
566 
567  int cntpml = 0;
568  int P1indx = 0, P2indx = 0, P3indx = 0;
569  Array<OneD, NekDouble> P1;
570  Array<OneD, NekDouble> P2;
571  Array<OneD, NekDouble> P3;
572  if (m_TestPML)
573  {
574  P1 = Array<OneD, NekDouble>(m_steps / m_checksteps);
575  P2 = Array<OneD, NekDouble>(m_steps / m_checksteps);
576  P3 = Array<OneD, NekDouble>(m_steps / m_checksteps);
577 
578  Array<OneD, NekDouble> x(nq);
579  Array<OneD, NekDouble> y(nq);
580  Array<OneD, NekDouble> z(nq);
581 
582  m_fields[0]->GetCoords(x, y, z);
583 
584  NekDouble Tol = 0.000001;
585  NekDouble rad;
586  for (int i = 0; i < nq; ++i)
587  {
588  rad = sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i]) * (y[i]));
589 
590  if (rad < Tol)
591  {
592  P1indx = i;
593  break;
594  }
595  }
596 
597  for (int i = 0; i < nq; ++i)
598  {
599  rad =
600  sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i] - 1.5) * (y[i] - 1.5));
601  if (rad < Tol)
602  {
603  P2indx = i;
604  break;
605  }
606  }
607 
608  for (int i = 0; i < nq; ++i)
609  {
610  rad =
611  sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i] - 3.0) * (y[i] - 3.0));
612  if (rad < Tol)
613  {
614  P3indx = i;
615  break;
616  }
617  }
618  }
619 
620  int indx;
621  while (step < m_steps || m_time < m_fintime - NekConstants::kNekZeroTol)
622  {
623 
624  timer.Start();
625  fields = m_intScheme->TimeIntegrate(step, m_timestep, m_ode);
626  timer.Stop();
627 
628  m_time += m_timestep;
629  elapsed = timer.TimePerTest(1);
630  intTime += elapsed;
631  cpuTime += elapsed;
632 
633  // Write out status information
634  if (m_session->GetComm()->GetRank() == 0 && !((step + 1) % m_infosteps))
635  {
636  std::cout << "Steps: " << std::setw(8) << std::left << step + 1
637  << " "
638  << "Time: " << std::setw(12) << std::left << m_time;
639 
640  std::stringstream ss;
641  ss << cpuTime / 60.0 << " min.";
642  std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
643  << std::endl;
644 
645  cpuTime = 0.0;
646  }
647 
648  switch (m_SourceType)
649  {
650  case ePointSource:
651  {
652  if (m_time <= m_PSduration)
653  {
654  Array<OneD, NekDouble> Impulse(nq);
655  Impulse = GaussianPulse(m_time, m_Psx, m_Psy, m_Psz,
657  Vmath::Vadd(nq, &Impulse[0], 1,
658  &fields[m_intVariables[2]][0], 1,
659  &fields[m_intVariables[2]][0], 1);
660  }
661  }
662  break;
663 
664  case ePlanarSource:
665  {
666  Array<OneD, NekDouble> Impulse(nq);
667  for (int i = 0; i < 3; ++i)
668  {
669  Impulse = GetIncidentField(i, m_time);
670  Vmath::Vmul(nq, m_SourceVector, 1, Impulse, 1, Impulse, 1);
671  Vmath::Vadd(nq, &Impulse[0], 1,
672  &fields[m_intVariables[i]][0], 1,
673  &fields[m_intVariables[i]][0], 1);
674  }
675  }
676  break;
677 
678  default:
679  break;
680  }
681 
682  // Transform data into coefficient space
683  for (i = 0; i < nvariables; ++i)
684  {
685  m_fields[m_intVariables[i]]->SetPhys(fields[i]);
686  m_fields[m_intVariables[i]]->FwdTransLocalElmt(
687  fields[i], m_fields[m_intVariables[i]]->UpdateCoeffs());
688  m_fields[m_intVariables[i]]->SetPhysState(false);
689  }
690  // for (i = 0; i < nq; ++i)
691  // std::cout << m_fields[0][0][i] <<std::endl;
692 
693  // Write out checkpoint files
694  if ((m_checksteps && step && !((step + 1) % m_checksteps)) ||
695  doCheckTime)
696  {
697  indx = (step + 1) / m_checksteps;
698  TimeSeries[indx] = m_time;
699 
701  {
702  Checkpoint_TotalFieldOutput(nchk, m_time, fields);
703  Checkpoint_TotPlotOutput(nchk, m_time, fields);
704  }
705  Checkpoint_PlotOutput(nchk, fields);
706  Checkpoint_EDFluxOutput(nchk, m_time, fields);
707  Checkpoint_EnergyOutput(nchk, m_time, fields);
708  Checkpoint_Output(nchk++);
709 
710  Energy[indx] = ComputeEnergyDensity(fields);
711 
712  std::cout << "|EHr|: F1 = " << RootMeanSquare(fields[0])
713  << ", F2 = " << RootMeanSquare(fields[1])
714  << ", F3 = " << RootMeanSquare(fields[2])
715  << ", Energy = " << Energy[indx] << std::endl;
716  if (nfields > 3)
717  {
718  std::cout << "|DBr|: D1 = " << RootMeanSquare(fields[3])
719  << ", D2 = " << RootMeanSquare(fields[4])
720  << ", D3 = " << RootMeanSquare(fields[5])
721  << std::endl;
722 
723  int nTraceNumPoints = GetTraceNpoints();
724  int totbdryexp =
725  m_fields[0]->GetBndCondExpansions()[0]->GetExpSize();
726  int npts = m_fields[0]
727  ->GetBndCondExpansions()[0]
728  ->GetExp(0)
729  ->GetNumPoints(0);
730 
731  Array<OneD, NekDouble> x0(nq);
732  Array<OneD, NekDouble> x1(nq);
733  Array<OneD, NekDouble> x2(nq);
734 
735  m_fields[0]->GetCoords(x0, x1, x2);
736 
737  Array<OneD, NekDouble> E1Fwd(nTraceNumPoints);
738  Array<OneD, NekDouble> E2Fwd(nTraceNumPoints);
739  Array<OneD, NekDouble> H3Fwd(nTraceNumPoints);
740 
741  m_fields[0]->ExtractTracePhys(fields[0], E1Fwd);
742  m_fields[0]->ExtractTracePhys(fields[1], E2Fwd);
743  m_fields[0]->ExtractTracePhys(fields[2], H3Fwd);
744 
745  int id2, cnt = 0;
746  NekDouble E1atPECloc, E1atPEC = 0.0;
747  NekDouble E2atPECloc, E2atPEC = 0.0;
748  NekDouble H3atPECloc, H3atPEC = 0.0;
749 
750  Array<OneD, NekDouble> E1Fwdloc(npts);
751  Array<OneD, NekDouble> E2Fwdloc(npts);
752  Array<OneD, NekDouble> H3Fwdloc(npts);
753 
754  for (int e = 0; e < totbdryexp; ++e)
755  {
756  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
757  m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
758  cnt + e));
759 
760  Vmath::Vcopy(npts, &E1Fwd[id2], 1, &E1Fwdloc[0], 1);
761  Vmath::Vcopy(npts, &E2Fwd[id2], 1, &E2Fwdloc[0], 1);
762  Vmath::Vcopy(npts, &H3Fwd[id2], 1, &H3Fwdloc[0], 1);
763 
764  E1atPECloc = Vmath::Vamax(npts, E1Fwdloc, 1);
765  E2atPECloc = Vmath::Vamax(npts, E2Fwdloc, 1);
766  H3atPECloc = Vmath::Vamax(npts, H3Fwdloc, 1);
767 
768  if (E1atPEC < E1atPECloc)
769  {
770  E1atPEC = E1atPECloc;
771  }
772 
773  if (E2atPEC < E2atPECloc)
774  {
775  E2atPEC = E2atPECloc;
776  }
777 
778  if (H3atPEC < H3atPECloc)
779  {
780  H3atPEC = H3atPECloc;
781  }
782  }
783 
784  std::cout << "At PEC, Max. E1 = " << E1atPEC
785  << ", E2 = " << E2atPEC << ", H3 = " << H3atPEC
786  << std::endl;
787  }
788 
789  if (m_SourceType == ePointSource)
790  {
791  Ezantipod[cntap++] = fields[2][indxantipod];
792  }
793 
794  if (m_TestPML)
795  {
796  P1[cntpml] = fields[2][P1indx];
797  P2[cntpml] = fields[2][P2indx];
798  P3[cntpml] = fields[2][P3indx];
799  cntpml++;
800  }
801  doCheckTime = false;
802  }
803 
804  // Step advance
805  ++step;
806  }
807 
808  // Print out summary statistics
809  if (m_session->GetComm()->GetRank() == 0)
810  {
811  std::cout << "Time-integration : " << intTime << "s" << std::endl;
812 
813  std::cout << "TimeSeries = " << std::endl;
814  for (int i = 0; i < m_steps / m_checksteps; ++i)
815  {
816  std::cout << TimeSeries[i] << ", ";
817  }
818  std::cout << std::endl << std::endl;
819 
820  std::cout << "Energy Density = " << std::endl;
821  for (int i = 0; i < m_steps / m_checksteps; ++i)
822  {
823  std::cout << Energy[i] << ", ";
824  }
825  std::cout << std::endl << std::endl;
826 
828  {
830  }
831 
832  if (m_SourceType == ePointSource)
833  {
834  std::cout << "Ez at antipod = " << std::endl;
835  for (int i = 0; i < m_steps / m_checksteps; ++i)
836  {
837  std::cout << Ezantipod[i] << ", ";
838  }
839  std::cout << std::endl << std::endl;
840  }
841 
842  if (m_TestPML)
843  {
844  std::cout << "P1 = " << std::endl;
845  for (int i = 0; i < m_steps / m_checksteps; ++i)
846  {
847  std::cout << P1[i] << ", ";
848  }
849  std::cout << std::endl << std::endl;
850 
851  std::cout << "P2 = " << std::endl;
852  for (int i = 0; i < m_steps / m_checksteps; ++i)
853  {
854  std::cout << P2[i] << ", ";
855  }
856  std::cout << std::endl << std::endl;
857 
858  std::cout << "P3 = " << std::endl;
859  for (int i = 0; i < m_steps / m_checksteps; ++i)
860  {
861  std::cout << P3[i] << ", ";
862  }
863  std::cout << std::endl << std::endl;
864  }
865  }
866 
867  for (i = 0; i < nvariables; ++i)
868  {
869  m_fields[m_intVariables[i]]->SetPhys(fields[i]);
870  m_fields[m_intVariables[i]]->SetPhysState(true);
871  }
872 
873  for (i = 0; i < nvariables; ++i)
874  {
875  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
876  m_fields[i]->UpdateCoeffs());
877  }
878 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
@ ePointSource
Definition: MMFMaxwell.h:60
@ ePlanarSource
Definition: MMFMaxwell.h:61
NekDouble m_Gaussianradius
Definition: MMFMaxwell.h:123
void Checkpoint_TotalFieldOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
void Checkpoint_EnergyOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
void Checkpoint_TotPlotOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
void Checkpoint_PlotOutput(const int n, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
Array< OneD, NekDouble > m_SourceVector
Definition: MMFMaxwell.h:121
void Printout_SurfaceCurrent(Array< OneD, Array< OneD, NekDouble >> &fields, const int time)
SourceType m_SourceType
Definition: MMFMaxwell.h:79
int m_PrintoutSurfaceCurrent
Definition: MMFMaxwell.h:107
Array< OneD, NekDouble > GaussianPulse(const NekDouble time, const NekDouble Psx, const NekDouble Psy, const NekDouble Psz, const NekDouble Gaussianradius)
void Checkpoint_EDFluxOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
NekDouble ComputeEnergyDensity(Array< OneD, Array< OneD, NekDouble >> &fields)
SOLVER_UTILS_EXPORT int GetTraceNpoints()
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
NekDouble m_fintime
Finish time of the simulation.
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
NekDouble m_checktime
Time between checkpoints.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
int m_steps
Number of steps to take.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT NekDouble RootMeanSquare(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2344
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
int m_infosteps
Number of time steps between outputting status information.
static const NekDouble kNekZeroTol
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max.
Definition: Vmath.cpp:999

References ASSERTL0, Checkpoint_EDFluxOutput(), Checkpoint_EnergyOutput(), Nektar::SolverUtils::EquationSystem::Checkpoint_Output(), Checkpoint_PlotOutput(), Checkpoint_TotalFieldOutput(), Checkpoint_TotPlotOutput(), ComputeEnergyDensity(), ePlanarSource, ePointSource, Nektar::SolverUtils::eScatField2D, GaussianPulse(), Nektar::SolverUtils::MMFSystem::GetIncidentField(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::GetTraceNpoints(), Nektar::NekConstants::kNekZeroTol, Nektar::SolverUtils::EquationSystem::m_checksteps, Nektar::SolverUtils::EquationSystem::m_checktime, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_fintime, m_Gaussianradius, Nektar::SolverUtils::UnsteadySystem::m_infosteps, Nektar::SolverUtils::UnsteadySystem::m_intScheme, Nektar::SolverUtils::UnsteadySystem::m_intVariables, Nektar::SolverUtils::UnsteadySystem::m_ode, m_PrintoutSurfaceCurrent, m_PSduration, m_Psx, m_Psy, m_Psz, Nektar::SolverUtils::EquationSystem::m_session, m_SourceType, m_SourceVector, Nektar::SolverUtils::EquationSystem::m_steps, Nektar::SolverUtils::MMFSystem::m_TestMaxwellType, m_TestPML, Nektar::SolverUtils::EquationSystem::m_time, Nektar::SolverUtils::EquationSystem::m_timestep, Printout_SurfaceCurrent(), Nektar::LibUtilities::rad(), Nektar::SolverUtils::MMFSystem::RootMeanSquare(), tinysimd::sqrt(), Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), Nektar::LibUtilities::Timer::TimePerTest(), Vmath::Vadd(), Vmath::Vamax(), Vmath::Vcopy(), and Vmath::Vmul().

◆ v_EvaluateExactSolution()

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1424 of file MMFMaxwell.cpp.

1427 {
1428  int nq = m_fields[0]->GetNpoints();
1429  outfield = Array<OneD, NekDouble>(nq);
1430 
1431  switch (m_TestMaxwellType)
1432  {
1434  {
1435  outfield = TestMaxwell1D(time, field);
1436  }
1437  break;
1438 
1441  {
1442  outfield = TestMaxwell2DPEC(time, field, m_PolType);
1443  }
1444  break;
1445 
1447  {
1448  outfield = TestMaxwell2DPMC(time, field, m_PolType);
1449  }
1450  break;
1451 
1453  {
1454  outfield = TestMaxwellSphere(time, m_freq, field);
1455  }
1456  break;
1457 
1458  default:
1459  {
1460  outfield = Array<OneD, NekDouble>(nq, 0.0);
1461  }
1462  break;
1463  }
1464 }
Array< OneD, NekDouble > TestMaxwell2DPMC(const NekDouble time, unsigned int field, const SolverUtils::PolType Polarization)
Array< OneD, NekDouble > TestMaxwell1D(const NekDouble time, unsigned int field)

References Nektar::SolverUtils::eMaxwell1D, Nektar::SolverUtils::eMaxwellSphere, Nektar::SolverUtils::eTestMaxwell2DPEC, Nektar::SolverUtils::eTestMaxwell2DPECAVGFLUX, Nektar::SolverUtils::eTestMaxwell2DPMC, Nektar::SolverUtils::EquationSystem::m_fields, m_freq, Nektar::SolverUtils::MMFSystem::m_PolType, Nektar::SolverUtils::MMFSystem::m_TestMaxwellType, TestMaxwell1D(), TestMaxwell2DPEC(), TestMaxwell2DPMC(), and TestMaxwellSphere().

◆ v_GenerateSummary()

void Nektar::MMFMaxwell::v_GenerateSummary ( SolverUtils::SummaryList s)
protectedvirtual

Print Summary.

Reimplemented from Nektar::SolverUtils::MMFSystem.

Definition at line 3067 of file MMFMaxwell.cpp.

3068 {
3071  s, "TestMaxwellType",
3073  SolverUtils::AddSummaryItem(s, "PolType",
3075  SolverUtils::AddSummaryItem(s, "IncType",
3077 
3078  if (m_varepsilon[0] * m_varepsilon[1] * m_varepsilon[2] > 1.0)
3079  {
3080  SolverUtils::AddSummaryItem(s, "varepsilon1", m_varepsilon[0]);
3081  SolverUtils::AddSummaryItem(s, "varepsilon2", m_varepsilon[1]);
3082  SolverUtils::AddSummaryItem(s, "varepsilon3", m_varepsilon[2]);
3083  }
3084 
3085  if (m_mu[0] * m_mu[1] * m_mu[2] > 1.0)
3086  {
3087  SolverUtils::AddSummaryItem(s, "mu1", m_mu[0]);
3088  SolverUtils::AddSummaryItem(s, "mu2", m_mu[1]);
3089  SolverUtils::AddSummaryItem(s, "mu3", m_mu[2]);
3090  }
3091 
3092  if (m_boundaryforSF > 0)
3093  {
3094  SolverUtils::AddSummaryItem(s, "boundarySF", m_boundaryforSF);
3095  }
3096 
3097  if (m_ElemtGroup1 > 0)
3098  {
3099  SolverUtils::AddSummaryItem(s, "CloakNlayer", m_CloakNlayer);
3100  SolverUtils::AddSummaryItem(s, "ElemtGroup1", m_ElemtGroup1);
3101  }
3102 
3103  SolverUtils::AddSummaryItem(s, "AddRotation", m_AddRotation);
3104 
3105  if (m_AddPML > 0)
3106  {
3107  SolverUtils::AddSummaryItem(s, "AddPML", m_AddPML);
3108  SolverUtils::AddSummaryItem(s, "PMLelement", m_PMLelement);
3109  SolverUtils::AddSummaryItem(s, "RecPML", m_RecPML);
3110  SolverUtils::AddSummaryItem(s, "PMLorder", m_PMLorder);
3111  SolverUtils::AddSummaryItem(s, "PMLthickness", m_PMLthickness);
3112  SolverUtils::AddSummaryItem(s, "PMLstart", m_PMLstart);
3113  SolverUtils::AddSummaryItem(s, "PMLmaxsigma", m_PMLmaxsigma);
3114  }
3115 
3116  if (m_SourceType)
3117  {
3118  SolverUtils::AddSummaryItem(s, "SourceType",
3123  SolverUtils::AddSummaryItem(s, "PSduration", m_PSduration);
3124  SolverUtils::AddSummaryItem(s, "Gaussianradius", m_Gaussianradius);
3125  }
3126 
3127  if (m_CloakType)
3128  {
3130  SolverUtils::AddSummaryItem(s, "DispersiveCloak", m_DispersiveCloak);
3131  SolverUtils::AddSummaryItem(s, "CloakNlayer", m_CloakNlayer);
3132  SolverUtils::AddSummaryItem(s, "Cloakraddelta", m_Cloakraddelta);
3133  }
3134 }
const char *const CloakTypeMap[]
Definition: MMFMaxwell.h:51
const char *const SourceTypeMap[]
Definition: MMFMaxwell.h:65
NekDouble m_CloakNlayer
Definition: MMFMaxwell.h:115
CloakType m_CloakType
Definition: MMFMaxwell.h:78
NekDouble m_PMLmaxsigma
Definition: MMFMaxwell.h:139
NekDouble m_PMLstart
Definition: MMFMaxwell.h:139
NekDouble m_PMLthickness
Definition: MMFMaxwell.h:139
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
Definition: MMFSystem.cpp:2469
const char *const IncTypeMap[]
Definition: MMFSystem.h:139
const char *const TestMaxwellTypeMap[]
Definition: MMFSystem.h:111
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
const char *const PolTypeMap[]
Definition: MMFSystem.h:126

References Nektar::SolverUtils::AddSummaryItem(), CloakTypeMap, Nektar::SolverUtils::IncTypeMap, m_AddPML, m_AddRotation, m_boundaryforSF, m_CloakNlayer, m_Cloakraddelta, m_CloakType, m_DispersiveCloak, m_ElemtGroup1, m_Gaussianradius, Nektar::SolverUtils::MMFSystem::m_IncType, m_mu, m_PMLelement, m_PMLmaxsigma, m_PMLorder, m_PMLstart, m_PMLthickness, Nektar::SolverUtils::MMFSystem::m_PolType, m_PSduration, m_Psx, m_Psy, m_Psz, m_RecPML, m_SourceType, Nektar::SolverUtils::MMFSystem::m_TestMaxwellType, m_varepsilon, Nektar::SolverUtils::PolTypeMap, SourceTypeMap, Nektar::SolverUtils::TestMaxwellTypeMap, and Nektar::SolverUtils::MMFSystem::v_GenerateSummary().

◆ v_InitObject()

void Nektar::MMFMaxwell::v_InitObject ( bool  DeclareFields = true)
virtual

Initialise the object.

Initialisation object for the unsteady linear advection equation.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 65 of file MMFMaxwell.cpp.

66 {
67  // Call to the initialisation object
68  UnsteadySystem::v_InitObject(DeclareFields);
69 
70  int nq = m_fields[0]->GetNpoints();
71  int shapedim = m_fields[0]->GetShapeDimension();
72 
73  m_session->LoadParameter("ElemtGroup0", m_ElemtGroup0, 0);
74  m_session->LoadParameter("ElemtGroup1", m_ElemtGroup1, 0);
75  m_session->LoadParameter("boundaryforSF", m_boundaryforSF, 0);
76  m_session->LoadParameter("PrintoutSurfaceCurrent", m_PrintoutSurfaceCurrent,
77  0);
78 
79  m_session->LoadParameter("AddRotation", m_AddRotation, 0);
80 
81  m_session->LoadParameter("NoInc", m_NoInc, 0);
82 
83  // PML parameters
84  m_session->LoadParameter("TestPML", m_TestPML, 0);
85  m_session->LoadParameter("PMLelement", m_PMLelement, 0);
86  m_session->LoadParameter("RecPML", m_RecPML, 0);
87 
88  m_session->LoadParameter("AddPML", m_AddPML, 0);
89  if (m_AddPML == 1)
90  {
92  }
93  m_session->LoadParameter("PMLorder", m_PMLorder, 3);
94 
95  m_session->LoadParameter("PMLthickness", m_PMLthickness, 0.0);
96  m_session->LoadParameter("PMLstart", m_PMLstart, 0.0);
97  m_session->LoadParameter("PMLmaxsigma", m_PMLmaxsigma, 100.0);
98 
99  // Point Source parmaters
100  m_session->LoadParameter("Psx", m_Psx, 0.0);
101  m_session->LoadParameter("Psy", m_Psy, 0.0);
102  m_session->LoadParameter("Psz", m_Psz, 0.0);
103  m_session->LoadParameter("PSduration", m_PSduration, 1.0);
104  m_session->LoadParameter("Gaussianradius", m_Gaussianradius, 1.0);
105 
106  // Cloaking parameter
107  m_session->LoadParameter("CloakNlayer", m_CloakNlayer, 5);
108  m_session->LoadParameter("Cloakraddelta", m_Cloakraddelta, 0.0);
109 
110  m_varepsilon = Array<OneD, NekDouble>(m_spacedim);
111  m_session->LoadParameter("varepsilon1", m_varepsilon[0], 1.0);
112  m_session->LoadParameter("varepsilon2", m_varepsilon[1], 1.0);
113  m_session->LoadParameter("varepsilon3", m_varepsilon[2], 1.0);
114  m_n1 = sqrt(m_varepsilon[0]);
115  m_n2 = sqrt(m_varepsilon[1]);
116  m_n3 = sqrt(m_varepsilon[2]);
117 
118  m_mu = Array<OneD, NekDouble>(m_spacedim);
119  m_session->LoadParameter("mu1", m_mu[0], 1.0);
120  m_session->LoadParameter("mu2", m_mu[1], 1.0);
121  m_session->LoadParameter("mu3", m_mu[2], 1.0);
122 
123  Array<OneD, Array<OneD, NekDouble>> Anisotropy(shapedim);
124  for (int j = 0; j < shapedim; ++j)
125  {
126  Anisotropy[j] = Array<OneD, NekDouble>(nq, 1.0);
127  }
128 
129  // Add Rectangular PML
130  MMFSystem::MMFInitObject(Anisotropy, m_RecPML);
131 
132  // Compute the cross producted MF
134 
135  m_session->LoadParameter("Frequency", m_freq, 1.0);
136 
137  // Define TestMaxwellType
138  if (m_session->DefinesSolverInfo("TESTMAXWELLTYPE"))
139  {
140  std::string TestMaxwellTypeStr =
141  m_session->GetSolverInfo("TESTMAXWELLTYPE");
142  for (int i = 0; i < (int)SolverUtils::SIZE_TestMaxwellType; ++i)
143  {
144  if (boost::iequals(SolverUtils::TestMaxwellTypeMap[i],
145  TestMaxwellTypeStr))
146  {
148  break;
149  }
150  }
151  }
152 
153  else
154  {
156  }
157 
158  // Define Polarization
159  if (m_session->DefinesSolverInfo("POLTYPE"))
160  {
161  std::string PolTypeStr = m_session->GetSolverInfo("POLTYPE");
162  for (int i = 0; i < (int)SolverUtils::SIZE_PolType; ++i)
163  {
164  if (boost::iequals(SolverUtils::PolTypeMap[i], PolTypeStr))
165  {
167  break;
168  }
169  }
170  }
171  else
172  {
174  }
175 
176  // Define Incident wave Type
177  if (m_session->DefinesSolverInfo("INCTYPE"))
178  {
179  std::string IncTypeStr = m_session->GetSolverInfo("INCTYPE");
180  for (int i = 0; i < (int)SolverUtils::SIZE_IncType; ++i)
181  {
182  if (boost::iequals(SolverUtils::IncTypeMap[i], IncTypeStr))
183  {
185  break;
186  }
187  }
188  }
189  else
190  {
192  }
193 
194  // Define Cloak Type
195  if (m_session->DefinesSolverInfo("CLOAKTYPE"))
196  {
197  std::string CloakTypeStr = m_session->GetSolverInfo("CLOAKTYPE");
198  for (int i = 0; i < (int)SIZE_CloakType; ++i)
199  {
200  if (boost::iequals(CloakTypeMap[i], CloakTypeStr))
201  {
202  m_CloakType = (CloakType)i;
203  break;
204  }
205  }
206  }
207  else
208  {
209  m_CloakType = (CloakType)0;
210  }
211 
212  // Define Source Type
213  if (m_session->DefinesSolverInfo("SOURCETYPE"))
214  {
215  std::string SourceTypeStr = m_session->GetSolverInfo("SOURCETYPE");
216  for (int i = 0; i < (int)SIZE_SourceType; ++i)
217  {
218  if (boost::iequals(SourceTypeMap[i], SourceTypeStr))
219  {
221  break;
222  }
223  }
224  }
225  else
226  {
228  }
229 
230  // Compute n_timesMFFwd and m_times_timesMFFwd
231  ComputeNtimesMF();
232 
233  // Compute vaepsilon and mu vector (m_epsveci, m_muvec0);
234  m_epsvec = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
235  m_muvec = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
236  for (int k = 0; k < m_spacedim; ++k)
237  {
238  m_epsvec[k] = Array<OneD, NekDouble>(nq, 1.0);
239  m_muvec[k] = Array<OneD, NekDouble>(nq, 1.0);
240  }
241 
242  Array<OneD, NekDouble> radvec(nq);
243  m_DispersiveCloak = false;
244  switch (m_CloakType)
245  {
246  case eOpticalCloak:
247  {
248  radvec = ComputeRadCloak();
250  }
251  break;
252 
253  case eOpticalConstCloak:
254  {
255  radvec = ComputeRadCloak(m_CloakNlayer);
257 
258  std::cout << "*** rad = [ " << Vmath::Vmax(nq, radvec, 1) << " , "
259  << Vmath::Vmin(nq, radvec, 1) << " ) " << std::endl;
260  }
261  break;
262 
264  {
265  m_DispersiveCloak = true;
266  m_wp2Tol = 0.01;
267  radvec = ComputeRadCloak();
269 
270  std::cout << "*** rad = [ " << Vmath::Vmax(nq, radvec, 1) << " , "
271  << Vmath::Vmin(nq, radvec, 1) << " ) " << std::endl;
272  std::cout << "*** wp2 = [ " << Vmath::Vmax(nq, m_wp2, 1) << " , "
273  << Vmath::Vmin(nq, m_wp2, 1) << " ) " << std::endl;
274  }
275  break;
276 
277  case eMicroWaveCloak:
278  {
279  radvec = ComputeRadCloak();
281  }
282  break;
283 
284  default:
285  {
287  }
288  break;
289  }
290 
291  NekDouble eps1min, eps1max, eps2min, eps2max, eps3min, eps3max;
292  NekDouble mu1min, mu1max, mu2min, mu2max, mu3min, mu3max;
293 
294  eps1min = Vmath::Vmin(nq, m_epsvec[0], 1);
295  eps3min = Vmath::Vmin(nq, m_epsvec[2], 1);
296  eps1max = Vmath::Vmax(nq, m_epsvec[0], 1);
297  eps3max = Vmath::Vmax(nq, m_epsvec[2], 1);
298 
299  if (m_DispersiveCloak)
300  {
301  Array<OneD, NekDouble> realepsr(nq);
302  Vmath::Sadd(nq, -m_wp2Tol, m_wp2, 1, realepsr, 1);
303  Vmath::Smul(nq, 1.0 / (m_Incfreq * m_Incfreq), realepsr, 1, realepsr,
304  1);
305  Vmath::Neg(nq, realepsr, 1);
306  Vmath::Sadd(nq, 1.0, realepsr, 1, realepsr, 1);
307 
308  eps2min = Vmath::Vmin(nq, realepsr, 1);
309  eps2max = Vmath::Vmax(nq, realepsr, 1);
310  }
311 
312  else
313  {
314  eps2min = Vmath::Vmin(nq, m_epsvec[1], 1);
315  eps2max = Vmath::Vmax(nq, m_epsvec[1], 1);
316  }
317 
318  mu1min = Vmath::Vmin(nq, m_muvec[0], 1);
319  mu2min = Vmath::Vmin(nq, m_muvec[1], 1);
320  mu3min = Vmath::Vmin(nq, m_muvec[2], 1);
321  mu1max = Vmath::Vmax(nq, m_muvec[0], 1);
322  mu2max = Vmath::Vmax(nq, m_muvec[1], 1);
323  mu3max = Vmath::Vmax(nq, m_muvec[2], 1);
324 
325  std::cout << "muvec0 = " << RootMeanSquare(m_muvec[0])
326  << ", muvec1 = " << RootMeanSquare(m_muvec[1]) << std::endl;
327 
328  std::cout << "*** epsvec1 = [ " << eps1min << " , " << eps1max
329  << " ], epsvec2 = [ " << eps2min << " , " << eps2max
330  << " ], epsvec3 = [ " << eps3min << " , " << eps3max << " ] "
331  << std::endl;
332  std::cout << "*** muvec1 = [ " << mu1min << " , " << mu1max
333  << " ], muvec2 = [ " << mu2min << " , " << mu2max
334  << " ], muvec3 = [ " << mu3min << " , " << mu3max << " ] "
335  << std::endl;
336 
337  NekDouble dtFactor;
338  switch (m_PolType)
339  {
340  // eTransMagnetic
342  {
343  dtFactor = mu1min * eps3min;
344  if (mu1min > mu2min)
345  {
346  dtFactor = mu2min * eps3min;
347  }
348  }
349  break;
350 
352  {
353  dtFactor = eps1min * mu3min;
354  if (eps1min > eps2min)
355  {
356  dtFactor = eps2min * mu3min;
357  }
358  }
359  break;
360 
361  default:
362  {
363  dtFactor = 1.0;
364  }
365  break;
366  }
367  std::cout << "*** dt factor proportional to varepsilon * mu is " << dtFactor
368  << std::endl;
369 
370  // Compute m_Zim and m_Yim
372 
373  // Compute m_epsvecminus1 and m_muminus1
374  m_negepsvecminus1 = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
375  m_negmuvecminus1 = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
376  for (int k = 0; k < m_spacedim; ++k)
377  {
378  m_negepsvecminus1[k] = Array<OneD, NekDouble>(nq, 0.0);
379  m_negmuvecminus1[k] = Array<OneD, NekDouble>(nq, 0.0);
380 
381  if (!m_NoInc)
382  {
383  Vmath::Sadd(nq, -1.0, m_muvec[k], 1, m_negmuvecminus1[k], 1);
384  Vmath::Sadd(nq, -1.0, m_epsvec[k], 1, m_negepsvecminus1[k], 1);
385 
386  Vmath::Neg(nq, m_negmuvecminus1[k], 1);
387  Vmath::Neg(nq, m_negepsvecminus1[k], 1);
388  }
389  }
390 
391  eps1min = Vmath::Vmin(nq, m_negepsvecminus1[0], 1);
392  eps2min = Vmath::Vmin(nq, m_negepsvecminus1[1], 1);
393  eps3min = Vmath::Vmin(nq, m_negepsvecminus1[2], 1);
394  eps1max = Vmath::Vmax(nq, m_negepsvecminus1[0], 1);
395  eps2max = Vmath::Vmax(nq, m_negepsvecminus1[1], 1);
396  eps3max = Vmath::Vmax(nq, m_negepsvecminus1[2], 1);
397 
398  mu1min = Vmath::Vmin(nq, m_negmuvecminus1[0], 1);
399  mu2min = Vmath::Vmin(nq, m_negmuvecminus1[1], 1);
400  mu3min = Vmath::Vmin(nq, m_negmuvecminus1[2], 1);
401  mu1max = Vmath::Vmax(nq, m_negmuvecminus1[0], 1);
402  mu2max = Vmath::Vmax(nq, m_negmuvecminus1[1], 1);
403  mu3max = Vmath::Vmax(nq, m_negmuvecminus1[2], 1);
404 
405  std::cout << "*** negepsvecminus1 = [ " << eps1min << " , " << eps1max
406  << " ], negepsvecminus1 = [ " << eps2min << " , " << eps2max
407  << " ], negepsvecminus1 = [ " << eps3min << " , " << eps3max
408  << " ] " << std::endl;
409  std::cout << "*** negmuvecminus1 = [ " << mu1min << " , " << mu1max
410  << " ], negmuvecminus1 = [ " << mu2min << " , " << mu2max
411  << " ], negmuvecminus1 = [ " << mu3min << " , " << mu3max << " ] "
412  << std::endl;
413 
414  // Compute de^m/dt \cdot e^k
415  if (m_AddRotation)
416  {
417  m_coriolis = Array<OneD, NekDouble>(nq);
419 
421  }
422 
423  // Generate Sigma Block with thicknes of m_PMLthickness and m_PMLmax
424  if (m_AddPML > 0)
425  {
427  }
428 
429  // If explicit it computes RHS and PROJECTION for the time integration
431  {
434  }
435  // Otherwise it gives an error (no implicit integration)
436  else
437  {
438  ASSERTL0(false, "Implicit unsteady Advection not set up.");
439  }
440 }
SourceType
Definition: MMFMaxwell.h:58
@ SIZE_SourceType
Definition: MMFMaxwell.h:62
CloakType
Definition: MMFMaxwell.h:42
@ eOpticalCloak
Definition: MMFMaxwell.h:44
@ eMicroWaveCloak
Definition: MMFMaxwell.h:47
@ eOpticalDispersiveCloak
Definition: MMFMaxwell.h:46
@ eOpticalConstCloak
Definition: MMFMaxwell.h:45
@ SIZE_CloakType
Definition: MMFMaxwell.h:48
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void ComputeMaterialVector(Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
Array< OneD, Array< OneD, NekDouble > > m_CrossProductMF
Definition: MMFMaxwell.h:125
void GenerateSigmaPML(const NekDouble PMLthickness, const NekDouble PMLstart, const NekDouble PMLmaxsigma, Array< OneD, Array< OneD, NekDouble >> &SigmaPML)
Array< OneD, NekDouble > ComputeRadCloak(const int Nlayer=5)
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the projection.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the RHS.
Definition: MMFMaxwell.cpp:887
void ComputeMaterialMicroWaveCloak(const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
void ComputeMaterialOpticalCloak(const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec, const bool Dispersion=false)
Array< OneD, NekDouble > EvaluateCoriolis()
SOLVER_UTILS_EXPORT void ComputeZimYim(Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
Definition: MMFSystem.cpp:1108
SOLVER_UTILS_EXPORT void DeriveCrossProductMF(Array< OneD, Array< OneD, NekDouble >> &CrossProductMF)
Definition: MMFSystem.cpp:571
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:53
SOLVER_UTILS_EXPORT void ComputeNtimesMF()
Definition: MMFSystem.cpp:618
SOLVER_UTILS_EXPORT void Computedemdxicdote()
Definition: MMFSystem.cpp:1280
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
@ SIZE_TestMaxwellType
Length of enum list.
Definition: MMFSystem.h:108
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:384

References ASSERTL0, CloakTypeMap, Nektar::SolverUtils::MMFSystem::Computedemdxicdote(), ComputeMaterialMicroWaveCloak(), ComputeMaterialOpticalCloak(), ComputeMaterialVector(), Nektar::SolverUtils::MMFSystem::ComputeNtimesMF(), ComputeRadCloak(), Nektar::SolverUtils::MMFSystem::ComputeZimYim(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineProjection(), Nektar::SolverUtils::MMFSystem::DeriveCrossProductMF(), DoOdeProjection(), DoOdeRhs(), eMicroWaveCloak, eOpticalCloak, eOpticalConstCloak, eOpticalDispersiveCloak, Nektar::SolverUtils::eTransElectric, Nektar::SolverUtils::eTransMagnetic, EvaluateCoriolis(), GenerateSigmaPML(), Nektar::SolverUtils::IncTypeMap, m_AddPML, m_AddRotation, m_boundaryforSF, m_CloakNlayer, m_Cloakraddelta, m_CloakType, m_coriolis, m_CrossProductMF, m_DispersiveCloak, m_ElemtGroup0, m_ElemtGroup1, Nektar::SolverUtils::MMFSystem::m_epsvec, Nektar::SolverUtils::UnsteadySystem::m_explicitAdvection, Nektar::SolverUtils::EquationSystem::m_fields, m_freq, m_Gaussianradius, Nektar::SolverUtils::MMFSystem::m_Incfreq, Nektar::SolverUtils::MMFSystem::m_IncType, m_mu, Nektar::SolverUtils::MMFSystem::m_muvec, m_n1, m_n2, m_n3, Nektar::SolverUtils::MMFSystem::m_negepsvecminus1, Nektar::SolverUtils::MMFSystem::m_negmuvecminus1, m_NoInc, Nektar::SolverUtils::UnsteadySystem::m_ode, m_PMLelement, m_PMLmaxsigma, m_PMLorder, m_PMLstart, m_PMLthickness, Nektar::SolverUtils::MMFSystem::m_PolType, m_PrintoutSurfaceCurrent, m_PSduration, m_Psx, m_Psy, m_Psz, m_RecPML, Nektar::SolverUtils::EquationSystem::m_session, m_SigmaPML, m_SourceType, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::MMFSystem::m_TestMaxwellType, m_TestPML, m_varepsilon, m_wp2, m_wp2Tol, Nektar::SolverUtils::MMFSystem::MMFInitObject(), Vmath::Neg(), Nektar::SolverUtils::PolTypeMap, Nektar::SolverUtils::MMFSystem::RootMeanSquare(), Vmath::Sadd(), SIZE_CloakType, Nektar::SolverUtils::SIZE_IncType, Nektar::SolverUtils::SIZE_PolType, SIZE_SourceType, Nektar::SolverUtils::SIZE_TestMaxwellType, Vmath::Smul(), SourceTypeMap, tinysimd::sqrt(), Nektar::SolverUtils::TestMaxwellTypeMap, Nektar::SolverUtils::UnsteadySystem::v_InitObject(), Vmath::Vmax(), and Vmath::Vmin().

◆ v_SetInitialConditions()

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

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

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1335 of file MMFMaxwell.cpp.

1338 {
1339  boost::ignore_unused(domain);
1340 
1341  int nq = GetTotPoints();
1342  int nvar = m_fields.size();
1343 
1344  switch (m_TestMaxwellType)
1345  {
1347  {
1348  m_fields[0]->SetPhys(TestMaxwell1D(initialtime, 0));
1349  m_fields[1]->SetPhys(TestMaxwell1D(initialtime, 1));
1350  }
1351  break;
1352 
1355  {
1356  m_fields[0]->SetPhys(TestMaxwell2DPEC(initialtime, 0, m_PolType));
1357  m_fields[1]->SetPhys(TestMaxwell2DPEC(initialtime, 1, m_PolType));
1358  m_fields[2]->SetPhys(TestMaxwell2DPEC(initialtime, 2, m_PolType));
1359  }
1360  break;
1361 
1363  {
1364  m_fields[0]->SetPhys(TestMaxwell2DPMC(initialtime, 0, m_PolType));
1365  m_fields[1]->SetPhys(TestMaxwell2DPMC(initialtime, 1, m_PolType));
1366  m_fields[2]->SetPhys(TestMaxwell2DPMC(initialtime, 2, m_PolType));
1367  }
1368  break;
1369 
1372  {
1373  Array<OneD, NekDouble> Zeros(nq, 0.0);
1374 
1375  for (int i = 0; i < nvar; i++)
1376  {
1377  m_fields[i]->SetPhys(Zeros);
1378  }
1379  }
1380  break;
1381 
1383  {
1384  m_fields[0]->SetPhys(TestMaxwellSphere(initialtime, m_freq, 0));
1385  m_fields[1]->SetPhys(TestMaxwellSphere(initialtime, m_freq, 1));
1386  m_fields[2]->SetPhys(TestMaxwellSphere(initialtime, m_freq, 2));
1387  }
1388  break;
1389 
1391  {
1392  m_fields[2]->SetPhys(GaussianPulse(initialtime, m_Psx, m_Psy, m_Psz,
1393  m_Gaussianradius));
1394  }
1395  break;
1396 
1397  default:
1398  break;
1399  }
1400 
1401  // forward transform to fill the modal coeffs
1402  for (int i = 0; i < nvar; ++i)
1403  {
1404  m_fields[i]->SetPhysState(true);
1405  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1406  m_fields[i]->UpdateCoeffs());
1407  }
1408 
1409  if (dumpInitialConditions)
1410  {
1411  std::string outname = m_sessionName + "_initial.chk";
1412  WriteFld(outname);
1413 
1414  Array<OneD, Array<OneD, NekDouble>> fields(nvar);
1415  for (int i = 0; i < nvar; ++i)
1416  {
1417  fields[i] = m_fields[i]->GetPhys();
1418  }
1419 
1420  Checkpoint_PlotOutput(0, fields);
1421  }
1422 }

References Checkpoint_PlotOutput(), Nektar::SolverUtils::eELF2DSurface, Nektar::SolverUtils::eMaxwell1D, Nektar::SolverUtils::eMaxwellSphere, Nektar::SolverUtils::eScatField2D, Nektar::SolverUtils::eTestMaxwell2DPEC, Nektar::SolverUtils::eTestMaxwell2DPECAVGFLUX, Nektar::SolverUtils::eTestMaxwell2DPMC, Nektar::SolverUtils::eTotField2D, GaussianPulse(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_freq, m_Gaussianradius, Nektar::SolverUtils::MMFSystem::m_PolType, m_Psx, m_Psy, m_Psz, Nektar::SolverUtils::EquationSystem::m_sessionName, Nektar::SolverUtils::MMFSystem::m_TestMaxwellType, TestMaxwell1D(), TestMaxwell2DPEC(), TestMaxwell2DPMC(), TestMaxwellSphere(), and Nektar::SolverUtils::EquationSystem::WriteFld().

◆ WeakDGMaxwellDirDeriv()

void Nektar::MMFMaxwell::WeakDGMaxwellDirDeriv ( const Array< OneD, const Array< OneD, NekDouble >> &  InField,
Array< OneD, Array< OneD, NekDouble >> &  OutField,
const NekDouble  time = 0.0 
)
protected

Calculate weak DG advection in the form \( \langle\phi, \hat{F}\cdot n\rangle - (\nabla \phi \cdot F) \).

Parameters
InFieldFields.
OutFieldStorage for result.
NumericalFluxIncludesNormalDefault: true.
InFieldIsPhysSpaceDefault: false.
nvariablesNumber of fields.

Definition at line 1246 of file MMFMaxwell.cpp.

1249 {
1250  int i;
1251  int nq = GetNpoints();
1252  int ncoeffs = GetNcoeffs();
1253  int nTracePointsTot = GetTraceNpoints();
1254  int nvar = 3;
1255 
1256  Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
1257  for (i = 0; i < m_shapedim; ++i)
1258  {
1259  fluxvector[i] = Array<OneD, NekDouble>(nq);
1260  }
1261 
1262  Array<OneD, Array<OneD, NekDouble>> physfield(nvar);
1263  for (i = 0; i < nvar; ++i)
1264  {
1265  physfield[i] = InField[i];
1266  }
1267 
1268  Array<OneD, NekDouble> tmpc(ncoeffs);
1269  for (i = 0; i < nvar; ++i)
1270  {
1271  GetMaxwellFluxVector(i, physfield, fluxvector);
1272 
1273  OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
1274  for (int j = 0; j < m_shapedim; ++j)
1275  {
1276  // Directional derivation with respect to the j'th moving frame
1277  // tmp_j = ( \nabla \phi, fluxvector[j] \mathbf{e}^j )
1278  m_fields[i]->IProductWRTDirectionalDerivBase(m_CrossProductMF[j],
1279  fluxvector[j], tmpc);
1280  Vmath::Vadd(ncoeffs, &tmpc[0], 1, &OutField[i][0], 1,
1281  &OutField[i][0], 1);
1282  }
1283  }
1284 
1285  // V the numerical flux and add to the modal coeffs
1286  // if the NumericalFlux function does not include the
1287  // normal in the output
1288  Array<OneD, Array<OneD, NekDouble>> numfluxFwd(nvar);
1289  Array<OneD, Array<OneD, NekDouble>> numfluxBwd(nvar);
1290 
1291  for (i = 0; i < nvar; ++i)
1292  {
1293  numfluxFwd[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
1294  numfluxBwd[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
1295  }
1296 
1297  // Evaluate numerical flux in physical space which may in
1298  // general couple all component of vectors
1299  NumericalMaxwellFlux(physfield, numfluxFwd, numfluxBwd, time);
1300 
1301  // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
1302  for (i = 0; i < nvar; ++i)
1303  {
1304  Vmath::Neg(ncoeffs, OutField[i], 1);
1305  m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
1306  OutField[i]);
1307  m_fields[i]->SetPhysState(false);
1308  }
1309 }
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)
Definition: MMFSystem.cpp:1730

References Nektar::SolverUtils::MMFSystem::GetMaxwellFluxVector(), Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::SolverUtils::EquationSystem::GetTraceNpoints(), m_CrossProductMF, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Neg(), Nektar::SolverUtils::MMFSystem::NumericalMaxwellFlux(), and Vmath::Vadd().

Referenced by DoOdeRhs().

Friends And Related Function Documentation

◆ MemoryManager< MMFMaxwell >

friend class MemoryManager< MMFMaxwell >
friend

Definition at line 65 of file MMFMaxwell.h.

Member Data Documentation

◆ className

std::string Nektar::MMFMaxwell::className
static
Initial value:
=
"MMFMaxwell", MMFMaxwell::create, "MMFMaxwell equation.")
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: MMFMaxwell.h:83
EquationSystemFactory & GetEquationSystemFactory()

Name of class.

Definition at line 93 of file MMFMaxwell.h.

◆ m_AddPML

int Nektar::MMFMaxwell::m_AddPML
protected

Definition at line 109 of file MMFMaxwell.h.

Referenced by DoOdeRhs(), GenerateSigmaPML(), v_GenerateSummary(), and v_InitObject().

◆ m_AddRotation

int Nektar::MMFMaxwell::m_AddRotation
protected

Definition at line 112 of file MMFMaxwell.h.

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

◆ m_boundaryforSF

int Nektar::MMFMaxwell::m_boundaryforSF
protected

Definition at line 106 of file MMFMaxwell.h.

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

◆ m_Cloaking

bool Nektar::MMFMaxwell::m_Cloaking
protected

Definition at line 114 of file MMFMaxwell.h.

◆ m_CloakNlayer

NekDouble Nektar::MMFMaxwell::m_CloakNlayer
protected

Definition at line 115 of file MMFMaxwell.h.

Referenced by v_GenerateSummary(), and v_InitObject().

◆ m_Cloakraddelta

NekDouble Nektar::MMFMaxwell::m_Cloakraddelta
protected

◆ m_CloakType

CloakType Nektar::MMFMaxwell::m_CloakType

Definition at line 78 of file MMFMaxwell.h.

Referenced by v_GenerateSummary(), and v_InitObject().

◆ m_coriolis

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

Definition at line 144 of file MMFMaxwell.h.

Referenced by AddCoriolis(), and v_InitObject().

◆ m_CrossProductMF

Array<OneD, Array<OneD, NekDouble> > Nektar::MMFMaxwell::m_CrossProductMF
protected

Definition at line 125 of file MMFMaxwell.h.

Referenced by v_InitObject(), and WeakDGMaxwellDirDeriv().

◆ m_DispersiveCloak

bool Nektar::MMFMaxwell::m_DispersiveCloak

◆ m_ElemtGroup0

int Nektar::MMFMaxwell::m_ElemtGroup0
protected

Definition at line 104 of file MMFMaxwell.h.

Referenced by ComputeMaterialMicroWaveCloak(), and v_InitObject().

◆ m_ElemtGroup1

int Nektar::MMFMaxwell::m_ElemtGroup1
protected

◆ m_freq

NekDouble Nektar::MMFMaxwell::m_freq
protected

◆ m_Gaussianradius

NekDouble Nektar::MMFMaxwell::m_Gaussianradius
protected

Definition at line 123 of file MMFMaxwell.h.

Referenced by v_DoSolve(), v_GenerateSummary(), v_InitObject(), and v_SetInitialConditions().

◆ m_mu

Array<OneD, NekDouble> Nektar::MMFMaxwell::m_mu
protected

Definition at line 135 of file MMFMaxwell.h.

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

◆ m_n1

NekDouble Nektar::MMFMaxwell::m_n1
protected

Definition at line 133 of file MMFMaxwell.h.

Referenced by TestMaxwell1D(), and v_InitObject().

◆ m_n2

NekDouble Nektar::MMFMaxwell::m_n2
protected

Definition at line 133 of file MMFMaxwell.h.

Referenced by TestMaxwell1D(), and v_InitObject().

◆ m_n3

NekDouble Nektar::MMFMaxwell::m_n3
protected

Definition at line 133 of file MMFMaxwell.h.

Referenced by v_InitObject().

◆ m_NoInc

int Nektar::MMFMaxwell::m_NoInc
protected

Definition at line 142 of file MMFMaxwell.h.

Referenced by v_InitObject().

◆ m_PMLelement

int Nektar::MMFMaxwell::m_PMLelement
protected

Definition at line 138 of file MMFMaxwell.h.

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

◆ m_PMLmaxsigma

NekDouble Nektar::MMFMaxwell::m_PMLmaxsigma
protected

Definition at line 139 of file MMFMaxwell.h.

Referenced by v_GenerateSummary(), and v_InitObject().

◆ m_PMLorder

int Nektar::MMFMaxwell::m_PMLorder
protected

Definition at line 110 of file MMFMaxwell.h.

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

◆ m_PMLstart

NekDouble Nektar::MMFMaxwell::m_PMLstart
protected

Definition at line 139 of file MMFMaxwell.h.

Referenced by v_GenerateSummary(), and v_InitObject().

◆ m_PMLthickness

NekDouble Nektar::MMFMaxwell::m_PMLthickness
protected

Definition at line 139 of file MMFMaxwell.h.

Referenced by v_GenerateSummary(), and v_InitObject().

◆ m_PrintoutSurfaceCurrent

int Nektar::MMFMaxwell::m_PrintoutSurfaceCurrent
protected

Definition at line 107 of file MMFMaxwell.h.

Referenced by v_DoSolve(), and v_InitObject().

◆ m_PSduration

NekDouble Nektar::MMFMaxwell::m_PSduration
protected

Definition at line 123 of file MMFMaxwell.h.

Referenced by GaussianPulse(), v_DoSolve(), v_GenerateSummary(), and v_InitObject().

◆ m_Psx

NekDouble Nektar::MMFMaxwell::m_Psx
protected

Definition at line 122 of file MMFMaxwell.h.

Referenced by v_DoSolve(), v_GenerateSummary(), v_InitObject(), and v_SetInitialConditions().

◆ m_Psy

NekDouble Nektar::MMFMaxwell::m_Psy
protected

Definition at line 122 of file MMFMaxwell.h.

Referenced by v_DoSolve(), v_GenerateSummary(), v_InitObject(), and v_SetInitialConditions().

◆ m_Psz

NekDouble Nektar::MMFMaxwell::m_Psz
protected

Definition at line 122 of file MMFMaxwell.h.

Referenced by v_DoSolve(), v_GenerateSummary(), v_InitObject(), and v_SetInitialConditions().

◆ m_RecPML

int Nektar::MMFMaxwell::m_RecPML
protected

Definition at line 138 of file MMFMaxwell.h.

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

◆ m_SigmaPML

Array<OneD, Array<OneD, NekDouble> > Nektar::MMFMaxwell::m_SigmaPML
protected

Definition at line 140 of file MMFMaxwell.h.

Referenced by AddPML(), and v_InitObject().

◆ m_SourceType

SourceType Nektar::MMFMaxwell::m_SourceType

Definition at line 79 of file MMFMaxwell.h.

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

◆ m_SourceVector

Array<OneD, NekDouble> Nektar::MMFMaxwell::m_SourceVector
protected

Definition at line 121 of file MMFMaxwell.h.

Referenced by v_DoSolve().

◆ m_TestPML

int Nektar::MMFMaxwell::m_TestPML
protected

Definition at line 137 of file MMFMaxwell.h.

Referenced by v_DoSolve(), and v_InitObject().

◆ m_varepsilon

Array<OneD, NekDouble> Nektar::MMFMaxwell::m_varepsilon
protected

Definition at line 134 of file MMFMaxwell.h.

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

◆ m_wp2

Array<OneD, NekDouble> Nektar::MMFMaxwell::m_wp2
protected

Definition at line 119 of file MMFMaxwell.h.

Referenced by AddPML(), ComputeMaterialOpticalCloak(), DoOdeRhs(), and v_InitObject().

◆ m_wp2Tol

NekDouble Nektar::MMFMaxwell::m_wp2Tol
protected

Definition at line 118 of file MMFMaxwell.h.

Referenced by ComputeMaterialOpticalCloak(), and v_InitObject().