Nektar++
APE.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: APE.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2018 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: APE1/APE4 (Acoustic Perturbation Equations)
33//
34///////////////////////////////////////////////////////////////////////////////
35
36#include <iostream>
37
39
40namespace Nektar
41{
42
44 "APE", APE::create, "APE1/APE4 (Acoustic Perturbation Equations)");
45
48 : UnsteadySystem(pSession, pGraph), AcousticSystem(pSession, pGraph)
49{
50 m_ip = 0;
51 m_irho = -1;
52 m_iu = 1;
53
54 m_conservative = false;
55}
56
57/**
58 * @brief Initialization object for the APE class.
59 */
60void APE::v_InitObject(bool DeclareFields)
61{
62 AcousticSystem::v_InitObject(DeclareFields);
63
64 // Initialize basefield again
66 for (int i = 0; i < m_bf.size(); ++i)
67 {
69 }
70 GetFunction("Baseflow", m_fields[0], true)
71 ->Evaluate(m_bfNames, m_bf, m_time);
72
73 // Define the normal velocity fields
75 for (int i = 0; i < m_bfFwdBwd.size(); i++)
76 {
78 }
79
80 std::string riemName;
81 m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
82 if (boost::to_lower_copy(riemName) == "characteristics" ||
83 boost::to_lower_copy(riemName) == "apeupwind" ||
84 boost::to_lower_copy(riemName) == "upwind")
85 {
86 riemName = "APEUpwind";
87 }
88 if (boost::to_lower_copy(riemName) == "laxfriedrichs")
89 {
90 riemName = "APELaxFriedrichs";
91 }
93 riemName, m_session);
94 m_riemannSolver->SetVector("N", &APE::GetNormals, this);
95 m_riemannSolver->SetVector("basefieldFwdBwd", &APE::GetBasefieldFwdBwd,
96 this);
97 m_riemannSolver->SetAuxVec("vecLocs", &APE::GetVecLocs, this);
98
99 // Set up advection operator
100 std::string advName;
101 m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
104 m_advection->SetFluxVector(&APE::v_GetFluxVector, this);
105 m_advection->SetRiemannSolver(m_riemannSolver);
106 m_advection->InitObject(m_session, m_fields);
107
109 {
112 }
113 else
114 {
115 ASSERTL0(false, "Implicit APE not set up.");
116 }
117}
118
119/**
120 * @brief Return the flux vector for the APE equations.
121 *
122 * @param physfield Fields.
123 * @param flux Resulting flux. flux[eq][dir][pt]
124 */
126 const Array<OneD, Array<OneD, NekDouble>> &physfield,
128{
129 int nq = physfield[0].size();
130 Array<OneD, NekDouble> tmp1(nq);
131 Array<OneD, NekDouble> tmp2(nq);
132
133 ASSERTL1(flux[0].size() == m_spacedim,
134 "Dimension of flux array and velocity array do not match");
135
136 // F_{adv,p',j} = \bar{rho} \bar{c^2} u'_j + p' \bar{u}_j
137 for (int j = 0; j < m_spacedim; ++j)
138 {
139 Vmath::Zero(nq, flux[0][j], 1);
140
141 // construct \bar{rho} \bar{c^2} u'_j
142 Vmath::Vmul(nq, m_bf[0], 1, m_bf[1], 1, tmp1, 1);
143 Vmath::Vmul(nq, tmp1, 1, physfield[j + 1], 1, tmp1, 1);
144
145 // construct p' \bar{u}_j term
146 Vmath::Vmul(nq, physfield[0], 1, m_bf[j + 2], 1, tmp2, 1);
147
148 // add both terms
149 Vmath::Vadd(nq, tmp1, 1, tmp2, 1, flux[0][j], 1);
150 }
151
152 for (int i = 1; i < flux.size(); ++i)
153 {
154 ASSERTL1(flux[i].size() == m_spacedim,
155 "Dimension of flux array and velocity array do not match");
156
157 // F_{adv,u'_i,j} = (p'/ \bar{rho} + \bar{u}_k u'_k) \delta_{ij}
158 for (int j = 0; j < m_spacedim; ++j)
159 {
160 Vmath::Zero(nq, flux[i][j], 1);
161
162 if (i - 1 == j)
163 {
164 // contruct p'/ \bar{rho} term
165 Vmath::Vdiv(nq, physfield[0], 1, m_bf[1], 1, flux[i][j], 1);
166
167 // construct \bar{u}_k u'_k term
168 Vmath::Zero(nq, tmp1, 1);
169 for (int k = 0; k < m_spacedim; ++k)
170 {
171 Vmath::Vvtvp(nq, physfield[k + 1], 1, m_bf[k + 2], 1, tmp1,
172 1, tmp1, 1);
173 }
174
175 // add terms
176 Vmath::Vadd(nq, flux[i][j], 1, tmp1, 1, flux[i][j], 1);
177 }
178 }
179 }
180}
181
182/**
183 * @brief Outflow characteristic boundary conditions for compressible
184 * flow problems.
185 */
186void APE::v_RiemannInvariantBC(int bcRegion, int cnt,
189 Array<OneD, Array<OneD, NekDouble>> &physarray)
190{
191 int id1, id2, nBCEdgePts;
192 int nVariables = physarray.size();
193
194 const Array<OneD, const int> &traceBndMap = m_fields[0]->GetTraceBndMap();
195
196 int eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
197
198 for (int e = 0; e < eMax; ++e)
199 {
200 nBCEdgePts = m_fields[0]
201 ->GetBndCondExpansions()[bcRegion]
202 ->GetExp(e)
203 ->GetTotPoints();
204 id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
205 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt + e]);
206
207 // Calculate (v.n)
208 Array<OneD, NekDouble> Vn(nBCEdgePts, 0.0);
209 for (int i = 0; i < m_spacedim; ++i)
210 {
211 Vmath::Vvtvp(nBCEdgePts, &Fwd[m_iu + i][id2], 1,
212 &m_traceNormals[i][id2], 1, &Vn[0], 1, &Vn[0], 1);
213 }
214
215 // Calculate (v0.n)
216 Array<OneD, NekDouble> Vn0(nBCEdgePts, 0.0);
217 for (int i = 0; i < m_spacedim; ++i)
218 {
219 Vmath::Vvtvp(nBCEdgePts, &BfFwd[2 + i][id2], 1,
220 &m_traceNormals[i][id2], 1, &Vn0[0], 1, &Vn0[0], 1);
221 }
222
223 for (int i = 0; i < nBCEdgePts; ++i)
224 {
225 NekDouble c = sqrt(BfFwd[0][id2 + i]);
226
227 // LODI
228 NekDouble h1, h2;
229
230 // outgoing
231 if (Vn0[i] - c > 0)
232 {
233 // u/2 - p/(2*rho0*sqr(c0sq))
234 h1 = Vn[i] / 2 -
235 Fwd[m_ip][id2 + i] / (2 * BfFwd[1][id2 + i] * c);
236 }
237 // incoming
238 else
239 {
240 h1 = 0.0;
241 }
242
243 // outgoing
244 if (Vn0[i] + c > 0)
245 {
246 // u/2 + p/(2*rho0*sqr(c0sq))
247 h2 = Vn[i] / 2 +
248 Fwd[m_ip][id2 + i] / (2 * BfFwd[1][id2 + i] * c);
249 }
250 // incoming
251 else
252 {
253 h2 = 0.0;
254 }
255
256 // compute primitive variables
257 // p = rho0*sqr(c0sq) * (h2 - h1)
258 Fwd[m_ip][id2 + i] = BfFwd[1][id2 + i] * c * (h2 - h1);
259 // u = h1 + h2
260 NekDouble VnNew = h1 + h2;
261
262 // adjust velocity pert. according to new value
263 for (int j = 0; j < m_spacedim; ++j)
264 {
265 Fwd[m_iu + j][id2 + i] =
266 Fwd[m_iu + j][id2 + i] +
267 (VnNew - Vn[i]) * m_traceNormals[j][id2 + i];
268 }
269 }
270
271 // Copy boundary adjusted values into the boundary expansion
272 for (int i = 0; i < nVariables; ++i)
273 {
274 Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
275 &(m_fields[i]
276 ->GetBndCondExpansions()[bcRegion]
277 ->UpdatePhys())[id1],
278 1);
279 }
280 }
281}
282
283} // 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
APE(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition: APE.cpp:46
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 APE equations.
Definition: APE.cpp:125
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: APE.h:52
void v_InitObject(bool DeclareFields=true) override
Initialization object for the APE class.
Definition: APE.cpp:60
static std::string className
Name of class.
Definition: APE.h:63
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: APE.cpp:186
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
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.
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
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 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 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