38 #include <boost/algorithm/string.hpp>
47 string LinearSWE::className =
49 "LinearSWE", LinearSWE::create,
50 "Linear shallow water equation in primitive variables.");
69 ASSERTL0(
false,
"Implicit SWE not set up.");
92 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
104 m_session->LoadSolverInfo(
"UpwindType", riemName,
"NoSolver");
107 ASSERTL0(
false,
"LinearHLL only valid for constant depth");
128 int nTracePointsTot =
m_fields[0]->GetTrace()->GetTotPoints();
160 ASSERTL0(
false,
"Unsupported projection type.");
190 m_fields[0]->IProductWRTBase(tmp,mod);
191 m_fields[0]->MultiplyByElmtInvMass(mod,mod);
198 m_fields[0]->IProductWRTBase(tmp,mod);
199 m_fields[0]->MultiplyByElmtInvMass(mod,mod);
218 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
231 int nvariables = inarray.num_elements();
253 for(i = 0; i < nvariables; ++i)
280 fluxvector(nvariables);
282 for (i = 0; i < nvariables; ++i)
285 for(j = 0; j < ndim; ++j)
301 for(i = 0; i < nvariables; ++i)
322 ASSERTL0(
false,
"Unknown projection scheme for the NonlinearSWE");
333 int nvariables = inarray.num_elements();
344 for(i = 0; i < nvariables; ++i)
358 for(i = 0; i < nvariables; ++i)
360 m_fields[i]->FwdTrans(inarray[i],coeffs);
361 m_fields[i]->BwdTrans_IterPerExp(coeffs,outarray[i]);
366 ASSERTL0(
false,
"Unknown projection scheme");
378 int nvariables =
m_fields.num_elements();
385 for (
int i = 0; i < nvariables; ++i)
388 m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
392 for(
int n = 0; n <
m_fields[0]->GetBndConditions().num_elements(); ++n)
395 if (boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
"Wall"))
401 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
403 for (
int i = 0; i < nvariables; ++i)
406 m_fields[i]->EvaluateBoundaryConditions(time, varName);
409 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
424 int nvariables = physarray.num_elements();
428 int e, id1, id2,
npts;
430 for (e = 0; e <
m_fields[0]->GetBndCondExpansions()[bcRegion]
433 npts =
m_fields[0]->GetBndCondExpansions()[bcRegion]->
434 GetExp(e)->GetTotPoints();
435 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->
437 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
439 GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
468 for (i = 0; i < nvariables; ++i)
471 &(
m_fields[i]->GetBndCondExpansions()[bcRegion]->
472 UpdatePhys())[id1], 1);
482 int nvariables = physarray.num_elements();
486 int e, id1, id2,
npts;
488 for(e = 0; e <
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize(); ++e)
490 npts =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExp(e)->GetNumPoints(0);
491 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e) ;
492 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
m_fields[0]->GetTraceMap()->GetBndCondCoeffsToGlobalCoeffsMap(cnt+e));
525 ASSERTL0(
false,
"3D not implemented for Shallow Water Equations");
528 ASSERTL0(
false,
"Illegal expansion dimension");
534 for (i = 0; i < nvariables; ++i)
536 Vmath::Vcopy(npts,&Fwd[i][id2], 1,&(
m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1],1);
548 int nq =
m_fields[0]->GetTotPoints();
572 Vmath::Vadd(nq, flux[i+1][i], 1, tmp, 1, flux[i+1][i], 1);
582 if(physin.get() == physout.get())
586 for (
int i = 0; i < 3; ++i)
608 Vmath::Vdiv(nq,physin[1],1,physin[0],1,physout[1],1);
611 Vmath::Vdiv(nq,physin[2],1,physin[0],1,physout[2],1);
636 if(physin.get() == physout.get())
640 for (
int i = 0; i < 3; ++i)
651 Vmath::Vmul(nq,physout[0],1,tmp[1],1,physout[1],1);
654 Vmath::Vmul(nq,physout[0],1,tmp[2],1,physout[2],1);
663 Vmath::Vmul(nq,physout[0],1,physin[1],1,physout[1],1);
666 Vmath::Vmul(nq,physout[0],1,physin[2],1,physout[2],1);
698 const int npts = physfield[0].num_elements();
710 if (
m_session->DefinesSolverInfo(
"UpwindType"))
712 std::string UpwindType;
713 UpwindType =
m_session->GetSolverInfo(
"UpwindType");
714 if (UpwindType ==
"LinearAverage")
718 if (UpwindType ==
"LinearHLL")
Array< OneD, NekDouble > m_coriolis
Coriolis force.
const Array< OneD, NekDouble > & GetDepthFwd()
virtual void v_InitObject()
Init object for UnsteadySystem class.
#define ASSERTL0(condition, msg)
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
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.
Array< OneD, NekDouble > m_dFwd
Still water depth traces.
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
virtual void v_PrimitiveToConservative()
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
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()
void AddCoriolis(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
SolverUtils::AdvectionSharedPtr m_advection
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.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
SOLVER_UTILS_EXPORT int GetTotPoints()
virtual void v_ConservativeToPrimitive()
void ConservativeToPrimitive()
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
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 WallBoundary(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary condition.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
virtual void v_InitObject()
Init object for UnsteadySystem class.
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).
void CopyBoundaryTrace(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
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.
Array< OneD, NekDouble > m_dBwd
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
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[]
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
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
SOLVER_UTILS_EXPORT int GetNcoeffs()
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
bool m_constantDepth
Indicates if constant depth case.
const Array< OneD, NekDouble > & GetDepthBwd()
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
void Zero(int n, T *x, const int incx)
Zero vector.
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.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.