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) override
 Initialise the object. More...
 
virtual void v_DoSolve () override
 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 ()
 
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &s) override
 Print a summary of time stepping parameters. More...
 
SOLVER_UTILS_EXPORT void MMFInitObject (const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
 
SOLVER_UTILS_EXPORT void CopyBoundaryTrace (const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const BoundaryCopyType BDCopyType, const int var=0, const std::string btype="NoUserDefined")
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
virtual SOLVER_UTILS_EXPORT ~UnsteadySystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Calculate the larger time-step mantaining the problem stable. More...
 
SOLVER_UTILS_EXPORT void SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeSharedPtrGetTimeIntegrationScheme ()
 Returns the time integration scheme. More...
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeOperatorsGetTimeIntegrationSchemeOperators ()
 Returns the time integration scheme operators. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void InitObject (bool DeclareField=true)
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise (bool dumpInitialConditions=true)
 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 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 NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDoubleErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write checkpoint file of custom data fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write input fields to the given filename. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
 Input field data from the given file to multiple domains. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, std::vector< std::string > &fieldStr, Array< OneD, Array< OneD, NekDouble > > &coeffs)
 Output a field. Input field data into array from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, MultiRegions::ExpListSharedPtr &pField, std::string &pFieldName)
 Output a field. Input field data into ExpList from the given file. More...
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary. More...
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMapUpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetFinalTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, const Array< OneD, const NekDouble > &input)
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT int GetInfoSteps ()
 
SOLVER_UTILS_EXPORT void SetInfoSteps (int num)
 
SOLVER_UTILS_EXPORT void SetIterationNumberPIT (int num)
 
SOLVER_UTILS_EXPORT void SetWindowNumberPIT (int num)
 
SOLVER_UTILS_EXPORT Array< OneD, const Array< OneD, NekDouble > > GetTraceNormals ()
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
SOLVER_UTILS_EXPORT bool NegatedOp ()
 Identify if operator is negated in DoSolve. More...
 
SOLVER_UTILS_EXPORT bool ParallelInTime ()
 Check if solver use Parallel-in-Time. 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
 

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) override
 Print Summary. More...
 
virtual void v_SetInitialConditions (const NekDouble initialtime, bool dumpInitialConditions, const int domain) override
 
virtual void v_EvaluateExactSolution (unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
 
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...
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoSolve () override
 Solves an unsteady problem. More...
 
virtual SOLVER_UTILS_EXPORT void v_PrintStatusInformation (const int step, const NekDouble cpuTime)
 Print Status Information. More...
 
virtual SOLVER_UTILS_EXPORT void v_PrintSummaryStatistics (const NekDouble intTime)
 Print Summary Statistics. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true) override
 Sets up initial conditions. More...
 
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &s) override
 Print a summary of time stepping parameters. 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)
 
virtual SOLVER_UTILS_EXPORT bool v_UpdateTimeStepCheck ()
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. More...
 
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...
 
