39 #include <boost/algorithm/string.hpp>
49 string NonlinearPeregrine::className =
51 "NonlinearPeregrine", NonlinearPeregrine::create,
52 "Nonlinear Peregrine equations in conservative variables.");
54 NonlinearPeregrine::NonlinearPeregrine(
68 if (
m_session->DefinesSolverInfo(
"PROBLEMTYPE"))
71 std::string ProblemTypeStr =
m_session->GetSolverInfo(
"PROBLEMTYPE");
93 ASSERTL0(
false,
"Implicit Peregrine not set up.");
98 if (
m_session->DefinesParameter(
"ConstDepth"))
104 ASSERTL0(
false,
"Constant Depth not specified");
114 "Continuous projection type not supported for Peregrine.");
128 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
137 m_session->LoadSolverInfo(
"UpwindType", riemName,
"NoSolver");
160 ASSERTL0(
false,
"Unsupported projection type.");
190 m_fields[0]->IProductWRTBase(tmp, mod);
191 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
193 Vmath::Vadd(nq, tmp, 1, outarray[1], 1, outarray[1], 1);
198 m_fields[0]->IProductWRTBase(tmp, mod);
199 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
201 Vmath::Vadd(nq, tmp, 1, outarray[2], 1, outarray[2], 1);
209 Vmath::Vadd(nq, tmp, 1, outarray[1], 1, outarray[1], 1);
214 Vmath::Vadd(nq, tmp, 1, outarray[2], 1, outarray[2], 1);
218 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
244 m_fields[0]->IProductWRTBase(tmp, mod);
245 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
247 Vmath::Vadd(nq, tmp, 1, outarray[i + 1], 1, outarray[i + 1], 1);
258 Vmath::Vadd(nq, tmp, 1, outarray[i + 1], 1, outarray[i + 1], 1);
263 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
274 int nvariables = inarray.num_elements();
287 for (i = 0; i < nvariables; ++i)
305 for (i = 0; i < nvariables - 1; ++i)
325 "Variable depth not supported for the Peregrine "
339 m_fields[0]->IProductWRTBase(outarray[1], modarray[1]);
340 m_fields[0]->IProductWRTBase(outarray[2], modarray[2]);
350 m_fields[0]->MultiplyByElmtInvMass(modarray[1], modarray[1]);
351 m_fields[0]->MultiplyByElmtInvMass(modarray[2], modarray[2]);
363 for (i = 0; i < 2; ++i)
403 m_fields[0]->IProductWRTDerivBase(0, physfield[0], coeffsfield[0]);
404 m_fields[0]->IProductWRTDerivBase(1, physfield[1], coeffsfield[1]);
405 Vmath::Vadd(ncoeffs, coeffsfield[0], 1, coeffsfield[1], 1,
412 m_fields[0]->AddTraceIntegral(upwindX[0], upwindY[0],
414 m_fields[0]->MultiplyByElmtInvMass(coeffsfield[0], coeffsfield[0]);
415 m_fields[0]->BwdTrans(coeffsfield[0], physfield[0]);
417 Vmath::Smul(nq, -invgamma, physfield[0], 1, physfield[0], 1);
445 Vmath::Smul(nq, gamma, physfield[0], 1, physfield[0], 1);
450 m_fields[0]->IProductWRTDerivBase(0, physfield[0], coeffsfield[0]);
451 m_fields[1]->IProductWRTDerivBase(1, physfield[0], coeffsfield[1]);
462 m_fields[0]->AddTraceIntegral(upwindX[0], uptemp,
464 m_fields[0]->AddTraceIntegral(uptemp, upwindY[0],
468 Vmath::Vadd(ncoeffs, f1, 1, coeffsfield[0], 1, modarray[1], 1);
469 Vmath::Vadd(ncoeffs, f2, 1, coeffsfield[1], 1, modarray[2], 1);
471 m_fields[1]->MultiplyByElmtInvMass(modarray[1], modarray[1]);
472 m_fields[2]->MultiplyByElmtInvMass(modarray[2], modarray[2]);
474 m_fields[1]->BwdTrans(modarray[1], outarray[1]);
475 m_fields[2]->BwdTrans(modarray[2], outarray[2]);
484 ASSERTL0(
false,
"Unknown projection scheme for the Peregrine");
487 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
498 int nvariables = inarray.num_elements();
508 for (i = 0; i < nvariables; ++i)
523 for (i = 0; i < nvariables; ++i)
525 m_fields[i]->FwdTrans(inarray[i], coeffs);
526 m_fields[i]->BwdTrans_IterPerExp(coeffs, outarray[i]);
531 ASSERTL0(
false,
"Unknown projection scheme");
542 int nvariables =
m_fields.num_elements();
549 for (
int i = 0; i < nvariables; ++i)
552 m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
556 for (
int n = 0; n <
m_fields[0]->GetBndConditions().num_elements(); ++n)
560 if (boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
"Wall"))
566 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
568 for (
int i = 0; i < nvariables; ++i)
570 m_fields[i]->EvaluateBoundaryConditions(time);
573 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
586 int nvariables = physarray.num_elements();
590 int e, id1, id2,
npts;
592 m_fields[0]->GetBndCondExpansions()[bcRegion];
593 for (e = 0; e < bcexp->GetExpSize(); ++e)
595 npts = bcexp->GetExp(e)->GetTotPoints();
596 id1 = bcexp->GetPhys_Offset(e);
597 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
598 m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(
608 &tmp[0], 1, &tmp[0], 1);
618 &Fwd[1 + i][id2], 1, &Fwd[1 + i][id2], 1);
622 for (i = 0; i < nvariables; ++i)
624 bcexp =
m_fields[i]->GetBndCondExpansions()[bcRegion];
625 Vmath::Vcopy(npts, &Fwd[i][id2], 1, &(bcexp->UpdatePhys())[id1], 1);
642 int e, id1, id2,
npts;
644 m_fields[0]->GetBndCondExpansions()[bcRegion];
646 for (e = 0; e < bcexp->GetExpSize();
649 npts = bcexp->GetExp(e)->GetNumPoints(0);
650 id1 = bcexp->GetPhys_Offset(e);
651 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
652 m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(
671 &tmp_n[0], 1, &tmp_n[0], 1);
676 &tmp_t[0], 1, &tmp_t[0], 1);
685 &Fwd[1][id2], 1, &Fwd[1][id2], 1);
690 &Fwd[2][id2], 1, &Fwd[2][id2], 1);
695 "3D not implemented for Shallow Water Equations");
698 ASSERTL0(
false,
"Illegal expansion dimension");
702 for (i = 0; i < nvariables; ++i)
704 bcexp =
m_fields[i]->GetBndCondExpansions()[bcRegion];
705 Vmath::Vcopy(npts, &Fwd[i][id2], 1, &(bcexp->UpdatePhys())[id1], 1);
716 int nq =
m_fields[0]->GetTotPoints();
732 Vmath::Vmul(nq, physfield[0], 1, physfield[0], 1, tmp, 1);
740 Vmath::Vmul(nq, velocity[j], 1, physfield[i + 1], 1,
745 Vmath::Vadd(nq, flux[i + 1][i], 1, tmp, 1, flux[i + 1][i], 1);
756 if (physin.get() == physout.get())
760 for (
int i = 0; i < 3; ++i)
771 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
774 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
782 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
785 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
813 if (physin.get() == physout.get())
817 for (
int i = 0; i < 3; ++i)
828 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
831 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
840 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
843 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
877 const int npts = physfield[0].num_elements();
881 Vmath::Vdiv(npts, physfield[1 + i], 1, physfield[0], 1, velocity[i], 1);
902 for (
int j = 0; j < nq; j++)
904 (
m_fields[3]->UpdatePhys())[j] = fce[j];
934 for (i = 0; i < 2; ++i)
945 m_fields[1]->GetFwdBwdTracePhys(inarray[0], Fwd[0], Bwd[0]);
946 m_fields[2]->GetFwdBwdTracePhys(inarray[1], Fwd[1], Bwd[1]);
951 for (i = 0; i < nTraceNumPoints; ++i)
953 numfluxX[i] = 0.5 * (Fwd[0][i] + Bwd[0][i]);
954 numfluxY[i] = 0.5 * (Fwd[1][i] + Bwd[1][i]);
966 for (
int n = 0; n <
m_fields[0]->GetBndConditions().num_elements(); ++n)
970 if (boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
"Wall"))
976 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
978 ASSERTL0(
false,
"time-dependent BC not implemented for Boussinesq");
980 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
998 for (
int i = 0; i < nvariables; ++i)
1001 m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
1006 int e, id1, id2,
npts;
1008 m_fields[0]->GetBndCondExpansions()[bcRegion];
1009 for (e = 0; e < bcexp->GetExpSize(); ++e)
1011 npts = bcexp->GetExp(e)->GetTotPoints();
1012 id1 = bcexp->GetPhys_Offset(e);
1013 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
1014 m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(
1021 ASSERTL0(
false,
"1D not yet implemented for Boussinesq");
1032 &tmp_n[0], 1, &tmp_n[0], 1);
1037 &tmp_t[0], 1, &tmp_t[0], 1);
1046 &Fwd[0][id2], 1, &Fwd[0][id2], 1);
1051 &Fwd[1][id2], 1, &Fwd[1][id2], 1);
1055 ASSERTL0(
false,
"3D not implemented for Boussinesq equations");
1058 ASSERTL0(
false,
"Illegal expansion dimension");
1062 bcexp =
m_fields[1]->GetBndCondExpansions()[bcRegion];
1063 Vmath::Vcopy(npts, &Fwd[0][id2], 1, &(bcexp->UpdatePhys())[id1], 1);
1065 bcexp =
m_fields[2]->GetBndCondExpansions()[bcRegion];
1066 Vmath::Vcopy(npts, &Fwd[1][id2], 1, &(bcexp->UpdatePhys())[id1], 1);
1077 for (
int n = 0; n <
m_fields[0]->GetBndConditions().num_elements(); ++n)
1081 if(boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
"Wall"))
1086 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
1091 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize() - 1;
1104 m_fields[0]->ExtractTracePhys(inarray, z);
1108 int e, id1, id2,
npts;
1110 m_fields[0]->GetBndCondExpansions()[bcRegion];
1112 for (e = 0; e < bcexp->GetExpSize(); ++e)
1114 npts = bcexp->GetExp(e)->GetTotPoints();
1115 id1 = bcexp->GetPhys_Offset(e);
1116 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
1117 m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(
1121 bcexp =
m_fields[1]->GetBndCondExpansions()[bcRegion];
1122 Vmath::Vcopy(npts, &z[id2], 1, &(bcexp->UpdatePhys())[id1], 1);
1147 m_fields[1]->GetFwdBwdTracePhys(physfield, Fwd[0], Bwd[0]);
1152 for (i = 0; i < nTraceNumPoints; ++i)
1154 outX[i] = 0.5 * (Fwd[0][i] + Bwd[0][i]);
1155 outY[i] = 0.5 * (Fwd[0][i] + Bwd[0][i]);
1180 for (
int i = 0; i < nq; i++)
1182 (
m_fields[0]->UpdatePhys())[i] = amp * pow((1.0 / cosh(
1183 sqrt(0.75 * (amp / (d * d * d))) *
1184 (A * (x0[i] + x_offset) - C * time))), 2.0);
1185 (
m_fields[1]->UpdatePhys())[i] = (amp / d) * pow((1.0 / cosh(
1186 sqrt(0.75 * (amp / (d * d * d))) *
1187 (A * (x0[i] + x_offset) - C * time)
1188 )), 2.0) * sqrt(
m_g * d);
1198 for (
int i = 0; i < 4; ++i)
1212 bool dumpInitialConditions,
1230 if (dumpInitialConditions)
Array< OneD, NekDouble > m_coriolis
Coriolis force.
void SetBoundaryConditionsContVariables(Array< OneD, NekDouble > &inarray, NekDouble time)
#define ASSERTL0(condition, msg)
void NumericalFluxForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &numfluxX, Array< OneD, NekDouble > &numfluxY)
void PrimitiveToConservative()
Base class for unsteady solvers.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Array< OneD, NekDouble > m_depth
Still water depth.
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
void WallBoundaryForcing(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &inarray)
std::vector< std::pair< std::string, std::string > > SummaryList
int m_expdim
Expansion dimension.
ProblemType m_problemType
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
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
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
SolverUtils::AdvectionSharedPtr m_advection
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
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 SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
const Array< OneD, NekDouble > & GetDepth()
SOLVER_UTILS_EXPORT int GetTotPoints()
StdRegions::ConstFactorMap m_factors
virtual void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Set the initial conditions.
void ConservativeToPrimitive()
void NumericalFluxConsVariables(Array< OneD, NekDouble > &physfield, Array< OneD, NekDouble > &outX, Array< OneD, NekDouble > &outY)
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
Array< OneD, Array< OneD, NekDouble > > m_bottomSlope
virtual ~NonlinearPeregrine()
problem type selector
virtual void v_InitObject()
Init object for UnsteadySystem class.
void WallBoundary(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary condition.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
virtual void v_ConservativeToPrimitive()
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
virtual void v_InitObject()
Init object for UnsteadySystem class.
const char *const ProblemTypeMap[]
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
RiemannSolverFactory & GetRiemannSolverFactory()
void WCESolve(Array< OneD, NekDouble > &fce, NekDouble lambda)
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
int m_spacedim
Spatial dimension (>= expansion dim).
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
void Neg(int n, T *x, const int incx)
Negate x = -x.
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
void SetBoundaryConditionsForcing(Array< OneD, Array< OneD, NekDouble > > &inarray, NekDouble time)
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
EquationSystemFactory & GetEquationSystemFactory()
void AddVariableDepth(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
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.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector plus vector): z = w*x - y
virtual void v_PrimitiveToConservative()
SOLVER_UTILS_EXPORT int GetNcoeffs()
First order Laitone solitary wave.
bool m_constantDepth
Indicates if constant depth case.
void WallBoundaryContVariables(int bcRegion, int cnt, Array< OneD, NekDouble > &inarray)
void AddCoriolis(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
NekDouble m_g
Acceleration of gravity.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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 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.
static FlagList NullFlagList
An empty flag list.
void LaitoneSolitaryWave(NekDouble amp, NekDouble d, NekDouble time, NekDouble x_offset)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.