35#include <boost/algorithm/string.hpp>
49 "Nonlinear shallow water equation in conservative variables.");
68 ASSERTL0(
false,
"Implicit SWE not set up.");
89 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
95 m_session->LoadSolverInfo(
"UpwindType", riemName,
"Average");
115 ASSERTL0(
false,
"Unsupported projection type.");
143 m_fields[0]->IProductWRTBase(tmp, mod);
144 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
146 Vmath::Vadd(nq, tmp, 1, outarray[1], 1, outarray[1], 1);
151 m_fields[0]->IProductWRTBase(tmp, mod);
152 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
154 Vmath::Vadd(nq, tmp, 1, outarray[2], 1, outarray[2], 1);
162 Vmath::Vadd(nq, tmp, 1, outarray[1], 1, outarray[1], 1);
167 Vmath::Vadd(nq, tmp, 1, outarray[2], 1, outarray[2], 1);
171 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
196 m_fields[0]->IProductWRTBase(tmp, mod);
197 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
199 Vmath::Vadd(nq, tmp, 1, outarray[i + 1], 1, outarray[i + 1], 1);
210 Vmath::Vadd(nq, tmp, 1, outarray[i + 1], 1, outarray[i + 1], 1);
215 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
226 int nvariables = inarray.size();
246 for (i = 0; i < nvariables; ++i)
279 for (i = 0; i < nvariables; ++i)
282 for (j = 0; j < ndim; ++j)
297 for (i = 0; i < nvariables; ++i)
300 fluxvector[i][0], tmp);
302 fluxvector[i][1], tmp1);
326 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
336 int nvariables = inarray.size();
344 if (inarray != outarray)
348 for (i = 0; i < nvariables; ++i)
363 for (i = 0; i < nvariables; ++i)
365 m_fields[i]->FwdTrans(inarray[i], coeffs);
366 m_fields[i]->BwdTrans(coeffs, outarray[i]);
371 ASSERTL0(
false,
"Unknown projection scheme");
388 for (
int i = 0; i < nvariables; ++i)
391 m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
395 for (
int n = 0; n <
m_fields[0]->GetBndConditions().size(); ++n)
398 if (
m_fields[0]->GetBndConditions()[n]->GetBoundaryConditionType() ==
405 if (boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
412 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
414 for (
int i = 0; i < nvariables; ++i)
417 m_fields[i]->EvaluateBoundaryConditions(time, varName);
420 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
433 int nvariables = physarray.size();
437 int e, id1, id2, npts;
443 ->GetBndCondExpansions()[bcRegion]
446 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
447 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
448 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
457 &tmp[0], 1, &tmp[0], 1);
467 &Fwd[1 + i][id2], 1, &Fwd[1 + i][id2], 1);
471 for (i = 0; i < nvariables; ++i)
475 ->GetBndCondExpansions()[bcRegion]
476 ->UpdatePhys())[id1],
488 int nvariables = physarray.size();
492 int e, id1, id2, npts;
498 ->GetBndCondExpansions()[bcRegion]
501 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
502 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
503 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
521 &tmp_n[0], 1, &tmp_n[0], 1);
526 &tmp_t[0], 1, &tmp_t[0], 1);
535 &Fwd[1][id2], 1, &Fwd[1][id2], 1);
540 &Fwd[2][id2], 1, &Fwd[2][id2], 1);
545 "3D not implemented for Shallow Water Equations");
548 ASSERTL0(
false,
"Illegal expansion dimension");
552 for (i = 0; i < nvariables; ++i)
556 ->GetBndCondExpansions()[bcRegion]
557 ->UpdatePhys())[id1],
569 int nq =
m_fields[0]->GetTotPoints();
585 Vmath::Vmul(nq, physfield[0], 1, physfield[0], 1, tmp, 1);
598 Vmath::Vadd(nq, flux[i + 1][i], 1, tmp, 1, flux[i + 1][i], 1);
608 if (physin.get() == physout.get())
612 for (
int i = 0; i < 3; ++i)
623 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
626 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
634 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
637 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
664 if (physin.get() == physout.get())
668 for (
int i = 0; i < 3; ++i)
679 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
682 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
690 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
693 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
725 const int npts = physfield[0].size();
#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 v_PrimitiveToConservative() override
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
static std::string className
Name of class.
void v_ConservativeToPrimitive() override
void AddCoriolis(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void v_GenerateSummary(SolverUtils::SummaryList &s) override
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 AddVariableDepth(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
void WallBoundary(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary condition.
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
NonlinearSWE(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
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.
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
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.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetExpSize()
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()
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.
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
MultiRegions::Direction const DirCartesianMap[]
@ eMixed_CG_Discontinuous
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
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 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.