SOLVER_UTILS_EXPORT void DoDummyProjection (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
 Perform dummy projection. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members. More...
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareFeld=true)
 Initialisation object for EquationSystem. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true)
 Virtual function for initialisation implementation. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoSolve ()
 Virtual function for solve implementation. 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_GenerateSummary (SummaryList &l)
 Virtual function for generating summary information. More...
 
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution (unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 
virtual SOLVER_UTILS_EXPORT void v_Output (void)
 
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure (void)
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp (void)
 Virtual function to identify if operator is negated in DoSolve. More...
 
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
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check. More...
 
std::vector< int > m_intVariables
 
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 
NekDouble m_CFLGrowth
 CFL growth rate. More...
 
NekDouble m_CFLEnd
 Maximun cfl in cfl growth. More...
 
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
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...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
std::ofstream m_errFile
 
NekDouble m_epsilon
 Diffusion coefficient. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
bool m_verbose
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtrm_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_fields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
NekDouble m_lastCheckTime
 
NekDouble m_TimeIncrementFactor
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_infosteps
 Number of time steps between outputting status information. More...
 
int m_iterPIT = 0
 Number of parallel-in-time time iteration. More...
 
int m_windowPIT = 0
 Index of windows for parallel-in-time time iteration. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_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 []
 
static std::string projectionTypeLookupIds []
 

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 2869 of file MMFMaxwell.cpp.

2871{
2872 int nq = physarray[0].size();
2873
2874 Array<OneD, NekDouble> tmp(nq);
2875
2876 int indx;
2877 for (int j = 0; j < m_shapedim; ++j)
2878 {
2879 if (j == 0)
2880 {
2881 indx = 2;
2882 }
2883
2884 else if (j == 1)
2885 {
2886 indx = 1;
2887 }
2888
2889 Vmath::Vmul(nq, m_coriolis, 1, physarray[indx], 1, tmp, 1);
2890
2891 switch (m_PolType)
2892 {
2894 {
2895 Vmath::Vmul(nq, m_muvec[indx], 1, tmp, 1, tmp, 1);
2896 }
2897 break;
2898
2900 {
2901 Vmath::Vmul(nq, m_epsvec[indx], 1, tmp, 1, tmp, 1);
2902 }
2903 break;
2904
2905 default:
2906 break;
2907 }
2908
2909 if (j == 1)
2910 {
2911 Vmath::Neg(nq, tmp, 1);
2912 }
2913
2914 Vmath::Vadd(nq, tmp, 1, outarray[j], 1, outarray[j], 1);
2915 }
2916}
Array< OneD, NekDouble > m_coriolis
Definition: MMFMaxwell.h:144
Array< OneD, Array< OneD, NekDouble > > m_muvec
Definition: MMFSystem.h:212
Array< OneD, Array< OneD, NekDouble > > m_epsvec
Definition: MMFSystem.h:211
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:207
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:513
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:354

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 1171 of file MMFMaxwell.cpp.

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

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 2425 of file MMFMaxwell.cpp.

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

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 2684 of file MMFMaxwell.cpp.

2687{
2688 boost::ignore_unused(time);
2689
2690 int nvar = m_fields.size();
2691 int nq = m_fields[0]->GetTotPoints();
2692 int ncoeffs = m_fields[0]->GetNcoeffs();
2693
2694 std::string outname = m_sessionName + "EDFlux_" +
2695 boost::lexical_cast<std::string>(n) + ".chk";
2696
2697 std::vector<std::string> variables(nvar);
2698 variables[0] = "EDFx";
2699 variables[1] = "EDFy";
2700 variables[2] = "EDFz";
2701 for (int i = 3; i < nvar; ++i)
2702 {
2703 variables[i] = m_boundaryConditions->GetVariable(i);
2704 }
2705
2706 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2707 for (int i = 0; i < nvar; ++i)
2708 {
2709 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2710 }
2711
2712 Array<OneD, NekDouble> tmp(nq);
2713
2714 // TE: H^3 (E^2 e^1 - E^1 e^2 )
2715 // TM: -E^3 (H^2 e^1 - H^1 e^2 )
2716 for (int k = 0; k < m_spacedim; ++k)
2717 {
2718 Vmath::Vmul(nq, &fieldphys[0][0], 1, &m_movingframes[1][k * nq], 1,
2719 &tmp[0], 1);
2720 Vmath::Vvtvm(nq, &fieldphys[1][0], 1, &m_movingframes[0][k * nq], 1,
2721 &tmp[0], 1, &tmp[0], 1);
2722
2723 Vmath::Vmul(nq, &fieldphys[2][0], 1, &tmp[0], 1, &tmp[0], 1);
2724
2726 {
2727 Vmath::Neg(nq, tmp, 1);
2728 }
2729
2730 m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2731 }
2732
2733 WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2734}
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:185

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 2736 of file MMFMaxwell.cpp.

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

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 2595 of file MMFMaxwell.cpp.

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

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 2562 of file MMFMaxwell.cpp.

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

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 2630 of file MMFMaxwell.cpp.

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

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 2185 of file MMFMaxwell.cpp.

2187{
2188 int nq = GetTotPoints();
2189 NekDouble energy;
2190
2191 Array<OneD, NekDouble> tmp(nq, 0.0);
2192
2193 for (int i = 0; i < 3; ++i)
2194 {
2195 Vmath::Vvtvp(nq, &fields[i][0], 1, &fields[i][0], 1, &tmp[0], 1,
2196 &tmp[0], 1);
2197 }
2198
2199 energy = 0.5 * (m_fields[0]->Integral(tmp));
2200 return energy;
2201}
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 2371 of file MMFMaxwell.cpp.

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

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

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 2203 of file MMFMaxwell.cpp.

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

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

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 1984 of file MMFMaxwell.cpp.

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

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 1320 of file MMFMaxwell.cpp.

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

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 889 of file MMFMaxwell.cpp.

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

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 2838 of file MMFMaxwell.cpp.

2839{
2840 int nq = GetTotPoints();
2841
2842 NekDouble m_Omega = 1.5486 * 0.000001;
2843
2844 NekDouble x0j, x1j, x2j;
2845 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2846
2847 Array<OneD, NekDouble> x(nq);
2848 Array<OneD, NekDouble> y(nq);
2849 Array<OneD, NekDouble> z(nq);
2850
2851 m_fields[0]->GetCoords(x, y, z);
2852
2853 Array<OneD, NekDouble> outarray(nq, 0.0);
2854 for (int j = 0; j < nq; ++j)
2855 {
2856 x0j = x[j];
2857 x1j = y[j];
2858 x2j = z[j];
2859
2860 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2861 cos_theta);
2862
2863 outarray[j] = 2.0 * m_Omega * sin_theta;
2864 }
2865
2866 return outarray;
2867}
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:778
std::vector< double > z(NPUPPER)

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

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 2801 of file MMFMaxwell.cpp.

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

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

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 2029 of file MMFMaxwell.cpp.

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

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(), Vmath::Vsub(), and Nektar::UnitTests::z().

