Nektar++
LinearSWE.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: LinearSWE.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: Linear Shallow water equations in primitive variables
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36#include <iostream>
37
40
41namespace Nektar
42{
43
44std::string LinearSWE::className =
46 "LinearSWE", LinearSWE::create,
47 "Linear shallow water equation in primitive variables.");
48
51 : ShallowWaterSystem(pSession, pGraph)
52{
53}
54
55void LinearSWE::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(&LinearSWE::GetFluxVector, this);
81
82 // Setting up Riemann solver for advection operator
83 m_session->LoadSolverInfo("UpwindType", riemName, "NoSolver");
84 if ((riemName == "LinearHLL") && (m_constantDepth != true))
85 {
86 ASSERTL0(false, "LinearHLL only valid for constant depth");
87 }
90 riemName, m_session);
91
92 // Setting up parameters for advection operator Riemann solver
93 m_riemannSolver->SetParam("gravity", &LinearSWE::GetGravity, this);
94 m_riemannSolver->SetAuxVec("vecLocs", &LinearSWE::GetVecLocs, this);
95 m_riemannSolver->SetVector("N", &LinearSWE::GetNormals, this);
96
97 // The numerical flux for linear SWE requires depth information
98 int nTracePointsTot = m_fields[0]->GetTrace()->GetTotPoints();
99 m_dFwd = Array<OneD, NekDouble>(nTracePointsTot);
100 m_dBwd = Array<OneD, NekDouble>(nTracePointsTot);
101 m_fields[0]->GetFwdBwdTracePhys(m_depth, m_dFwd, m_dBwd);
103 m_riemannSolver->SetScalar("depthFwd", &LinearSWE::GetDepthFwd,
104 this);
105 m_riemannSolver->SetScalar("depthBwd", &LinearSWE::GetDepthBwd,
106 this);
107
108 // Concluding initialisation of advection operators
109 m_advection->SetRiemannSolver(m_riemannSolver);
110 m_advection->InitObject(m_session, m_fields);
111 break;
112 }
113 default:
114 {
115 ASSERTL0(false, "Unsupported projection type.");
116 break;
117 }
118 }
119
123}
124
126 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
127 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
128{
129 int ndim = m_spacedim;
130 int nvariables = inarray.size();
131 int nq = GetTotPoints();
132
133 switch (m_projectionType)
134 {
136 {
137 //-------------------------------------------------------
138 // Compute the DG advection including the numerical flux
139 // by using SolverUtils/Advection
140 // Input and output in physical space
142 inarray, outarray, time);
143 //-------------------------------------------------------
144
145 //-------------------------------------------------------
146 // negate the outarray since moving terms to the rhs
147 for (int i = 0; i < nvariables; ++i)
148 {
149 Vmath::Neg(nq, outarray[i], 1);
150 }
151 //-------------------------------------------------------
152
153 //-------------------------------------------------
154 // Add "source terms"
155 // Input and output in physical space
156
157 // Coriolis forcing
158 if (m_coriolis.size() != 0)
159 {
160 AddCoriolis(inarray, outarray);
161 }
162 //-------------------------------------------------
163 }
164 break;
166 {
167 //-------------------------------------------------------
168 // Compute the fluxvector in physical space
170 nvariables);
171
172 for (int i = 0; i < nvariables; ++i)
173 {
174 fluxvector[i] = Array<OneD, Array<OneD, NekDouble>>(ndim);
175 for (int j = 0; j < ndim; ++j)
176 {
177 fluxvector[i][j] = Array<OneD, NekDouble>(nq);
178 }
179 }
180
181 LinearSWE::GetFluxVector(inarray, fluxvector);
182 //-------------------------------------------------------
183
184 //-------------------------------------------------------
185 // Take the derivative of the flux terms
186 // and negate the outarray since moving terms to the rhs
187 Array<OneD, NekDouble> tmp0(nq);
188 Array<OneD, NekDouble> tmp1(nq);
189
190 for (int i = 0; i < nvariables; ++i)
191 {
193 fluxvector[i][0], tmp0);
195 fluxvector[i][1], tmp1);
196 Vmath::Vadd(nq, tmp0, 1, tmp1, 1, outarray[i], 1);
197 Vmath::Neg(nq, outarray[i], 1);
198 }
199
200 //-------------------------------------------------
201 // Add "source terms"
202 // Input and output in physical space
203
204 // Coriolis forcing
205 if (m_coriolis.size() != 0)
206 {
207 AddCoriolis(inarray, outarray);
208 }
209 //-------------------------------------------------
210 }
211 break;
212 default:
213 ASSERTL0(false, "Unknown projection scheme for the LinearSWE");
214 break;
215 }
216}
217
219{
221 if (m_session->DefinesSolverInfo("UpwindType"))
222 {
223 std::string UpwindType;
224 UpwindType = m_session->GetSolverInfo("UpwindType");
225 if (UpwindType == "LinearAverage")
226 {
227 SolverUtils::AddSummaryItem(s, "Riemann Solver", "Linear Average");
228 }
229 else if (UpwindType == "LinearHLL")
230 {
231 SolverUtils::AddSummaryItem(s, "Riemann Solver", "Linear HLL");
232 }
233 }
234 SolverUtils::AddSummaryItem(s, "Variables", "eta should be in field[0]");
235 SolverUtils::AddSummaryItem(s, "", "u should be in field[1]");
236 SolverUtils::AddSummaryItem(s, "", "v should be in field[2]");
237}
238
239// Physfield in primitive Form
241 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
243{
244 int nq = m_fields[0]->GetTotPoints();
245
246 NekDouble g = m_g;
247
248 // Flux vector for the mass equation
249 for (int i = 0; i < m_spacedim; ++i)
250 {
251 Vmath::Vmul(nq, m_depth, 1, physfield[i + 1], 1, flux[0][i], 1);
252 }
253
254 // Flux vector for the momentum equations
255 for (int i = 0; i < m_spacedim; ++i)
256 {
257 for (int j = 0; j < m_spacedim; ++j)
258 {
259 // must zero fluxes as not initialised to zero in AdvectionWeakDG
260 Vmath::Zero(nq, flux[i + 1][j], 1);
261 }
262
263 // Add (g eta) to appropriate field
264 Vmath::Smul(nq, g, physfield[0], 1, flux[i + 1][i], 1);
265 }
266}
267
268/**
269 * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
270 * \f$ h\mathbf{v} \f$.
271 *
272 * @param physfield Velocity field.
273 * @param velocity Velocity field.
274 */
276 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
278{
279 const int npts = physfield[0].size();
280
281 for (int i = 0; i < m_spacedim; ++i)
282 {
283 Vmath::Vcopy(npts, physfield[1 + i], 1, velocity[i], 1);
284 }
285}
286
289{
290 int cnt = 0;
291
292 // loop over Boundary Regions
293 for (int bcRegion = 0; bcRegion < m_fields[0]->GetBndConditions().size();
294 ++bcRegion)
295 {
296 if (m_fields[0]
297 ->GetBndConditions()[bcRegion]
298 ->GetBoundaryConditionType() == SpatialDomains::ePeriodic)
299 {
300 continue;
301 }
302
303 // Copy the forward trace of the field to the backward trace
304 int id2, npts;
305
306 for (int e = 0;
307 e < m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
308 ++e)
309 {
310 npts = m_fields[0]
311 ->GetBndCondExpansions()[bcRegion]
312 ->GetExp(e)
313 ->GetTotPoints();
314 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
315 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt +
316 e));
317
318 Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
319 }
320
321 cnt += m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
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: LinearSWE.h:60
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print a summary of time stepping parameters.
Definition: LinearSWE.cpp:218
LinearSWE(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition: LinearSWE.cpp:49
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: LinearSWE.h:49
const Array< OneD, NekDouble > & GetDepthBwd()
Definition: LinearSWE.h:92
const Array< OneD, NekDouble > & GetDepthFwd()
Definition: LinearSWE.h:87
void v_DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time) override
Definition: LinearSWE.cpp:125
void v_InitObject(bool DeclareFields=true) override
Init object for UnsteadySystem class.
Definition: LinearSWE.cpp:55
void GetVelocityVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
Definition: LinearSWE.cpp:275
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Definition: LinearSWE.cpp:240
Array< OneD, NekDouble > m_dFwd
Still water depth traces.
Definition: LinearSWE.h:99
Array< OneD, NekDouble > m_dBwd
Definition: LinearSWE.h:100
void CopyBoundaryTrace(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: LinearSWE.cpp:287
Base class for unsteady solvers.
NekDouble m_g
Acceleration of gravity.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
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)
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()
int m_spacedim
Spatial dimension (>= expansion dim).
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetExpSize()
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 Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825