Nektar++
LEE.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: LEE.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2017 Kilian Lackhove
10// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
11// Department of Aeronautics, Imperial College London (UK), and Scientific
12// Computing and Imaging Institute, University of Utah (USA).
13//
14// Permission is hereby granted, free of charge, to any person obtaining a
15// copy of this software and associated documentation files (the "Software"),
16// to deal in the Software without restriction, including without limitation
17// the rights to use, copy, modify, merge, publish, distribute, sublicense,
18// and/or sell copies of the Software, and to permit persons to whom the
19// Software is furnished to do so, subject to the following conditions:
20//
21// The above copyright notice and this permission notice shall be included
22// in all copies or substantial portions of the Software.
23//
24// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30// DEALINGS IN THE SOFTWARE.
31//
32// Description: Linearized Euler Equations
33//
34///////////////////////////////////////////////////////////////////////////////
35
36#include <iostream>
37
39
40namespace Nektar
41{
43 "LEE", LEE::create, "Linearized Euler Equations");
44
47 : UnsteadySystem(pSession, pGraph), AcousticSystem(pSession, pGraph)
48{
49 m_ip = 0;
50 m_irho = 1;
51 m_iu = 2;
52
53 m_conservative = true;
54}
55
56/**
57 * @brief Initialization object for the LEE class.
58 */
59void LEE::v_InitObject(bool DeclareFields)
60{
61 AcousticSystem::v_InitObject(DeclareFields);
62
63 m_bfNames.push_back("gamma");
64
65 // Initialize basefield again
67 for (int i = 0; i < m_bf.size(); ++i)
68 {
70 }
71 GetFunction("Baseflow", m_fields[0], true)
72 ->Evaluate(m_bfNames, m_bf, m_time);
73
74 // Define the normal velocity fields
76 for (int i = 0; i < m_bfFwdBwd.size(); i++)
77 {
79 }
80
81 std::string riemName;
82 m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
83 if (boost::to_lower_copy(riemName) == "characteristics" ||
84 boost::to_lower_copy(riemName) == "leeupwind" ||
85 boost::to_lower_copy(riemName) == "upwind")
86 {
87 riemName = "LEEUpwind";
88 }
89 if (boost::to_lower_copy(riemName) == "laxfriedrichs")
90 {
91 riemName = "LEELaxFriedrichs";
92 }
94 riemName, m_session);
95 m_riemannSolver->SetVector("N", &LEE::GetNormals, this);
96 m_riemannSolver->SetVector("basefieldFwdBwd", &LEE::GetBasefieldFwdBwd,
97 this);
98 m_riemannSolver->SetAuxVec("vecLocs", &LEE::GetVecLocs, this);
99
100 // Set up advection operator
101 std::string advName;
102 m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
105 m_advection->SetFluxVector(&LEE::v_GetFluxVector, this);
106 m_advection->SetRiemannSolver(m_riemannSolver);
107 m_advection->InitObject(m_session, m_fields);
108
110 {
113 }
114 else
115 {
116 ASSERTL0(false, "Implicit LEE not set up.");
117 }
118}
119
122{
123 int nq = GetTotPoints();
124
126 for (int i = 0; i < m_spacedim + 2; ++i)
127 {
128 if (i == 1)
129 {
130 // skip rho
131 continue;
132 }
133
134 linTerm[i] = Array<OneD, NekDouble>(nq);
135 }
136
140
141 Array<OneD, NekDouble> gammaMinOne(nq);
142 Vmath::Sadd(nq, -1.0, gamma, 1, gammaMinOne, 1);
143
145 Vmath::Vmul(nq, c0sq, 1, rho0, 1, p0, 1);
146 Vmath::Vdiv(nq, p0, 1, gamma, 1, p0, 1);
147
149 for (int i = 0; i < m_spacedim; ++i)
150 {
151 u0[i] = m_bf[2 + i];
152 }
153
154 Array<OneD, const NekDouble> p = inarray[0];
155 Array<OneD, const NekDouble> rho = inarray[1];
156
158 for (int i = 0; i < m_spacedim; ++i)
159 {
160 ru[i] = inarray[2 + i];
161 }
162
163 Array<OneD, NekDouble> grad(nq);
164 Array<OneD, NekDouble> tmp1(nq);
165
166 // p
167 {
168 Vmath::Zero(nq, linTerm[m_ip], 1);
169 // (1-gamma) ( ru_j / rho0 * dp0/dx_j - p * du0_j/dx_j )
170 for (int j = 0; j < m_spacedim; ++j)
171 {
172 // ru_j / rho0 * dp0/dx_j
173 m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j], p0, grad);
174 Vmath::Vmul(nq, grad, 1, ru[j], 1, tmp1, 1);
175 Vmath::Vdiv(nq, tmp1, 1, rho0, 1, tmp1, 1);
176 // p * du0_j/dx_j - ru_j / rho0 * dp0/dx_j
177 m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j], u0[j],
178 grad);
179 Vmath::Vvtvm(nq, grad, 1, p, 1, tmp1, 1, tmp1, 1);
180 // (gamma-1) (p * du0_j/dx_j - ru_j / rho0 * dp0/dx_j)
181 Vmath::Vvtvp(nq, gammaMinOne, 1, tmp1, 1, linTerm[m_ip], 1,
182 linTerm[m_ip], 1);
183 }
184 }
185
186 // rho has no linTerm
187
188 // ru_i
189 for (int i = 0; i < m_spacedim; ++i)
190 {
191 Vmath::Zero(nq, linTerm[m_iu + i], 1);
192 // du0_i/dx_j * (u0_j * rho + ru_j)
193 for (int j = 0; j < m_spacedim; ++j)
194 {
195 // d u0_i / d x_j
196 m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j], u0[i],
197 grad);
198 // u0_j * rho + ru_j
199 Vmath::Vvtvp(nq, u0[j], 1, rho, 1, ru[j], 1, tmp1, 1);
200 // du0_i/dx_j * (u0_j * rho + ru_j)
201 Vmath::Vvtvp(nq, grad, 1, tmp1, 1, linTerm[m_iu + i], 1,
202 linTerm[m_iu + i], 1);
203 }
204 }
205
207 for (int i = 0; i < m_spacedim + 2; ++i)
208 {
209 if (i == 1)
210 {
211 // skip rho
212 continue;
213 }
214
215 m_fields[0]->FwdTrans(linTerm[i], tmpC);
216 m_fields[0]->BwdTrans(tmpC, linTerm[i]);
217
218 Vmath::Vadd(nq, outarray[i], 1, linTerm[i], 1, outarray[i], 1);
219 }
220}
221
222/**
223 * @brief Return the flux vector for the LEE equations.
224 *
225 * @param physfield Fields.
226 * @param flux Resulting flux. flux[eq][dir][pt]
227 */
229 const Array<OneD, Array<OneD, NekDouble>> &physfield,
231{
232 int nq = physfield[0].size();
233
234 ASSERTL1(flux[0].size() == m_spacedim,
235 "Dimension of flux array and velocity array do not match");
236
239 for (int i = 0; i < m_spacedim; ++i)
240 {
241 u0[i] = m_bf[2 + i];
242 }
243
245 Array<OneD, const NekDouble> rho = physfield[m_irho];
247 for (int i = 0; i < m_spacedim; ++i)
248 {
249 ru[i] = physfield[m_iu + i];
250 }
251
252 // F_{adv,p',j} = c0^2 * ru_j + u0_j * p
253 for (int j = 0; j < m_spacedim; ++j)
254 {
255 int i = 0;
256 Vmath::Vvtvvtp(nq, c0sq, 1, ru[j], 1, u0[j], 1, p, 1, flux[i][j], 1);
257 }
258
259 // F_{adv,rho',j} = u0_j * rho' + ru_j
260 for (int j = 0; j < m_spacedim; ++j)
261 {
262 int i = 1;
263 // u0_j * rho' + ru_j
264 Vmath::Vvtvp(nq, u0[j], 1, rho, 1, ru[j], 1, flux[i][j], 1);
265 }
266
267 for (int i = 0; i < m_spacedim; ++i)
268 {
269 // F_{adv,u'_i,j} = ru_i * u0_j + delta_ij * p
270 for (int j = 0; j < m_spacedim; ++j)
271 {
272 // ru_i * u0_j
273 Vmath::Vmul(nq, ru[i], 1, u0[j], 1, flux[m_iu + i][j], 1);
274
275 // kronecker delta
276 if (i == j)
277 {
278 // delta_ij + p
279 Vmath::Vadd(nq, p, 1, flux[m_iu + i][j], 1, flux[m_iu + i][j],
280 1);
281 }
282 }
283 }
284}
285
286/**
287 * @brief Outflow characteristic boundary conditions for compressible
288 * flow problems.
289 */
290void LEE::v_RiemannInvariantBC(int bcRegion, int cnt,
293 Array<OneD, Array<OneD, NekDouble>> &physarray)
294{
295 int id1, id2, nBCEdgePts;
296 int nVariables = physarray.size();
297
298 const Array<OneD, const int> &traceBndMap = m_fields[0]->GetTraceBndMap();
299
300 int eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
301
302 for (int e = 0; e < eMax; ++e)
303 {
304 nBCEdgePts = m_fields[0]
305 ->GetBndCondExpansions()[bcRegion]
306 ->GetExp(e)
307 ->GetTotPoints();
308 id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
309 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt + e]);
310
311 // Calculate (v.n)
312 Array<OneD, NekDouble> RVn(nBCEdgePts, 0.0);
313 for (int i = 0; i < m_spacedim; ++i)
314 {
315 Vmath::Vvtvp(nBCEdgePts, &Fwd[m_iu + i][id2], 1,
316 &m_traceNormals[i][id2], 1, &RVn[0], 1, &RVn[0], 1);
317 }
318
319 // Calculate (v0.n)
320 Array<OneD, NekDouble> RVn0(nBCEdgePts, 0.0);
321 for (int i = 0; i < m_spacedim; ++i)
322 {
323 Vmath::Vvtvp(nBCEdgePts, &BfFwd[2 + i][id2], 1,
324 &m_traceNormals[i][id2], 1, &RVn0[0], 1, &RVn0[0], 1);
325 }
326
327 for (int i = 0; i < nBCEdgePts; ++i)
328 {
329 NekDouble c = sqrt(BfFwd[0][id2 + i]);
330
331 NekDouble h1, h4, h5;
332
333 if (RVn0[i] > 0)
334 {
335 // rho - p / c^2
336 h1 = Fwd[m_irho][id2 + i] - Fwd[m_ip][id2 + i] / (c * c);
337 }
338 else
339 {
340 h1 = 0.0;
341 }
342
343 if (RVn0[i] - c > 0)
344 {
345 // ru / 2 - p / (2*c)
346 h4 = RVn[i] / 2 - Fwd[m_ip][id2 + i] / (2 * c);
347 }
348 else
349 {
350 h4 = 0.0;
351 }
352
353 if (RVn0[i] + c > 0)
354 {
355 // ru / 2 + p / (2*c)
356 h5 = RVn[i] / 2 + Fwd[m_ip][id2 + i] / (2 * c);
357 }
358 else
359 {
360 h5 = 0.0;
361 }
362
363 // compute conservative variables
364 // p = c0*(h5-h4)
365 // rho = h1 + (h5-h4)/c0
366 // ru = h4+h5
367 Fwd[m_ip][id2 + i] = c * (h5 - h4);
368 Fwd[m_irho][id2 + i] = h1 + (h5 - h4) / c;
369 NekDouble RVnNew = h4 + h5;
370
371 // adjust velocity pert. according to new value
372 // here we just omit the wall parallel velocity components, i.e
373 // setting them to zero. This is equivalent to setting the two
374 // vorticity characteristics h2 and h3 to zero. Mathematically,
375 // this is only legitimate for incoming characteristics. However,
376 // as h2 and h3 are convected by the flow, the value we precribe at
377 // an the boundary for putgoing characteristics does not matter.
378 // This implementation saves a few operations and is more robust
379 // for mixed in/outflow boundaries and at the boundaries edges.
380 for (int j = 0; j < m_spacedim; ++j)
381 {
382 Fwd[m_iu + j][id2 + i] = RVnNew * m_traceNormals[j][id2 + i];
383 }
384 }
385
386 // Copy boundary adjusted values into the boundary expansion
387 for (int i = 0; i < nVariables; ++i)
388 {
389 Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
390 &(m_fields[i]
391 ->GetBndCondExpansions()[bcRegion]
392 ->UpdatePhys())[id1],
393 1);
394 }
395 }
396}
397
398} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the right-hand side.
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
std::vector< std::string > m_bfNames
bool m_conservative
we are dealing with a conservative formualtion
SolverUtils::AdvectionSharedPtr m_advection
int m_ip
indices of the fields
Array< OneD, Array< OneD, NekDouble > > m_bfFwdBwd
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection and call the method for imposing the boundary conditions in case of discontinu...
Array< OneD, Array< OneD, NekDouble > > m_bf
void v_InitObject(bool DeclareFields=true) override
Initialization object for the AcousticSystem class.
const Array< OneD, const Array< OneD, NekDouble > > & GetBasefieldFwdBwd()
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: LEE.h:52
void v_GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux) override
Return the flux vector for the LEE equations.
Definition: LEE.cpp:228
void v_InitObject(bool DeclareFields=true) override
Initialization object for the LEE class.
Definition: LEE.cpp:59
LEE(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition: LEE.cpp:45
static std::string className
Name of class.
Definition: LEE.h:63
void v_AddLinTerm(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
Definition: LEE.cpp:120
void v_RiemannInvariantBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &BfFwd, Array< OneD, Array< OneD, NekDouble > > &physarray) override
Outflow characteristic boundary conditions for compressible flow problems.
Definition: LEE.cpp:290
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)
int m_spacedim
Spatial dimension (>= expansion dim).
NekDouble m_time
Current time of simulation.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNcoeffs()
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 GetTotPoints()
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
Base class for unsteady solvers.
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
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:87
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
EquationSystemFactory & GetEquationSystemFactory()
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
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 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
Definition: Vmath.hpp:366
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 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
Definition: Vmath.hpp:381
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 Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.hpp:439
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285