Referenced by v_InitObject().

◆ print_MMF()

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

Definition at line 3138 of file MMFMaxwell.cpp.

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

References tinysimd::sqrt().

◆ Printout_SurfaceCurrent()

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

Definition at line 1906 of file MMFMaxwell.cpp.

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

Referenced by v_DoSolve().

◆ TestMaxwell1D()

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

Definition at line 1469 of file MMFMaxwell.cpp.

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

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 1575 of file MMFMaxwell.cpp.

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

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 1706 of file MMFMaxwell.cpp.

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

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 1792 of file MMFMaxwell.cpp.

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

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

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

◆ v_DoSolve()

void Nektar::MMFMaxwell::v_DoSolve ( void  )
overridevirtual

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_checksteps ? m_steps / m_checksteps + 1 : 0;
493
494 Array<OneD, NekDouble> TimeSeries(Ntot ? Ntot : 1);
495 Array<OneD, NekDouble> Energy(Ntot ? Ntot : 1);
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>(
514
515 Array<OneD, NekDouble> x(nq);
516 Array<OneD, NekDouble> y(nq);
517 Array<OneD, NekDouble> z(nq);
518
519 m_fields[0]->GetCoords(x, y, z);
520
521 NekDouble Tol = 0.000001;
523 for (i = 0; i < nq; ++i)
524 {
525 rad = sqrt((x[i] + m_Psx) * (x[i] + m_Psx) +
526 (y[i] + m_Psy) * (y[i] + m_Psy) +
527 (z[i] + m_Psz) * (z[i] + m_Psz));
528 std::cout << "rad" << rad << std::endl;
529 if (rad < Tol)
530 {
531 indxantipod = i;
532 break;
533 }
534 }
535 }
536 break;
537
538 case ePlanarSource:
539 {
540 m_SourceVector = Array<OneD, NekDouble>(nq, 0.0);
541
542 Array<OneD, NekDouble> x(nq);
543 Array<OneD, NekDouble> y(nq);
544 Array<OneD, NekDouble> z(nq);
545
546 m_fields[0]->GetCoords(x, y, z);
547
548 NekDouble Tol = 0.000001;
550 for (i = 0; i < nq; ++i)
551 {
552 rad = sqrt((x[i] - m_Psx) * (x[i] - m_Psx));
553 if (rad < Tol)
554 {
555 m_SourceVector[i] = 1.0;
556 }
557 }
558
559 std::cout << "*** Area of Planar Source = "
560 << m_fields[0]->Integral(m_SourceVector) << std::endl;
561 }
562 break;
563
564 default:
565 break;
566 }
567
568 int cntpml = 0;
569 int P1indx = 0, P2indx = 0, P3indx = 0;
570 Array<OneD, NekDouble> P1;
571 Array<OneD, NekDouble> P2;
572 Array<OneD, NekDouble> P3;
573 if (m_TestPML)
574 {
575 P1 = Array<OneD, NekDouble>(m_checksteps ? m_steps / m_checksteps : 1);
576 P2 = Array<OneD, NekDouble>(m_checksteps ? m_steps / m_checksteps : 1);
577 P3 = Array<OneD, NekDouble>(m_checksteps ? m_steps / m_checksteps : 1);
578
579 Array<OneD, NekDouble> x(nq);
580 Array<OneD, NekDouble> y(nq);
581 Array<OneD, NekDouble> z(nq);
582
583 m_fields[0]->GetCoords(x, y, z);
584
585 NekDouble Tol = 0.000001;
587 for (int i = 0; i < nq; ++i)
588 {
589 rad = sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i]) * (y[i]));
590
591 if (rad < Tol)
592 {
593 P1indx = i;
594 break;
595 }
596 }
597
598 for (int i = 0; i < nq; ++i)
599 {
600 rad =
601 sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i] - 1.5) * (y[i] - 1.5));
602 if (rad < Tol)
603 {
604 P2indx = i;
605 break;
606 }
607 }
608
609 for (int i = 0; i < nq; ++i)
610 {
611 rad =
612 sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i] - 3.0) * (y[i] - 3.0));
613 if (rad < Tol)
614 {
615 P3indx = i;
616 break;
617 }
618 }
619 }
620
621 int indx;
623 {
624
625 timer.Start();
626 fields = m_intScheme->TimeIntegrate(step, m_timestep);
627 timer.Stop();
628
630 elapsed = timer.TimePerTest(1);
631 intTime += elapsed;
632 cpuTime += elapsed;
633
634 // Write out status information
635 if (m_infosteps && !((step + 1) % m_infosteps) &&
636 m_session->GetComm()->GetRank() == 0)
637 {
638 std::cout << "Steps: " << std::setw(8) << std::left << step + 1
639 << " "
640 << "Time: " << std::setw(12) << std::left << m_time;
641
642 std::stringstream ss;
643 ss << cpuTime / 60.0 << " min.";
644 std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
645 << std::endl;
646
647 cpuTime = 0.0;
648 }
649
650 switch (m_SourceType)
651 {
652 case ePointSource:
653 {
654 if (m_time <= m_PSduration)
655 {
656 Array<OneD, NekDouble> Impulse(nq);
657 Impulse = GaussianPulse(m_time, m_Psx, m_Psy, m_Psz,
659 Vmath::Vadd(nq, &Impulse[0], 1,
660 &fields[m_intVariables[2]][0], 1,
661 &fields[m_intVariables[2]][0], 1);
662 }
663 }
664 break;
665
666 case ePlanarSource:
667 {
668 Array<OneD, NekDouble> Impulse(nq);
669 for (int i = 0; i < 3; ++i)
670 {
671 Impulse = GetIncidentField(i, m_time);
672 Vmath::Vmul(nq, m_SourceVector, 1, Impulse, 1, Impulse, 1);
673 Vmath::Vadd(nq, &Impulse[0], 1,
674 &fields[m_intVariables[i]][0], 1,
675 &fields[m_intVariables[i]][0], 1);
676 }
677 }
678 break;
679
680 default:
681 break;
682 }
683
684 // Transform data into coefficient space
685 for (i = 0; i < nvariables; ++i)
686 {
687 m_fields[m_intVariables[i]]->SetPhys(fields[i]);
688 m_fields[m_intVariables[i]]->FwdTransLocalElmt(
689 fields[i], m_fields[m_intVariables[i]]->UpdateCoeffs());
690 m_fields[m_intVariables[i]]->SetPhysState(false);
691 }
692 // for (i = 0; i < nq; ++i)
693 // std::cout << m_fields[0][0][i] <<std::endl;
694
695 // Write out checkpoint files
696 if ((m_checksteps && step && !((step + 1) % m_checksteps)) ||
697 doCheckTime)
698 {
699 indx = (step + 1) / m_checksteps;
700 TimeSeries[indx] = m_time;
701
703 {
704 Checkpoint_TotalFieldOutput(nchk, m_time, fields);
705 Checkpoint_TotPlotOutput(nchk, m_time, fields);
706 }
707 Checkpoint_PlotOutput(nchk, fields);
708 Checkpoint_EDFluxOutput(nchk, m_time, fields);
709 Checkpoint_EnergyOutput(nchk, m_time, fields);
710 Checkpoint_Output(nchk++);
711
712 Energy[indx] = ComputeEnergyDensity(fields);
713
714 std::cout << "|EHr|: F1 = " << RootMeanSquare(fields[0])
715 << ", F2 = " << RootMeanSquare(fields[1])
716 << ", F3 = " << RootMeanSquare(fields[2])
717 << ", Energy = " << Energy[indx] << std::endl;
718 if (nfields > 3)
719 {
720 std::cout << "|DBr|: D1 = " << RootMeanSquare(fields[3])
721 << ", D2 = " << RootMeanSquare(fields[4])
722 << ", D3 = " << RootMeanSquare(fields[5])
723 << std::endl;
724
725 int nTraceNumPoints = GetTraceNpoints();
726 int totbdryexp =
727 m_fields[0]->GetBndCondExpansions()[0]->GetExpSize();
728 int npts = m_fields[0]
729 ->GetBndCondExpansions()[0]
730 ->GetExp(0)
731 ->GetNumPoints(0);
732
733 Array<OneD, NekDouble> x0(nq);
734 Array<OneD, NekDouble> x1(nq);
735 Array<OneD, NekDouble> x2(nq);
736
737 m_fields[0]->GetCoords(x0, x1, x2);
738
739 Array<OneD, NekDouble> E1Fwd(nTraceNumPoints);
740 Array<OneD, NekDouble> E2Fwd(nTraceNumPoints);
741 Array<OneD, NekDouble> H3Fwd(nTraceNumPoints);
742
743 m_fields[0]->ExtractTracePhys(fields[0], E1Fwd);
744 m_fields[0]->ExtractTracePhys(fields[1], E2Fwd);
745 m_fields[0]->ExtractTracePhys(fields[2], H3Fwd);
746
747 int id2, cnt = 0;
748 NekDouble E1atPECloc, E1atPEC = 0.0;
749 NekDouble E2atPECloc, E2atPEC = 0.0;
750 NekDouble H3atPECloc, H3atPEC = 0.0;
751
752 Array<OneD, NekDouble> E1Fwdloc(npts);
753 Array<OneD, NekDouble> E2Fwdloc(npts);
754 Array<OneD, NekDouble> H3Fwdloc(npts);
755
756 for (int e = 0; e < totbdryexp; ++e)
757 {
758 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
759 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
760 cnt + e));
761
762 Vmath::Vcopy(npts, &E1Fwd[id2], 1, &E1Fwdloc[0], 1);
763 Vmath::Vcopy(npts, &E2Fwd[id2], 1, &E2Fwdloc[0], 1);
764 Vmath::Vcopy(npts, &H3Fwd[id2], 1, &H3Fwdloc[0], 1);
765
766 E1atPECloc = Vmath::Vamax(npts, E1Fwdloc, 1);
767 E2atPECloc = Vmath::Vamax(npts, E2Fwdloc, 1);
768 H3atPECloc = Vmath::Vamax(npts, H3Fwdloc, 1);
769
770 if (E1atPEC < E1atPECloc)
771 {
772 E1atPEC = E1atPECloc;
773 }
774
775 if (E2atPEC < E2atPECloc)
776 {
777 E2atPEC = E2atPECloc;
778 }
779
780 if (H3atPEC < H3atPECloc)
781 {
782 H3atPEC = H3atPECloc;
783 }
784 }
785
786 std::cout << "At PEC, Max. E1 = " << E1atPEC
787 << ", E2 = " << E2atPEC << ", H3 = " << H3atPEC
788 << std::endl;
789 }
790
792 {
793 Ezantipod[cntap++] = fields[2][indxantipod];
794 }
795
796 if (m_TestPML)
797 {
798 P1[cntpml] = fields[2][P1indx];
799 P2[cntpml] = fields[2][P2indx];
800 P3[cntpml] = fields[2][P3indx];
801 cntpml++;
802 }
803 doCheckTime = false;
804 }
805
806 // Step advance
807 ++step;
808 }
809
810 // Print out summary statistics
811 if (m_checksteps && m_session->GetComm()->GetRank() == 0)
812 {
813 std::cout << "Time-integration : " << intTime << "s" << std::endl;
814
815 std::cout << "TimeSeries = " << std::endl;
816 for (int i = 0; i < m_steps / m_checksteps; ++i)
817 {
818 std::cout << TimeSeries[i] << ", ";
819 }
820 std::cout << std::endl << std::endl;
821
822 std::cout << "Energy Density = " << std::endl;
823 for (int i = 0; i < m_steps / m_checksteps; ++i)
824 {
825 std::cout << Energy[i] << ", ";
826 }
827 std::cout << std::endl << std::endl;
828
830 {
832 }
833
835 {
836 std::cout << "Ez at antipod = " << std::endl;
837 for (int i = 0; i < m_steps / m_checksteps; ++i)
838 {
839 std::cout << Ezantipod[i] << ", ";
840 }
841 std::cout << std::endl << std::endl;
842 }
843
844 if (m_TestPML)
845 {
846 std::cout << "P1 = " << std::endl;
847 for (int i = 0; i < m_steps / m_checksteps; ++i)
848 {
849 std::cout << P1[i] << ", ";
850 }
851 std::cout << std::endl << std::endl;
852
853 std::cout << "P2 = " << std::endl;
854 for (int i = 0; i < m_steps / m_checksteps; ++i)
855 {
856 std::cout << P2[i] << ", ";
857 }
858 std::cout << std::endl << std::endl;
859
860 std::cout << "P3 = " << std::endl;
861 for (int i = 0; i < m_steps / m_checksteps; ++i)
862 {
863 std::cout << P3[i] << ", ";
864 }
865 std::cout << std::endl << std::endl;
866 }
867 }
868
869 for (i = 0; i < nvariables; ++i)
870 {
871 m_fields[m_intVariables[i]]->SetPhys(fields[i]);
872 m_fields[m_intVariables[i]]->SetPhysState(true);
873 }
874
875 for (i = 0; i < nvariables; ++i)
876 {
877 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
878 m_fields[i]->UpdateCoeffs());
879 }
880}
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
@ ePointSource
Definition: MMFMaxwell.h:60
@ ePlanarSource
Definition: MMFMaxwell.h:61
void Checkpoint_EDFluxOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
NekDouble m_Gaussianradius
Definition: MMFMaxwell.h:123
Array< OneD, NekDouble > m_SourceVector
Definition: MMFMaxwell.h:121
NekDouble ComputeEnergyDensity(Array< OneD, Array< OneD, NekDouble > > &fields)
void Checkpoint_EnergyOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
void Checkpoint_TotalFieldOutput(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)
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 Printout_SurfaceCurrent(Array< OneD, Array< OneD, NekDouble > > &fields, const int time)
void Checkpoint_PlotOutput(const int n, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
SOLVER_UTILS_EXPORT int GetTraceNpoints()
NekDouble m_timestep
Time step size.
int m_infosteps
Number of time steps between outputting status information.
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:2327
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
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:994

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::EquationSystem::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(), Vmath::Vmul(), and Nektar::UnitTests::z().

◆ v_EvaluateExactSolution()

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1427 of file MMFMaxwell.cpp.

1430{
1431 int nq = m_fields[0]->GetNpoints();
1432 outfield = Array<OneD, NekDouble>(nq);
1433
1434 switch (m_TestMaxwellType)
1435 {
1437 {
1438 outfield = TestMaxwell1D(time, field);
1439 }
1440 break;
1441
1444 {
1445 outfield = TestMaxwell2DPEC(time, field, m_PolType);
1446 }
1447 break;
1448
1450 {
1451 outfield = TestMaxwell2DPMC(time, field, m_PolType);
1452 }
1453 break;
1454
1456 {
1457 outfield = TestMaxwellSphere(time, m_freq, field);
1458 }
1459 break;
1460
1461 default:
1462 {
1463 outfield = Array<OneD, NekDouble>(nq, 0.0);
1464 }
1465 break;
1466 }
1467}
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)
overrideprotectedvirtual

