38 #include <boost/algorithm/string.hpp>
47 string NonlinearSWE::className =
49 "NonlinearSWE", NonlinearSWE::create,
50 "Nonlinear shallow water equation in conservative variables.");
52 NonlinearSWE::NonlinearSWE(
69 ASSERTL0(
false,
"Implicit SWE not set up.");
92 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
104 m_session->LoadSolverInfo(
"UpwindType", riemName,
"Average");
146 ASSERTL0(
false,
"Unsupported projection type.");
176 m_fields[0]->IProductWRTBase(tmp,mod);
177 m_fields[0]->MultiplyByElmtInvMass(mod,mod);
184 m_fields[0]->IProductWRTBase(tmp,mod);
185 m_fields[0]->MultiplyByElmtInvMass(mod,mod);
204 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
231 m_fields[0]->IProductWRTBase(tmp,mod);
232 m_fields[0]->MultiplyByElmtInvMass(mod,mod);
234 Vmath::Vadd(nq,tmp,1,outarray[i+1],1,outarray[i+1],1);
245 Vmath::Vadd(nq,tmp,1,outarray[i+1],1,outarray[i+1],1);
250 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
263 int nvariables = inarray.num_elements();
285 for(i = 0; i < nvariables; ++i)
318 fluxvector(nvariables);
320 for (i = 0; i < nvariables; ++i)
323 for(j = 0; j < ndim; ++j)
340 for(i = 0; i < nvariables; ++i)
368 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
379 int nvariables = inarray.num_elements();
390 for(i = 0; i < nvariables; ++i)
404 for(i = 0; i < nvariables; ++i)
406 m_fields[i]->FwdTrans(inarray[i],coeffs);
407 m_fields[i]->BwdTrans_IterPerExp(coeffs,outarray[i]);
412 ASSERTL0(
false,
"Unknown projection scheme");
424 int nvariables =
m_fields.num_elements();
431 for (
int i = 0; i < nvariables; ++i)
434 m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
438 for (
int n = 0; n <
m_fields[0]->GetBndConditions().num_elements(); ++n)
442 if (boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
"Wall"))
448 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
450 for (
int i = 0; i < nvariables; ++i)
453 m_fields[i]->EvaluateBoundaryConditions(time, varName);
456 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
471 int nvariables = physarray.num_elements();
475 int e, id1, id2,
npts;
477 for (e = 0; e <
m_fields[0]->GetBndCondExpansions()[bcRegion]
480 npts =
m_fields[0]->GetBndCondExpansions()[bcRegion]->
481 GetExp(e)->GetTotPoints();
482 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->
484 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
486 GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
515 for (i = 0; i < nvariables; ++i)
518 &(
m_fields[i]->GetBndCondExpansions()[bcRegion]->
519 UpdatePhys())[id1], 1);
529 int nvariables = physarray.num_elements();
533 int e, id1, id2,
npts;
535 for(e = 0; e <
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize(); ++e)
537 npts =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExp(e)->GetNumPoints(0);
538 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e) ;
539 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
572 ASSERTL0(
false,
"3D not implemented for Shallow Water Equations");
575 ASSERTL0(
false,
"Illegal expansion dimension");
581 for (i = 0; i < nvariables; ++i)
583 Vmath::Vcopy(npts,&Fwd[i][id2], 1,&(
m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1],1);
595 int nq =
m_fields[0]->GetTotPoints();
611 Vmath::Vmul(nq, physfield[0], 1, physfield[0], 1, tmp, 1);
624 Vmath::Vadd(nq, flux[i+1][i], 1, tmp, 1, flux[i+1][i], 1);
634 if(physin.get() == physout.get())
638 for (
int i = 0; i < 3; ++i)
660 Vmath::Vdiv(nq,physin[1],1,physin[0],1,physout[1],1);
663 Vmath::Vdiv(nq,physin[2],1,physin[0],1,physout[2],1);
688 if(physin.get() == physout.get())
692 for (
int i = 0; i < 3; ++i)
703 Vmath::Vmul(nq,physout[0],1,tmp[1],1,physout[1],1);
706 Vmath::Vmul(nq,physout[0],1,tmp[2],1,physout[2],1);
715 Vmath::Vmul(nq,physout[0],1,physin[1],1,physout[1],1);
718 Vmath::Vmul(nq,physout[0],1,physin[2],1,physout[2],1);
750 const int npts = physfield[0].num_elements();
754 Vmath::Vdiv(npts, physfield[1+i], 1, physfield[0], 1,
Array< OneD, NekDouble > m_coriolis
Coriolis force.
#define ASSERTL0(condition, msg)
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
void PrimitiveToConservative()
virtual void v_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 DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
int m_expdim
Expansion dimension.
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
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
void WallBoundary(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary condition.
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 WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void AddVariableDepth(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
const Array< OneD, NekDouble > & GetDepth()
virtual void v_InitObject()
Init object for UnsteadySystem class.
SOLVER_UTILS_EXPORT int GetTotPoints()
void ConservativeToPrimitive()
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
Array< OneD, Array< OneD, NekDouble > > m_bottomSlope
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.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
virtual void v_InitObject()
Init object for UnsteadySystem class.
void AddCoriolis(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
RiemannSolverFactory & GetRiemannSolverFactory()
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
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 Neg(int n, T *x, const int incx)
Negate x = -x.
EquationSystemFactory & GetEquationSystemFactory()
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()
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
MultiRegions::Direction const DirCartesianMap[]
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
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_ConservativeToPrimitive()
SOLVER_UTILS_EXPORT int GetNcoeffs()
bool m_constantDepth
Indicates if constant depth case.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
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.
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.