Nektar++
EnforceEntropyPressure.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: EnforceEntropyPressure.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// License for the specific language governing rights and limitations under
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: Modified Riemann invariant boundary condition.
33// Enforce the entropy and pressure at the inflow boundary;
34// Enforce the pressure at the outflow boundary.
35// The input can be either VALUE or FILE.
36//
37///////////////////////////////////////////////////////////////////////////////
38
40
41namespace Nektar
42{
43
46 "EnforceEntropyPressure", EnforceEntropyPressure::create,
47 "Riemann boundary condition enforcing entropy and pressure.");
48
52 const Array<OneD, Array<OneD, NekDouble>> &pTraceNormals,
53 const Array<OneD, Array<OneD, NekDouble>> &pGridVelocity,
54 const int pSpaceDim, const int bcRegion, const int cnt)
55 : CFSBndCond(pSession, pFields, pTraceNormals, pGridVelocity, pSpaceDim,
56 bcRegion, cnt)
57{
58
60 m_fields[0]->GetBndCondExpansions()[m_bcRegion];
61
62 //-> Gather a list of index from trace to this boundary
63 m_npts = bndexp->GetTotPoints();
64
66
67 const Array<OneD, const int> &traceBndMap = m_fields[0]->GetTraceBndMap();
68
69 // Construct a map for the boundary to trace map for easy acess to
70 // phys space points
71 int cnt1 = 0;
72 for (int e = 0; e < bndexp->GetNumElmts(); ++e)
73 {
74 int nTracePts = bndexp->GetExp(e)->GetTotPoints();
75
76 int id =
77 m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[m_offset + e]);
78
79 // Loop on the points of the m_bcRegion
80 for (int i = 0; i < nTracePts; i++)
81 {
82 // the ith point in region e
83 m_bndToTraceMap[cnt1++] = id + i;
84 }
85 }
86
89
90 for (int i = 0; i < m_fields.size(); ++i)
91 {
92 m_bndPhys[i] =
93 m_fields[i]->GetBndCondExpansions()[m_bcRegion]->UpdatePhys();
94
95 BCvals[i] = Array<OneD, NekDouble>(m_npts);
96 Vmath::Vcopy(m_npts, m_bndPhys[i], 1, BCvals[i], 1);
97 }
98
99 // Set up boudnary required BCs
101 Vmath::Vcopy(m_npts, BCvals[0], 1, m_rhoBC, 1);
102
104 // Evaluate velocity on boundary
105 for (int i = 0; i < m_spacedim; ++i)
106 {
108 Vmath::Vcopy(m_npts, BCvals[i + 1], 1, m_velBC[i], 1);
109 Vmath::Vdiv(m_npts, m_velBC[i], 1, m_rhoBC, 1, m_velBC[i], 1);
110 }
112 m_varConv->GetPressure(BCvals, m_pBC);
113
114 // Computing the normal velocity for characteristics coming
115 // from outside the computational domain
117 for (int i = 0; i < m_spacedim; i++)
118 {
119 for (int j = 0; j < m_npts; ++j)
120 {
121 m_VnInf[j] += m_traceNormals[i][m_bndToTraceMap[j]] * m_velBC[i][j];
122 }
123 }
124}
125
128 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &physarray,
129 [[maybe_unused]] const NekDouble &time)
130{
131 int i, j;
132 int nDimensions = m_spacedim;
133
134 Array<OneD, Array<OneD, NekDouble>> FwdBnd(Fwd.size());
135 Array<OneD, Array<OneD, NekDouble>> bndPhys(Fwd.size());
136
137 // make a local copy of Fwd along boundary of interest
138 for (i = 0; i < Fwd.size(); ++i)
139 {
140 FwdBnd[i] = Array<OneD, NekDouble>(m_npts);
141 for (j = 0; j < m_npts; ++j)
142 {
143 FwdBnd[i][j] = Fwd[i][m_bndToTraceMap[j]];
144 }
145 }
146
147 // Computing the normal velocity for characteristics coming
148 // from inside the computational domain
150
151 for (i = 0; i < nDimensions; ++i)
152 {
153 for (j = 0; j < m_npts; ++j)
154 {
155 Vn[j] += m_traceNormals[i][m_bndToTraceMap[j]] * FwdBnd[i + 1][j];
156 }
157 }
158 // divide by density.
159 Vmath::Vdiv(m_npts, Vn, 1, FwdBnd[0], 1, Vn, 1);
160
161 // Get speed of sound
163 Array<OneD, NekDouble> soundSpeed(m_npts);
164
165 m_varConv->GetPressure(FwdBnd, pressure);
166 m_varConv->GetSoundSpeed(FwdBnd, soundSpeed);
167
168 // Get Mach. Note: it is computed by Vn/c
170 Vmath::Vdiv(m_npts, Vn, 1, soundSpeed, 1, Mach, 1);
171 Vmath::Vabs(m_npts, Mach, 1, Mach, 1);
172
173 // Auxiliary variables
174 Array<OneD, NekDouble> velBC(nDimensions, 0.0);
175
176 // L represents properties outside boundary
177 // R represents properties inside boundary (numerical state)
178 NekDouble rhoL, uL, pL;
179 NekDouble EBC, rR, cstar, pstar, rhostar, ustar; // vn
180
181 NekDouble gamMinOne = m_gamma - 1.0;
182 NekDouble twoOverGamMinOne = 2.0 / gamMinOne;
183 NekDouble gamInv = 1.0 / m_gamma;
184
185 // Loop on m_bcRegions
186 for (int pnt = 0; pnt < m_npts; ++pnt)
187 {
188 // Impose inflow Riemann invariant
189 if (Vn[pnt] <= 0.0)
190 {
191 // Subsonic flows
192 if (Mach[pnt] < 1.00)
193 {
194 // right characteristic
195 rR = -Vn[pnt] - sqrt(m_gamma * pressure[pnt] / FwdBnd[0][pnt]) *
196 twoOverGamMinOne;
197 // vn = -m_VnInf[pnt]; //vn BC
198
199 // fix rhostar and pstar to be the input values
200 // compute ustar using left-pointing characteristic line IR^-
201 pstar = m_pBC[pnt];
202 rhostar = m_rhoBC[pnt];
203 cstar = sqrt(m_gamma * pstar / rhostar);
204 ustar = rR + cstar * twoOverGamMinOne;
205
206 // add supplement equation that rhoL=rhostar
207 // then pL=pstar, according to IL^0
208 // and uL=ustar, according to IL^+
209 rhoL = rhostar;
210 pL = pstar;
211 uL = ustar;
212
213 // std subsnoic inflow
214 // rhoL = m_rhoBC[pnt];
215 // uL = -m_VnInf[pnt];
216 // pL = m_pBC[pnt];
217 }
218 else // Supersonic inflow
219 {
220 // all characteristics are from left so just impose
221 // star state to left values
222 // Note: m_vnInf is the negative of the normal velocity
223 // across boundary
224 rhoL = m_rhoBC[pnt];
225 uL = -m_VnInf[pnt];
226 pL = m_pBC[pnt];
227 }
228
229 // Boundary energy
230 EBC = pL * twoOverGamMinOne * 0.5;
231
232 // evaluate the different between the left state normal
233 // velocity and that from the desired condition (note
234 // m_VnInf is using an outwards normal definition.
235 NekDouble VnDiff = uL + m_VnInf[pnt];
236
237 // Boundary velocities & Kinite energy
238 // Note: normals are negated since they point outwards in
239 // the domain
240
241 // Note: Can just use the BC values directly!!
242 for (j = 0; j < nDimensions; ++j)
243 {
244 // Set velocity to the desired conditions modified to
245 // take account of the normal state for Riemann
246 // problem. (Negative accounts for outwards normal definition)
247 // velBC[j] = m_velBC[j][pnt];
248 velBC[j] = m_velBC[j][pnt] -
249 VnDiff * m_traceNormals[j][m_bndToTraceMap[pnt]];
250
251 EBC += 0.5 * rhoL * velBC[j] * velBC[j];
252 }
253
254 //-------------------------------------------------------------------------
255 // Impose Left hand Riemann Invariant boundary conditions
256 m_bndPhys[0][pnt] = rhoL;
257 for (j = 0; j < nDimensions; ++j)
258 {
259 m_bndPhys[j + 1][pnt] = rhoL * velBC[j];
260 }
261 m_bndPhys[nDimensions + 1][pnt] = EBC;
262 }
263 else // Outflow
264 {
265
266 // Note: Allowing the switch can cause worse convergence in this
267 // type BC.
268 // So improve it later.
269 if (Mach[pnt] < 1.00)
270 {
271 // subsonic outflow: fix pstar
272 rR = -Vn[pnt] - sqrt(m_gamma * pressure[pnt] / FwdBnd[0][pnt]) *
273 twoOverGamMinOne;
274 // vn = -m_VnInf[pnt];
275
276 pstar = m_pBC[pnt];
277 rhostar = FwdBnd[0][pnt] * pow((pstar / pressure[pnt]), gamInv);
278 cstar = sqrt(m_gamma * pstar / rhostar);
279 ustar = rR + cstar * twoOverGamMinOne;
280
281 rhoL = rhostar;
282 uL = ustar;
283 pL = pstar;
284 }
285 else
286 {
287 // supersonic outflow
288 // Just set to imposed state and let Riemann BC dictate values
289 rhoL = m_rhoBC[pnt];
290 uL = -m_VnInf[pnt];
291 pL = m_pBC[pnt];
292 }
293
294 // Boundary energy
295 EBC = pL * twoOverGamMinOne * 0.5;
296
297 // Boundary velocities & Kinite energy
298 // Note: normals are negated since they point outwards in
299 // the domain
300 for (j = 0; j < nDimensions; ++j)
301 {
302 velBC[j] = -1.0 * uL * m_traceNormals[j][m_bndToTraceMap[pnt]];
303 EBC += 0.5 * rhoL * velBC[j] * velBC[j];
304 }
305
306 // Impose Left hand Riemann Invariant boundary conditions
307 m_bndPhys[0][pnt] = rhoL;
308 for (j = 0; j < nDimensions; ++j)
309 {
310 m_bndPhys[j + 1][pnt] = rhoL * velBC[j];
311 }
312 m_bndPhys[nDimensions + 1][pnt] = EBC;
313 }
314 }
315}
316
317} // namespace Nektar
Encapsulates the user-defined boundary conditions for compressible flow solver.
Definition: CFSBndCond.h:72
int m_spacedim
Space dimension.
Definition: CFSBndCond.h:98
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Trace normals.
Definition: CFSBndCond.h:94
NekDouble m_gamma
Parameters of the flow.
Definition: CFSBndCond.h:105
int m_bcRegion
Id of the boundary region.
Definition: CFSBndCond.h:113
VariableConverterSharedPtr m_varConv
Auxiliary object to convert variables.
Definition: CFSBndCond.h:100
int m_offset
Offset.
Definition: CFSBndCond.h:115
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array of fields.
Definition: CFSBndCond.h:92
static std::string className
Name of the class.
void v_Apply(Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray, const NekDouble &time) override
Array< OneD, Array< OneD, NekDouble > > m_velBC
static CFSBndCondSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &pTraceNormals, const Array< OneD, Array< OneD, NekDouble > > &pGridVelocity, const int pSpaceDim, const int bcRegion, const int cnt)
Creates an instance of this class.
Array< OneD, NekDouble > m_VnInf
Reference normal velocity.
Array< OneD, NekDouble > m_rhoBC
Array< OneD, Array< OneD, NekDouble > > m_bndPhys
EnforceEntropyPressure(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &pTraceNormals, const Array< OneD, Array< OneD, NekDouble > > &pGridVelocity, const int pSpaceDim, const int bcRegion, const int cnt)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
CFSBndCondFactory & GetCFSBndCondFactory()
Declaration of the boundary condition factory singleton.
Definition: CFSBndCond.cpp:40
double NekDouble
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.hpp:352
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
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285