Print Summary.

Reimplemented from Nektar::SolverUtils::MMFSystem.

Definition at line 3070 of file MMFMaxwell.cpp.

3071{
3074 s, "TestMaxwellType",
3076 SolverUtils::AddSummaryItem(s, "PolType",
3078 SolverUtils::AddSummaryItem(s, "IncType",
3080
3081 if (m_varepsilon[0] * m_varepsilon[1] * m_varepsilon[2] > 1.0)
3082 {
3083 SolverUtils::AddSummaryItem(s, "varepsilon1", m_varepsilon[0]);
3084 SolverUtils::AddSummaryItem(s, "varepsilon2", m_varepsilon[1]);
3085 SolverUtils::AddSummaryItem(s, "varepsilon3", m_varepsilon[2]);
3086 }
3087
3088 if (m_mu[0] * m_mu[1] * m_mu[2] > 1.0)
3089 {
3090 SolverUtils::AddSummaryItem(s, "mu1", m_mu[0]);
3091 SolverUtils::AddSummaryItem(s, "mu2", m_mu[1]);
3092 SolverUtils::AddSummaryItem(s, "mu3", m_mu[2]);
3093 }
3094
3095 if (m_boundaryforSF > 0)
3096 {
3098 }
3099
3100 if (m_ElemtGroup1 > 0)
3101 {
3102 SolverUtils::AddSummaryItem(s, "CloakNlayer", m_CloakNlayer);
3103 SolverUtils::AddSummaryItem(s, "ElemtGroup1", m_ElemtGroup1);
3104 }
3105
3106 SolverUtils::AddSummaryItem(s, "AddRotation", m_AddRotation);
3107
3108 if (m_AddPML > 0)
3109 {
3111 SolverUtils::AddSummaryItem(s, "PMLelement", m_PMLelement);
3114 SolverUtils::AddSummaryItem(s, "PMLthickness", m_PMLthickness);
3116 SolverUtils::AddSummaryItem(s, "PMLmaxsigma", m_PMLmaxsigma);
3117 }
3118
3119 if (m_SourceType)
3120 {
3121 SolverUtils::AddSummaryItem(s, "SourceType",
3126 SolverUtils::AddSummaryItem(s, "PSduration", m_PSduration);
3127 SolverUtils::AddSummaryItem(s, "Gaussianradius", m_Gaussianradius);
3128 }
3129
3130 if (m_CloakType)
3131 {
3133 SolverUtils::AddSummaryItem(s, "DispersiveCloak", m_DispersiveCloak);
3134 SolverUtils::AddSummaryItem(s, "CloakNlayer", m_CloakNlayer);
3135 SolverUtils::AddSummaryItem(s, "Cloakraddelta", m_Cloakraddelta);
3136 }
3137}
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) override
Print a summary of time stepping parameters.
Definition: MMFSystem.cpp:2452
const char *const IncTypeMap[]
Definition: MMFSystem.h:138
const char *const TestMaxwellTypeMap[]
Definition: MMFSystem.h:110
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:125

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

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
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 {
203 break;
204 }
205 }
206 }
207 else
208 {
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
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
254 {
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
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 DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
Definition: MMFMaxwell.cpp:889
Array< OneD, Array< OneD, NekDouble > > m_CrossProductMF
Definition: MMFMaxwell.h:125
Array< OneD, NekDouble > ComputeRadCloak(const int Nlayer=5)
void ComputeMaterialMicroWaveCloak(const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec)
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
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 ComputeMaterialVector(Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec)
void GenerateSigmaPML(const NekDouble PMLthickness, const NekDouble PMLstart, const NekDouble PMLmaxsigma, Array< OneD, Array< OneD, NekDouble > > &SigmaPML)
Array< OneD, NekDouble > EvaluateCoriolis()
SOLVER_UTILS_EXPORT void DeriveCrossProductMF(Array< OneD, Array< OneD, NekDouble > > &CrossProductMF)
Definition: MMFSystem.cpp:554
SOLVER_UTILS_EXPORT void ComputeZimYim(Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec)
Definition: MMFSystem.cpp:1091
SOLVER_UTILS_EXPORT void ComputeNtimesMF()
Definition: MMFSystem.cpp:601
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 Computedemdxicdote()
Definition: MMFSystem.cpp:1263
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
@ SIZE_TestMaxwellType
Length of enum list.
Definition: MMFSystem.h:107
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add scalar y = alpha + x.
Definition: Vmath.cpp:379

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

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 1338 of file MMFMaxwell.cpp.

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

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 1248 of file MMFMaxwell.cpp.

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

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().