35#include <boost/algorithm/string.hpp>
49 "Linear shallow water equation in primitive variables.");
68 ASSERTL0(
false,
"Implicit SWE not set up.");
89 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
95 m_session->LoadSolverInfo(
"UpwindType", riemName,
"NoSolver");
98 ASSERTL0(
false,
"LinearHLL only valid for constant depth");
110 int nTracePointsTot =
m_fields[0]->GetTrace()->GetTotPoints();
127 ASSERTL0(
false,
"Unsupported projection type.");
155 m_fields[0]->IProductWRTBase(tmp, mod);
156 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
158 Vmath::Vadd(nq, tmp, 1, outarray[1], 1, outarray[1], 1);
163 m_fields[0]->IProductWRTBase(tmp, mod);
164 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
166 Vmath::Vadd(nq, tmp, 1, outarray[2], 1, outarray[2], 1);
174 Vmath::Vadd(nq, tmp, 1, outarray[1], 1, outarray[1], 1);
179 Vmath::Vadd(nq, tmp, 1, outarray[2], 1, outarray[2], 1);
183 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
194 int nvariables = inarray.size();
214 for (i = 0; i < nvariables; ++i)
241 for (i = 0; i < nvariables; ++i)
244 for (j = 0; j < ndim; ++j)
259 for (i = 0; i < nvariables; ++i)
262 fluxvector[i][0], tmp);
264 fluxvector[i][1], tmp1);
282 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
292 int nvariables = inarray.size();
300 if (inarray != outarray)
304 for (i = 0; i < nvariables; ++i)
319 for (i = 0; i < nvariables; ++i)
321 m_fields[i]->FwdTrans(inarray[i], coeffs);
322 m_fields[i]->BwdTrans(coeffs, outarray[i]);
327 ASSERTL0(
false,
"Unknown projection scheme");
344 for (
int i = 0; i < nvariables; ++i)
347 m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
351 for (
int n = 0; n <
m_fields[0]->GetBndConditions().size(); ++n)
353 if (
m_fields[0]->GetBndConditions()[n]->GetBoundaryConditionType() ==
360 if (boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
367 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
369 for (
int i = 0; i < nvariables; ++i)
372 m_fields[i]->EvaluateBoundaryConditions(time, varName);
375 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
388 int nvariables = physarray.size();
392 int e, id1, id2, npts;
398 ->GetBndCondExpansions()[bcRegion]
401 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
402 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
403 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
412 &tmp[0], 1, &tmp[0], 1);
422 &Fwd[1 + i][id2], 1, &Fwd[1 + i][id2], 1);
426 for (i = 0; i < nvariables; ++i)
430 ->GetBndCondExpansions()[bcRegion]
431 ->UpdatePhys())[id1],
443 int nvariables = physarray.size();
447 int e, id1, id2, npts;
453 ->GetBndCondExpansions()[bcRegion]
456 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
457 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
458 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
476 &tmp_n[0], 1, &tmp_n[0], 1);
481 &tmp_t[0], 1, &tmp_t[0], 1);
490 &Fwd[1][id2], 1, &Fwd[1][id2], 1);
495 &Fwd[2][id2], 1, &Fwd[2][id2], 1);
500 "3D not implemented for Shallow Water Equations");
503 ASSERTL0(
false,
"Illegal expansion dimension");
507 for (i = 0; i < nvariables; ++i)
511 ->GetBndCondExpansions()[bcRegion]
512 ->UpdatePhys())[id1],
524 int nq =
m_fields[0]->GetTotPoints();
549 Vmath::Vadd(nq, flux[i + 1][i], 1, tmp, 1, flux[i + 1][i], 1);
559 if (physin.get() == physout.get())
563 for (
int i = 0; i < 3; ++i)
574 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
577 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
585 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
588 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
616 if (physin.get() == physout.get())
620 for (
int i = 0; i < 3; ++i)
631 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
634 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
642 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
645 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
677 const int npts = physfield[0].size();
688 if (
m_session->DefinesSolverInfo(
"UpwindType"))
#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_GenerateSummary(SolverUtils::SummaryList &s) override
Print a summary of time stepping parameters.
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
LinearSWE(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
const Array< OneD, NekDouble > & GetDepthBwd()
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
const Array< OneD, NekDouble > & GetDepthFwd()
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
void v_ConservativeToPrimitive() override
void v_InitObject(bool DeclareFields=true) override
Init object for UnsteadySystem class.
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
static std::string className
Name of class.
void v_PrimitiveToConservative() override
Array< OneD, NekDouble > m_dFwd
Still water depth traces.
Array< OneD, NekDouble > m_dBwd
void WallBoundary(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary condition.
void AddCoriolis(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Base class for unsteady solvers.
NekDouble m_g
Acceleration of gravity.
void CopyBoundaryTrace(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
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 Zero(int n, T *x, const int incx)
Zero 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.