Nektar++
NonlinearSWE.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: NonlinearSWE.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Nonlinear Shallow water equations in conservative variables
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36#include <iostream>
37
40
41namespace Nektar
42{
43
44std::string NonlinearSWE::className =
46 "NonlinearSWE", NonlinearSWE::create,
47 "Nonlinear shallow water equation in conservative variables.");
48
51 : ShallowWaterSystem(pSession, pGraph)
52{
53}
54
55void NonlinearSWE::v_InitObject(bool DeclareFields)
56{
58
59 // Type of advection class to be used
60 switch (m_projectionType)
61 {
62 // Continuous field
64 {
65 // Do nothing
66 break;
67 }
68 // Discontinuous field
70 {
71 std::string advName;
72 std::string diffName;
73 std::string riemName;
74
75 //---------------------------------------------------------------
76 // Setting up advection and diffusion operators
77 m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
79 advName, advName);
80 m_advection->SetFluxVector(&NonlinearSWE::GetFluxVector, this);
81
82 // Setting up Riemann solver for advection operator
83 m_session->LoadSolverInfo("UpwindType", riemName, "Average");
86 riemName, m_session);
87
88 // Setting up parameters for advection operator Riemann solver
89 m_riemannSolver->SetParam("gravity", &NonlinearSWE::GetGravity,
90 this);
91 m_riemannSolver->SetAuxVec("vecLocs", &NonlinearSWE::GetVecLocs,
92 this);
93 m_riemannSolver->SetVector("N", &NonlinearSWE::GetNormals, this);
94 m_riemannSolver->SetScalar("depth", &NonlinearSWE::GetDepth, this);
95
96 // Concluding initialisation of advection operators
97 m_advection->SetRiemannSolver(m_riemannSolver);
98 m_advection->InitObject(m_session, m_fields);
99 break;
100 }
101 default:
102 {
103 ASSERTL0(false, "Unsupported projection type.");
104 break;
105 }
106 }
107
111}
112
114 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
115 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
116{
117 int ndim = m_spacedim;
118 int nvariables = inarray.size();
119 int nq = GetTotPoints();
120
121 switch (m_projectionType)
122 {
124 {
125 //-------------------------------------------------------
126 // Compute the DG advection including the numerical flux
127 // by using SolverUtils/Advection
128 // Input and output in physical space
130 inarray, outarray, time);
131 //-------------------------------------------------------
132
133 //-------------------------------------------------------
134 // negate the outarray since moving terms to the rhs
135 for (int i = 0; i < nvariables; ++i)
136 {
137 Vmath::Neg(nq, outarray[i], 1);
138 }
139 //-------------------------------------------------------
140
141 //-------------------------------------------------
142 // Add "source terms"
143 // Input and output in physical space
144
145 // Coriolis forcing
146 if (m_coriolis.size() != 0)
147 {
148 AddCoriolis(inarray, outarray);
149 }
150
151 // Variable Depth
152 if (m_constantDepth != true)
153 {
154 AddVariableDepth(inarray, outarray);
155 }
156 //-------------------------------------------------
157 }
158 break;
160 {
161 //-------------------------------------------------------
162 // Compute the fluxvector in physical space
164 nvariables);
165
166 for (int i = 0; i < nvariables; ++i)
167 {
168 fluxvector[i] = Array<OneD, Array<OneD, NekDouble>>(ndim);
169 for (int j = 0; j < ndim; ++j)
170 {
171 fluxvector[i][j] = Array<OneD, NekDouble>(nq);
172 }
173 }
174
175 NonlinearSWE::GetFluxVector(inarray, fluxvector);
176 //-------------------------------------------------------
177
178 //-------------------------------------------------------
179 // Take the derivative of the flux terms
180 // and negate the outarray since moving terms to the rhs
181 Array<OneD, NekDouble> tmp0(nq);
182 Array<OneD, NekDouble> tmp1(nq);
183
184 for (int i = 0; i < nvariables; ++i)
185 {
187 fluxvector[i][0], tmp0);
189 fluxvector[i][1], tmp1);
190 Vmath::Vadd(nq, tmp0, 1, tmp1, 1, outarray[i], 1);
191 Vmath::Neg(nq, outarray[i], 1);
192 }
193
194 //-------------------------------------------------
195 // Add "source terms"
196 // Input and output in physical space
197
198 // Coriolis forcing
199 if (m_coriolis.size() != 0)
200 {
201 AddCoriolis(inarray, outarray);
202 }
203
204 // Variable Depth
205 if (m_constantDepth != true)
206 {
207 AddVariableDepth(inarray, outarray);
208 }
209 //-------------------------------------------------
210 }
211 break;
212 default:
213 ASSERTL0(false, "Unknown projection scheme for the NonlinearSWE");
214 break;
215 }
216}
217
219{
221 SolverUtils::AddSummaryItem(s, "Variables", "h should be in field[0]");
222 SolverUtils::AddSummaryItem(s, "", "hu should be in field[1]");
223 SolverUtils::AddSummaryItem(s, "", "hv should be in field[2]");
224}
225
226// Physfield in conservative Form
228 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
230{
231 int nq = m_fields[0]->GetTotPoints();
232
233 NekDouble g = m_g;
235
236 // Flux vector for the mass equation
237 for (int i = 0; i < m_spacedim; ++i)
238 {
239 velocity[i] = Array<OneD, NekDouble>(nq);
240 Vmath::Vcopy(nq, physfield[i + 1], 1, flux[0][i], 1);
241 }
242
243 GetVelocityVector(physfield, velocity);
244
245 // Put (0.5 g h h) in tmp
247 Vmath::Vmul(nq, physfield[0], 1, physfield[0], 1, tmp, 1);
248 Vmath::Smul(nq, 0.5 * g, tmp, 1, tmp, 1);
249
250 // Flux vector for the momentum equations
251 for (int i = 0; i < m_spacedim; ++i)
252 {
253 for (int j = 0; j < m_spacedim; ++j)
254 {
255 Vmath::Vmul(nq, velocity[j], 1, physfield[i + 1], 1, flux[i + 1][j],
256 1);
257 }
258
259 // Add (0.5 g h h) to appropriate field
260 Vmath::Vadd(nq, flux[i + 1][i], 1, tmp, 1, flux[i + 1][i], 1);
261 }
262}
263
264/**
265 * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
266 * \f$ h\mathbf{v} \f$.
267 *
268 * @param physfield Momentum field.
269 * @param velocity Velocity field.
270 */
272 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
274{
275 const int npts = physfield[0].size();
276
277 for (int i = 0; i < m_spacedim; ++i)
278 {
279 Vmath::Vdiv(npts, physfield[1 + i], 1, physfield[0], 1, velocity[i], 1);
280 }
281}
282
283// physarray contains the conservative variables
285 const Array<OneD, const Array<OneD, NekDouble>> &physarray,
287{
288 int ncoeffs = GetNcoeffs();
289 int nq = GetTotPoints();
290
292 Array<OneD, NekDouble> mod(ncoeffs);
293
294 switch (m_projectionType)
295 {
297 {
298 for (int i = 0; i < m_spacedim; ++i)
299 {
300 Vmath::Vmul(nq, m_bottomSlope[i], 1, physarray[0], 1, tmp, 1);
301 Vmath::Smul(nq, m_g, tmp, 1, tmp, 1);
302 m_fields[0]->IProductWRTBase(tmp, mod);
303 m_fields[0]->MultiplyByElmtInvMass(mod, mod);
304 m_fields[0]->BwdTrans(mod, tmp);
305 Vmath::Vadd(nq, tmp, 1, outarray[i + 1], 1, outarray[i + 1], 1);
306 }
307 }
308 break;
310 {
311 for (int i = 0; i < m_spacedim; ++i)
312 {
313 Vmath::Vmul(nq, m_bottomSlope[i], 1, physarray[0], 1, tmp, 1);
314 Vmath::Smul(nq, m_g, tmp, 1, tmp, 1);
315 Vmath::Vadd(nq, tmp, 1, outarray[i + 1], 1, outarray[i + 1], 1);
316 }
317 }
318 break;
319 default:
320 ASSERTL0(false, "Unknown projection scheme for the NonlinearSWE");
321 break;
322 }
323}
324
325} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
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 DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
static std::string className
Name of class.
Definition: NonlinearSWE.h:60
void GetVelocityVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print a summary of time stepping parameters.
void AddVariableDepth(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
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.
Definition: NonlinearSWE.h:49
NonlinearSWE(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void v_InitObject(bool DeclareFields=true) override
Init object for UnsteadySystem class.
void v_DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time) override
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()
void AddCoriolis(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
SolverUtils::AdvectionSharedPtr m_advection
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.
void DoImplicitSolve(const Array< OneD, const Array< OneD, NekDouble > > &inpnts, Array< OneD, Array< OneD, NekDouble > > &outpnt, const NekDouble time, const NekDouble lambda)
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()
int m_spacedim
Spatial dimension (>= expansion dim).
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNcoeffs()
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTotPoints()
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:87
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
double NekDouble
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.
Definition: Vmath.hpp:72
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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.
Definition: Vmath.hpp:180
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
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.
Definition: Vmath.hpp:126
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825