39#include <boost/algorithm/string.hpp>
52 "Nonlinear Peregrine equations in conservative variables.");
69 if (
m_session->DefinesSolverInfo(
"PROBLEMTYPE"))
72 std::string ProblemTypeStr =
m_session->GetSolverInfo(
"PROBLEMTYPE");
94 ASSERTL0(
false,
"Implicit Peregrine not set up.");
99 if (
m_session->DefinesParameter(
"ConstDepth"))
105 ASSERTL0(
false,
"Constant Depth not specified");
115 "Continuous projection type not supported for Peregrine.");
129 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
138 m_session->LoadSolverInfo(
"UpwindType", riemName,
"NoSolver");
161 ASSERTL0(
false,
"Unsupported projection type.");
189 m_fields[0]->IProductWRTBase(tmp, mod);
190 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
192 Vmath::Vadd(nq, tmp, 1, outarray[1], 1, outarray[1], 1);
197 m_fields[0]->IProductWRTBase(tmp, mod);
198 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
200 Vmath::Vadd(nq, tmp, 1, outarray[2], 1, outarray[2], 1);
208 Vmath::Vadd(nq, tmp, 1, outarray[1], 1, outarray[1], 1);
213 Vmath::Vadd(nq, tmp, 1, outarray[2], 1, outarray[2], 1);
217 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
241 m_fields[0]->IProductWRTBase(tmp, mod);
242 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
244 Vmath::Vadd(nq, tmp, 1, outarray[i + 1], 1, outarray[i + 1], 1);
255 Vmath::Vadd(nq, tmp, 1, outarray[i + 1], 1, outarray[i + 1], 1);
260 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
270 int nvariables = inarray.size();
282 for (i = 0; i < 2; ++i)
299 for (i = 0; i < nvariables - 1; ++i)
319 "Variable depth not supported for the Peregrine "
333 m_fields[0]->IProductWRTBase(outarray[1], modarray[0]);
334 m_fields[0]->IProductWRTBase(outarray[2], modarray[1]);
344 m_fields[0]->MultiplyByElmtInvMass(modarray[0], modarray[0]);
345 m_fields[0]->MultiplyByElmtInvMass(modarray[1], modarray[1]);
357 for (i = 0; i < 2; ++i)
395 m_fields[0]->IProductWRTDerivBase(0, physfield[0], coeffsfield[0]);
396 m_fields[0]->IProductWRTDerivBase(1, physfield[1], coeffsfield[1]);
397 Vmath::Vadd(ncoeffs, coeffsfield[0], 1, coeffsfield[1], 1,
406 m_fields[0]->AddTraceIntegral(normflux, coeffsfield[0]);
407 m_fields[0]->MultiplyByElmtInvMass(coeffsfield[0], coeffsfield[0]);
408 m_fields[0]->BwdTrans(coeffsfield[0], physfield[0]);
410 Vmath::Smul(nq, -invgamma, physfield[0], 1, physfield[0], 1);
438 Vmath::Smul(nq, gamma, physfield[0], 1, physfield[0], 1);
443 m_fields[0]->IProductWRTDerivBase(0, physfield[0], coeffsfield[0]);
444 m_fields[0]->IProductWRTDerivBase(1, physfield[0], coeffsfield[1]);
456 m_fields[0]->AddTraceIntegral(normflux, coeffsfield[0]);
459 m_fields[0]->AddTraceIntegral(normflux, coeffsfield[1]);
462 Vmath::Vadd(ncoeffs, f1, 1, coeffsfield[0], 1, modarray[0], 1);
463 Vmath::Vadd(ncoeffs, f2, 1, coeffsfield[1], 1, modarray[1], 1);
465 m_fields[1]->MultiplyByElmtInvMass(modarray[0], modarray[0]);
466 m_fields[2]->MultiplyByElmtInvMass(modarray[1], modarray[1]);
468 m_fields[1]->BwdTrans(modarray[0], outarray[1]);
469 m_fields[2]->BwdTrans(modarray[1], outarray[2]);
478 ASSERTL0(
false,
"Unknown projection scheme for the Peregrine");
481 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
491 int nvariables = inarray.size();
499 if (inarray != outarray)
503 for (i = 0; i < nvariables; ++i)
519 for (i = 0; i < nvariables; ++i)
521 m_fields[i]->FwdTrans(inarray[i], coeffs);
522 m_fields[i]->BwdTrans(coeffs, outarray[i]);
527 ASSERTL0(
false,
"Unknown projection scheme");
544 for (
int i = 0; i < nvariables - 1; ++i)
547 m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
551 for (
int n = 0; n <
m_fields[0]->GetBndConditions().size(); ++n)
555 if (boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
562 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
564 for (
int i = 0; i < nvariables; ++i)
566 m_fields[i]->EvaluateBoundaryConditions(time);
569 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
582 int nvariables = physarray.size();
586 int e, id1, id2, npts;
588 m_fields[0]->GetBndCondExpansions()[bcRegion];
589 for (e = 0; e < bcexp->GetExpSize(); ++e)
591 npts = bcexp->GetExp(e)->GetTotPoints();
592 id1 = bcexp->GetPhys_Offset(e);
593 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
594 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
603 &tmp[0], 1, &tmp[0], 1);
613 &Fwd[1 + i][id2], 1, &Fwd[1 + i][id2], 1);
617 for (i = 0; i < nvariables; ++i)
619 bcexp =
m_fields[i]->GetBndCondExpansions()[bcRegion];
620 Vmath::Vcopy(npts, &Fwd[i][id2], 1, &(bcexp->UpdatePhys())[id1], 1);
634 int e, id1, id2, npts;
636 m_fields[0]->GetBndCondExpansions()[bcRegion];
638 for (e = 0; e < bcexp->GetExpSize(); ++e)
640 npts = bcexp->GetExp(e)->GetNumPoints(0);
641 id1 = bcexp->GetPhys_Offset(e);
642 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
643 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
661 &tmp_n[0], 1, &tmp_n[0], 1);
666 &tmp_t[0], 1, &tmp_t[0], 1);
675 &Fwd[1][id2], 1, &Fwd[1][id2], 1);
680 &Fwd[2][id2], 1, &Fwd[2][id2], 1);
685 "3D not implemented for Shallow Water Equations");
688 ASSERTL0(
false,
"Illegal expansion dimension");
692 for (i = 0; i < nvariables; ++i)
694 bcexp =
m_fields[i]->GetBndCondExpansions()[bcRegion];
695 Vmath::Vcopy(npts, &Fwd[i][id2], 1, &(bcexp->UpdatePhys())[id1], 1);
706 int nq =
m_fields[0]->GetTotPoints();
722 Vmath::Vmul(nq, physfield[0], 1, physfield[0], 1, tmp, 1);
735 Vmath::Vadd(nq, flux[i + 1][i], 1, tmp, 1, flux[i + 1][i], 1);
745 if (physin.get() == physout.get())
749 for (
int i = 0; i < 3; ++i)
760 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
763 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
771 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
774 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
802 if (physin.get() == physout.get())
806 for (
int i = 0; i < 3; ++i)
817 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
820 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
828 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
831 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
863 const int npts = physfield[0].size();
886 for (
int j = 0; j < nq; j++)
888 (
m_fields[3]->UpdatePhys())[j] = fce[j];
915 for (i = 0; i < 2; ++i)
926 m_fields[1]->GetFwdBwdTracePhys(inarray[0], Fwd[0], Bwd[0]);
927 m_fields[2]->GetFwdBwdTracePhys(inarray[1], Fwd[1], Bwd[1]);
932 for (i = 0; i < nTraceNumPoints; ++i)
934 numfluxX[i] = 0.5 * (Fwd[0][i] + Bwd[0][i]);
935 numfluxY[i] = 0.5 * (Fwd[1][i] + Bwd[1][i]);
947 for (
int n = 0; n <
m_fields[0]->GetBndConditions().size(); ++n)
951 if (boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
958 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
960 ASSERTL0(
false,
"time-dependent BC not implemented for Boussinesq");
962 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
975 for (
int i = 0; i < nvariables; ++i)
978 m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
983 int e, id1, id2, npts;
985 m_fields[0]->GetBndCondExpansions()[bcRegion];
986 for (e = 0; e < bcexp->GetExpSize(); ++e)
988 npts = bcexp->GetExp(e)->GetTotPoints();
989 id1 = bcexp->GetPhys_Offset(e);
990 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
991 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
997 ASSERTL0(
false,
"1D not yet implemented for Boussinesq");
1008 &tmp_n[0], 1, &tmp_n[0], 1);
1013 &tmp_t[0], 1, &tmp_t[0], 1);
1022 &Fwd[0][id2], 1, &Fwd[0][id2], 1);
1027 &Fwd[1][id2], 1, &Fwd[1][id2], 1);
1031 ASSERTL0(
false,
"3D not implemented for Boussinesq equations");
1034 ASSERTL0(
false,
"Illegal expansion dimension");
1038 bcexp =
m_fields[1]->GetBndCondExpansions()[bcRegion];
1039 Vmath::Vcopy(npts, &Fwd[0][id2], 1, &(bcexp->UpdatePhys())[id1], 1);
1041 bcexp =
m_fields[2]->GetBndCondExpansions()[bcRegion];
1042 Vmath::Vcopy(npts, &Fwd[1][id2], 1, &(bcexp->UpdatePhys())[id1], 1);
1052 for (
int n = 0; n <
m_fields[0]->GetBndConditions().size(); ++n)
1056 if (boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
1062 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
1067 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize() - 1;
1078 m_fields[0]->ExtractTracePhys(inarray,
z);
1082 int e, id1, id2, npts;
1084 m_fields[0]->GetBndCondExpansions()[bcRegion];
1086 for (e = 0; e < bcexp->GetExpSize(); ++e)
1088 npts = bcexp->GetExp(e)->GetTotPoints();
1089 id1 = bcexp->GetPhys_Offset(e);
1090 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
1091 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1095 bcexp =
m_fields[1]->GetBndCondExpansions()[bcRegion];
1096 Vmath::Vcopy(npts, &
z[id2], 1, &(bcexp->UpdatePhys())[id1], 1);
1117 m_fields[1]->GetFwdBwdTracePhys(physfield, Fwd, Bwd);
1122 for (i = 0; i < nTraceNumPoints; ++i)
1124 outX[i] = 0.5 * (Fwd[i] + Bwd[i]);
1125 outY[i] = 0.5 * (Fwd[i] + Bwd[i]);
1147 for (
int i = 0; i < nq; i++)
1150 amp * pow((1.0 / cosh(
sqrt(0.75 * (amp / (
d *
d *
d))) *
1151 (
A * (x0[i] + x_offset) - C * time))),
1155 pow((1.0 / cosh(
sqrt(0.75 * (amp / (
d *
d *
d))) *
1156 (
A * (x0[i] + x_offset) - C * time))),
1168 for (
int i = 0; i < 4; ++i)
1180 NekDouble initialtime,
bool dumpInitialConditions,
1181 [[maybe_unused]]
const int domain)
1199 !
m_comm->IsParallelInTime())
1203 else if (dumpInitialConditions &&
m_nchk == 0 &&
m_comm->IsParallelInTime())
1206 if (!fs::is_directory(newdir))
1208 fs::create_directory(newdir);
1210 if (
m_comm->GetTimeComm()->GetRank() == 0)
#define ASSERTL0(condition, msg)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
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)
void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0) override
Set the initial conditions.
void v_PrimitiveToConservative() override
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void v_ConservativeToPrimitive() override
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
void SetBoundaryConditionsForcing(Array< OneD, Array< OneD, NekDouble > > &inarray, NekDouble time)
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
~NonlinearPeregrine() override
void WallBoundaryForcing(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &inarray)
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print a summary of time stepping parameters.
ProblemType m_problemType
void LaitoneSolitaryWave(NekDouble amp, NekDouble d, NekDouble time, NekDouble x_offset)
void AddVariableDepth(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
static std::string className
Name of class.
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
void WallBoundary(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary condition.
void NumericalFluxForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &numfluxX, Array< OneD, NekDouble > &numfluxY)
void SetBoundaryConditionsContVariables(Array< OneD, NekDouble > &inarray, NekDouble time)
void AddCoriolis(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
void WallBoundaryContVariables(int bcRegion, int cnt, Array< OneD, NekDouble > &inarray)
void NumericalFluxConsVariables(Array< OneD, NekDouble > &physfield, Array< OneD, NekDouble > &outX, Array< OneD, NekDouble > &outY)
StdRegions::ConstFactorMap m_factors
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
void v_InitObject(bool DeclareFields=true) override
Init object for UnsteadySystem class.
NonlinearPeregrine(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void WCESolve(Array< OneD, NekDouble > &fce, NekDouble lambda)
Base class for unsteady solvers.
NekDouble m_g
Acceleration of gravity.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
Array< OneD, Array< OneD, NekDouble > > m_bottomSlope
const Array< OneD, NekDouble > & GetDepth()
SolverUtils::AdvectionSharedPtr m_advection
void PrimitiveToConservative()
bool m_constantDepth
Indicates if constant depth case.
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print a summary of time stepping parameters.
Array< OneD, NekDouble > m_coriolis
Coriolis force.
Array< OneD, NekDouble > m_depth
Still water depth.
void v_InitObject(bool DeclareFields=true) override
Init object for UnsteadySystem class.
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
void ConservativeToPrimitive()
int m_spacedim
Spatial dimension (>= expansion dim).
int m_expdim
Expansion dimension.
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
LibUtilities::CommSharedPtr m_comm
Communicator.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
std::string m_sessionName
Name of the session.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
SOLVER_UTILS_EXPORT int GetNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
int m_nchk
Number of checkpoints written so far.
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT int GetTotPoints()
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
const std::vector< NekDouble > velocity
@ eMixed_CG_Discontinuous
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
std::vector< std::pair< std::string, std::string > > SummaryList
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::vector< double > z(NPUPPER)
std::vector< double > d(NPUPPER *NPUPPER)
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
const char *const ProblemTypeMap[]
@ SIZE_ProblemType
Length of enum list.
@ eSolitaryWave
First order Laitone solitary wave.
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.
void Neg(int n, T *x, const int incx)
Negate x = -x.
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
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.
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
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
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.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
scalarT< T > sqrt(scalarT< T > in)