MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
dgcartesiansyseqnacousticperturb.h
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
7#ifndef DGSYSEQNACOUSTICPERTURB_H_
8#define DGSYSEQNACOUSTICPERTURB_H_
9
10#include <algorithm>
11#include <array>
12#include <iterator>
13#include <limits>
14#include <numeric>
15#include "INCLUDE/maiaconstants.h"
16#include "INCLUDE/maiatypes.h"
17#include "UTIL/tensor.h"
18#include "dgcartesiansyseqn.h"
19#include "filter.h"
20
21// Only needed (really: possible to use) in single-threaded applications
22#ifndef _OPENMP
23#include "MEMORY/scratch.h"
24#endif
25
26template <MInt nDim>
27class DgSysEqnAcousticPerturb : public DgSysEqn<nDim, DgSysEqnAcousticPerturb<nDim>> {
28 // Declare parent a friend so that CRTP can access private methods/members
29 friend class DgSysEqn<nDim, DgSysEqnAcousticPerturb>;
30
31 // Member functions
32 public:
33 explicit DgSysEqnAcousticPerturb(MInt solverId);
34 void calcInitialCondition(const MFloat t, const MFloat* x, MFloat* const nodeVars, MFloat* u) const;
35 void calcFlux(const MFloat* const nodeVars, const MFloat* const q, const MInt noNodes1D, MFloat* const flux) const;
36 void calcSource(const MFloat* const nodeVars, const MFloat* const u, const MInt noNodes1D, const MFloat t,
37 const MFloat* const x, MFloat* const src) const;
38 void calcSpongeSource(const MFloat* nodeVars, const MFloat* u, const MInt noNodes1D, const MFloat* eta,
39 MFloat* src) const;
40 MFloat getTimeStep(const MFloat* nodeVars, const MFloat* const u, const MInt noNodes1D, const MFloat invJacobian,
41 const MInt sbpMode) const;
42 void calcRiemann(const MFloat* nodeVarsL, const MFloat* nodeVarsR, const MFloat* stateL, const MFloat* stateR,
43 const MInt noNodes1D, const MInt dirId, MFloat* flux) const;
44 void calcRiemannLaxFriedich(const MFloat* nodeVarsL, const MFloat* nodeVarsR, const MFloat* stateL,
45 const MFloat* stateR, const MInt noNodes1D, const MInt dirId, MFloat* flux) const;
46 void calcRiemannRoe(const MFloat* nodeVarsL, const MFloat* nodeVarsR, const MFloat* stateL, const MFloat* stateR,
47 const MInt noNodes1D, const MInt dirId, MFloat* flux) const;
48
49 void primToCons(const MFloat* prim, MFloat* cons) const;
50 void consToPrim(const MFloat* cons, MFloat* prim) const;
51
52 void getDefaultNodeVars(MFloat* const nodeVars) const;
53 void getDefaultNodeVarsBody(MFloat* const nodeVars) const;
54 MBool extendNodeVar(const MInt varId) const;
55
56 private:
57 void calcFlux1D(const MFloat* const nodeVars, const MFloat* const q, const MInt noNodes1D, const MInt dirId,
58 MFloat* const flux) const;
59
60 MFloat cflFactor(const MInt polyDeg) const;
61
62 // Member variables
63 private:
64 static const MString s_sysEqnName;
65
66 static const MInt s_noVariables = nDim + 1;
67 static const MString s_consVarNames[s_noVariables];
68 static const MString s_primVarNames[s_noVariables];
69
70 // Node vars: U0, V0, W0, RHO0, C0, DC0_DX, DC0_DY, DC0_DZ
71 static const MInt s_noNodeVars = 2 * nDim + 2;
72 static const MBool s_hasTimeDependentNodeVars = false;
73 static const MString s_nodeVarNames[s_noNodeVars];
74
75 std::array<MFloat, nDim> m_meanVelocity;
76 MFloat m_meanDensity;
77 MFloat m_meanSpeedOfSound;
78 MBool m_constantSpeedOfSound = false;
79 MFloat m_spongePressureInfy;
80 MFloat m_spongeSigma;
81 MBool m_compressibleSourceTerm = true;
82 MString m_dgIntegrationMethod;
83 MInt m_dgTimeIntegrationScheme;
84 static const MInt s_maxPolyDeg = 31;
85 std::array<std::array<MFloat, s_maxPolyDeg + 1>, 2> m_cflFactor;
86
87 public:
88 // Hold indices for primitive and conservative variables
89 struct CV;
90};
91
92// Initialize static member variables
93template <MInt nDim>
94const MString DgSysEqnAcousticPerturb<nDim>::s_sysEqnName = "DG_SYSEQN_ACOUSTICPERTURB";
95
96// For variable names, see dgsyseqnacousticperturb.cpp
97
98// Use these classes to access position of member variables
99template <>
100struct DgSysEqnAcousticPerturb<2>::CV {
101 static constexpr const MInt nDim = 2;
102
103 // Conservative variables
104 static constexpr const MInt U = 0;
105 static constexpr const MInt V = 1;
106 static constexpr const MInt UU[nDim] = {0, 1};
107 static constexpr const MInt P = nDim;
108
109 // Node variables
110 static constexpr const MInt U0 = 0;
111 static constexpr const MInt V0 = 1;
112 static constexpr const MInt UU0[nDim] = {0, 1};
113 static constexpr const MInt RHO0 = 2;
114 static constexpr const MInt C0 = 3;
115 static constexpr const MInt DC0_DX = 4;
116 static constexpr const MInt DC0_DY = 5;
117 static constexpr const MInt DC0[nDim] = {4, 5};
118};
119
120template <>
121struct DgSysEqnAcousticPerturb<3>::CV {
122 static constexpr const MInt nDim = 3;
123
124 // Conservative variables
125 static constexpr const MInt U = 0;
126 static constexpr const MInt V = 1;
127 static constexpr const MInt W = 2;
128 static constexpr const MInt UU[nDim] = {0, 1, 2};
129 static constexpr const MInt P = nDim;
130
131 // Node variables
132 static constexpr const MInt U0 = 0;
133 static constexpr const MInt V0 = 1;
134 static constexpr const MInt W0 = 2;
135 static constexpr const MInt UU0[nDim] = {0, 1, 2};
136 static constexpr const MInt RHO0 = 3;
137 static constexpr const MInt C0 = 4;
138 static constexpr const MInt DC0_DX = 5;
139 static constexpr const MInt DC0_DY = 6;
140 static constexpr const MInt DC0_DZ = 7;
141 static constexpr const MInt DC0[nDim] = {5, 6, 7};
142};
143
144
154template <MInt nDim>
155DgSysEqnAcousticPerturb<nDim>::DgSysEqnAcousticPerturb(MInt solverId)
156 : DgSysEqn<nDim, DgSysEqnAcousticPerturb>(solverId) {
157 // Read and set mean velocity
158 for(MInt i = 0; i < nDim; i++) {
159 m_meanVelocity[i] = 0.0;
160 m_meanVelocity[i] =
161 Context::getSolverProperty<MFloat>("meanVelocity", this->m_solverId, AT_, &m_meanVelocity[i], i);
162 }
163
176 m_meanDensity = 1.0;
177 m_meanDensity = Context::getSolverProperty<MFloat>("meanDensity", this->m_solverId, AT_, &m_meanDensity);
178
194 m_meanSpeedOfSound = 1.0;
195 m_meanSpeedOfSound =
196 Context::getSolverProperty<MFloat>("meanSpeedOfSound", this->m_solverId, AT_, &m_meanSpeedOfSound);
197
207 m_constantSpeedOfSound = false;
208 m_constantSpeedOfSound =
209 Context::getSolverProperty<MBool>("constantSpeedOfSound", this->m_solverId, AT_, &m_constantSpeedOfSound);
210
211 m_log << "APE default mean variables: ";
212 for(MInt i = 0; i < nDim; i++) {
213 m_log << "u_" << i << " = " << m_meanVelocity[i] << ", ";
214 }
215 m_log << "rho0 = " << m_meanDensity << ", c0 = " << m_meanSpeedOfSound
216 << " (c0 = constant: " << m_constantSpeedOfSound << ")" << std::endl;
217
230 m_spongePressureInfy = 0.0;
231 m_spongePressureInfy =
232 Context::getSolverProperty<MFloat>("spongePressureInfy", this->m_solverId, AT_, &m_spongePressureInfy);
233
246 m_spongeSigma = 1.0;
247 m_spongeSigma = Context::getSolverProperty<MFloat>("spongeSigma", this->m_solverId, AT_, &m_spongeSigma);
248
249 // Enables compressible source term (due to recasting the APE in conservative form)
250 m_compressibleSourceTerm = true;
251 m_compressibleSourceTerm =
252 Context::getSolverProperty<MBool>("compressibleSourceTerm", this->m_solverId, AT_, &m_compressibleSourceTerm);
253
254 const MString statusMsg = (m_compressibleSourceTerm) ? "enabled" : "disabled";
255 m_log << "APE: compressible source term " << statusMsg << std::endl;
256
257 // Get the quadrature method being used
258 MString dgIntegrationMethod = "DG_INTEGRATE_GAUSS";
259 m_dgIntegrationMethod =
260 Context::getSolverProperty<MString>("dgIntegrationMethod", this->m_solverId, AT_, &dgIntegrationMethod);
261
262 // Get the time integration scheme being used (USED WITH SWITCH)
263 MString dgTimeIntegrationScheme = "DG_TIMEINTEGRATION_CARPENTER_4_5";
264 m_dgTimeIntegrationScheme = string2enum(
265 Context::getSolverProperty<MString>("dgTimeIntegrationScheme", this->m_solverId, AT_, &dgTimeIntegrationScheme));
266
267 // cfl correction factor table - Obtained by running different test cases
268 // to get the maximum stable CFL number solution for polynomial degrees
269 // varying from 0 to 15. 1st row -> 2D cases & 2nd row -> 3D cases.
270 // Note: The cfl factors are not calibrated and thus not suitable to be used
271 // with RiemannRoe's Flux.
272 switch(m_dgTimeIntegrationScheme) {
273 case DG_TIMEINTEGRATION_CARPENTER_4_5: {
274 default:
275 if(m_dgIntegrationMethod == "DG_INTEGRATE_GAUSS") {
276 // Legendre-Gauss nodes
277 m_cflFactor = {{{{2.497559, 2.390137, 1.978760, 1.652832, 1.413574, 1.230469, 1.097168, 0.984375, 0.890625,
278 0.812988, 0.747070, 0.691406, 0.643066, 0.602051, 0.565430, 0.533203}},
279 {{2.976074, 2.574463, 2.055664, 1.697998, 1.444092, 1.251221, 1.123047, 1.003418, 0.904541,
280 0.827637, 0.758057, 0.700684, 0.650635, 0.607910, 0.570068, 0.537109}}}};
281 } else {
282 // Legendre-Gauss-Lobatto nodes
283 m_cflFactor = {{{{5.916748, 5.916748, 3.868408, 2.753906, 2.116699, 1.722412, 1.458984, 1.259766, 1.107422,
284 0.990234, 0.895020, 0.815918, 0.750000, 0.692871, 0.644531, 0.602051}},
285 {{7.445068, 7.445068, 4.113770, 2.825928, 2.176514, 1.760254, 1.502686, 1.293945, 1.131592,
286 1.008301, 0.906982, 0.823975, 0.756836, 0.698242, 0.649414, 0.606689}}}};
287 }
288 break;
289 }
290 case DG_TIMEINTEGRATION_TOULORGEC_4_8: {
291 // Using CFL factors used for CARPENTER 4/5
292 if(m_dgIntegrationMethod == "DG_INTEGRATE_GAUSS") {
293 // Legendre-Gauss nodes
294 m_cflFactor = {{{{2.497559, 2.390137, 1.978760, 1.652832, 1.413574, 1.230469, 1.097168, 0.984375, 0.890625,
295 0.812988, 0.747070, 0.691406, 0.643066, 0.602051, 0.565430, 0.533203}},
296 {{2.976074, 2.574463, 2.055664, 1.697998, 1.444092, 1.251221, 1.123047, 1.003418, 0.904541,
297 0.827637, 0.758057, 0.700684, 0.650635, 0.607910, 0.570068, 0.537109}}}};
298 } else {
299 // Legendre-Gauss-Lobatto nodes
300 m_cflFactor = {{{{5.916748, 5.916748, 3.868408, 2.753906, 2.116699, 1.722412, 1.458984, 1.259766, 1.107422,
301 0.990234, 0.895020, 0.815918, 0.750000, 0.692871, 0.644531, 0.602051}},
302 {{7.445068, 7.445068, 4.113770, 2.825928, 2.176514, 1.760254, 1.502686, 1.293945, 1.131592,
303 1.008301, 0.906982, 0.823975, 0.756836, 0.698242, 0.649414, 0.606689}}}};
304 }
305 break;
306 }
307 case DG_TIMEINTEGRATION_NIEGEMANN_4_14: {
308 if(m_dgIntegrationMethod == "DG_INTEGRATE_GAUSS") {
309 // Legendre-Gauss nodes
310 m_cflFactor = {{{{9.669311, 9.232117, 7.577270, 6.309753, 5.397094, 4.675781, 4.187561, 3.739929, 3.382751,
311 3.085876, 2.836548, 2.625488, 2.443419, 2.283385, 2.144226, 2.021301}},
312 {{13.695678, 10.427733, 7.956481, 6.546325, 5.566406, 4.791747, 4.311645, 3.826904, 3.447692,
313 3.163574, 2.896850, 2.667236, 2.480529, 2.314696, 2.169739, 2.042175}}}};
314 } else {
315 // Legendre-Gauss-Lobatto nodes
316 m_cflFactor = {{{{18.998840, 18.998840, 12.540649, 9.887329, 7.777893, 6.397888, 5.569885, 4.792907, 4.213073,
317 3.758484, 3.392029, 3.090515, 2.840026, 2.625488, 2.442260, 2.283385}},
318 {{18.998840, 18.998840, 17.312683, 11.183837, 8.501526, 6.825806, 5.831970, 5.017882, 4.361510,
319 3.876770, 3.475524, 3.154296, 2.889892, 2.662597, 2.473572, 2.307739}}}};
320 }
321 break;
322 }
323 case DG_TIMEINTEGRATION_NIEGEMANN_4_13: {
324 // Using CFL factors used for CARPENTER 4/5
325 if(m_dgIntegrationMethod == "DG_INTEGRATE_GAUSS") {
326 // Legendre-Gauss nodes
327 m_cflFactor = {{{{2.497559, 2.390137, 1.978760, 1.652832, 1.413574, 1.230469, 1.097168, 0.984375, 0.890625,
328 0.812988, 0.747070, 0.691406, 0.643066, 0.602051, 0.565430, 0.533203}},
329 {{2.976074, 2.574463, 2.055664, 1.697998, 1.444092, 1.251221, 1.123047, 1.003418, 0.904541,
330 0.827637, 0.758057, 0.700684, 0.650635, 0.607910, 0.570068, 0.537109}}}};
331 } else {
332 // Legendre-Gauss-Lobatto nodes
333 m_cflFactor = {{{{5.916748, 5.916748, 3.868408, 2.753906, 2.116699, 1.722412, 1.458984, 1.259766, 1.107422,
334 0.990234, 0.895020, 0.815918, 0.750000, 0.692871, 0.644531, 0.602051}},
335 {{7.445068, 7.445068, 4.113770, 2.825928, 2.176514, 1.760254, 1.502686, 1.293945, 1.131592,
336 1.008301, 0.906982, 0.823975, 0.756836, 0.698242, 0.649414, 0.606689}}}};
337 }
338 break;
339 }
340 case DG_TIMEINTEGRATION_TOULORGEC_3_7: {
341 // Using CFL factors used for CARPENTER 4/5
342 if(m_dgIntegrationMethod == "DG_INTEGRATE_GAUSS") {
343 // Legendre-Gauss nodes
344 m_cflFactor = {{{{2.497559, 2.390137, 1.978760, 1.652832, 1.413574, 1.230469, 1.097168, 0.984375, 0.890625,
345 0.812988, 0.747070, 0.691406, 0.643066, 0.602051, 0.565430, 0.533203}},
346 {{2.976074, 2.574463, 2.055664, 1.697998, 1.444092, 1.251221, 1.123047, 1.003418, 0.904541,
347 0.827637, 0.758057, 0.700684, 0.650635, 0.607910, 0.570068, 0.537109}}}};
348 } else {
349 // Legendre-Gauss-Lobatto nodes
350 m_cflFactor = {{{{5.916748, 5.916748, 3.868408, 2.753906, 2.116699, 1.722412, 1.458984, 1.259766, 1.107422,
351 0.990234, 0.895020, 0.815918, 0.750000, 0.692871, 0.644531, 0.602051}},
352 {{7.445068, 7.445068, 4.113770, 2.825928, 2.176514, 1.760254, 1.502686, 1.293945, 1.131592,
353 1.008301, 0.906982, 0.823975, 0.756836, 0.698242, 0.649414, 0.606689}}}};
354 }
355 break;
356 }
357 case DG_TIMEINTEGRATION_TOULORGEF_4_8: {
358 // Using CFL factors used for CARPENTER 4/5
359 if(m_dgIntegrationMethod == "DG_INTEGRATE_GAUSS") {
360 // Legendre-Gauss nodes
361 m_cflFactor = {{{{2.497559, 2.390137, 1.978760, 1.652832, 1.413574, 1.230469, 1.097168, 0.984375, 0.890625,
362 0.812988, 0.747070, 0.691406, 0.643066, 0.602051, 0.565430, 0.533203}},
363 {{2.976074, 2.574463, 2.055664, 1.697998, 1.444092, 1.251221, 1.123047, 1.003418, 0.904541,
364 0.827637, 0.758057, 0.700684, 0.650635, 0.607910, 0.570068, 0.537109}}}};
365 } else {
366 // Legendre-Gauss-Lobatto nodes
367 m_cflFactor = {{{{5.916748, 5.916748, 3.868408, 2.753906, 2.116699, 1.722412, 1.458984, 1.259766, 1.107422,
368 0.990234, 0.895020, 0.815918, 0.750000, 0.692871, 0.644531, 0.602051}},
369 {{7.445068, 7.445068, 4.113770, 2.825928, 2.176514, 1.760254, 1.502686, 1.293945, 1.131592,
370 1.008301, 0.906982, 0.823975, 0.756836, 0.698242, 0.649414, 0.606689}}}};
371 }
372 break;
373 }
374 }
375}
376
392template <MInt nDim>
393void DgSysEqnAcousticPerturb<nDim>::calcInitialCondition(const MFloat t,
394 const MFloat* x,
395 MFloat* const nodeVars,
396 MFloat* u) const {
401 switch(this->m_initialCondition) {
402 case 0:
403 case 1:
404 {
405 const MFloat c = (this->m_initialCondition == 0) ? 0.0 : F1B2;
406 std::fill_n(u, s_noVariables, c);
407
408 // Initialize mean velocities, mean density, mean speed of sound and its
409 // derivatives
410 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
411 nodeVars[CV::RHO0] = m_meanDensity;
412 nodeVars[CV::C0] = m_meanSpeedOfSound;
413 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
414 break;
415 }
416
417 case 2: {
419 std::fill_n(u, s_noVariables, 0.0);
420
421 // freestream mean velocity in x-direction
422 const MFloat freestreamMa = 0.3;
423 const MFloat yPos = x[1];
424 MFloat u0 = 0.0;
425 // boundary layer profile from y=0 to y=1
426 // u0 = 0.3*(2*y - 2*y^3 + y^4)
427 if(yPos < 1.0 && yPos > 0.0) {
428 u0 = freestreamMa * (2 * yPos - 2 * pow(yPos, 3) + pow(yPos, 4));
429 } else if(yPos >= 1.0) {
430 u0 = freestreamMa;
431 }
432 nodeVars[CV::UU0[0]] = u0;
433 nodeVars[CV::UU0[1]] = 0.0;
434
435 // Mean density, mean speed of sound and its derivatives
436 nodeVars[CV::RHO0] = 1.0;
437 nodeVars[CV::C0] = 1.0;
438 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
439
440 break;
441 }
442
443 case 3: {
445 std::fill_n(u, s_noVariables, 0.0);
446
447 nodeVars[CV::UU0[0]] = m_meanVelocity[0];
448 nodeVars[CV::UU0[1]] = m_meanVelocity[1];
449
450 // Mean density, mean speed of sound and its derivatives
451 nodeVars[CV::RHO0] = 1.0;
452 nodeVars[CV::C0] = 1.0;
453 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
454
455 break;
456 }
457
458 case 5: {
465
466 // Set domain-specific values and user-defined properties
467 const MFloat frequency = this->m_initialNumberWaves; // = rate of change / oscillations
468 const MFloat A = F1B5; // = max. amplitude of the i.c.
469 const MFloat length = F2; // = charactersitic length
470 const MFloat a = F1B2; // = advection velocity
471 const MFloat c = F2; // = constant offset
472
473 // Calculate initialization value
474 const MFloat omega = F2 * PI * frequency / length;
475 const MFloat init = c + A * sin(omega * std::accumulate(x, x + nDim, -a * t));
476
477 // Intialize state variables
478 for(MInt i = 0; i < nDim; i++) {
479 u[CV::UU[i]] = init;
480 }
481 u[CV::P] = init * init;
482
483 // Initialize mean velocities, mean density, mean speed of sound and its
484 // derivatives
485 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
486 nodeVars[CV::RHO0] = m_meanDensity;
487 nodeVars[CV::C0] = m_meanSpeedOfSound;
488 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
489
490 break;
491 }
492
493 case 22: {
496 MFloat mach = 0.0;
497 // Note: mean velocity is used here as free stream mach number M=u/c_inf, not u/c_0
498 for(MInt i = 0; i < nDim; i++) {
499 mach += POW2(m_meanVelocity[i]);
500 }
501 mach = sqrt(mach);
502 const MFloat kappa = 1.4;
503
504 // Freestream values non-dim with stagnation state
505 const MFloat T_inf = 1.0 / (1.0 + (kappa - 1.0) * 0.5 * mach * mach);
506 const MFloat c_inf = sqrt(T_inf);
507 const MFloat rho_inf = pow(T_inf, 1.0 / (kappa - 1.0));
508
509 std::fill_n(u, s_noVariables, 0.0); // zero IC
510
511 nodeVars[CV::RHO0] = rho_inf;
512 nodeVars[CV::C0] = c_inf;
513 nodeVars[CV::UU0[0]] = m_meanVelocity[0] * c_inf;
514 nodeVars[CV::UU0[1]] = m_meanVelocity[1] * c_inf;
515 IF_CONSTEXPR(nDim == 3) { nodeVars[CV::UU0[2]] = m_meanVelocity[2] * c_inf; }
516 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
517
518 break;
519 }
520
521 case 100: {
525
526 u[CV::U] = 0.0;
527 u[CV::V] = 0.0;
528 // u[ CV::U ] = 0.04 * x[1] * exp( -0.04 * log(2.0) *
529 // (x[0] * x[0] + x[1]* x[1]) );
530 // u[ CV::V ] = -0.04 * ( x[0] - 67.0) *
531 // exp( -0.04 * log(2.0) * (x[0] * x[0] + x[1] * x[1] ) );
532
533 u[CV::P] = exp(-log(2.0) / 9.0 * (x[0] * x[0] + x[1] * x[1]));
534
535 // Initialize mean velocities, mean density, mean speed of sound and its
536 // derivatives
537 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
538 nodeVars[CV::RHO0] = m_meanDensity;
539 nodeVars[CV::C0] = m_meanSpeedOfSound;
540 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
541 break;
542 }
543
544 case 101: {
548
549 u[CV::U] = 0.0;
550 u[CV::V] = 0.0;
551 u[CV::P] = 2.0 - exp(-1.0 / 0.25 * (x[0] * x[0] + x[1] * x[1]));
552
553 // Initialize mean velocities, mean density, mean speed of sound and its
554 // derivatives
555 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
556 nodeVars[CV::RHO0] = m_meanDensity;
557 nodeVars[CV::C0] = m_meanSpeedOfSound;
558 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
559 break;
560 }
561
562 case 102: {
566
567 u[CV::U] = 0.0;
568 u[CV::V] = 0.0;
569 u[CV::P] = 2.0 - exp(-1.0 / 0.25 * ((x[0] + 7.5) * (x[0] + 7.5) + x[1] * x[1]));
570
571 // Initialize mean velocities, mean density, mean speed of sound and its
572 // derivatives
573 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
574 nodeVars[CV::RHO0] = m_meanDensity;
575 nodeVars[CV::C0] = m_meanSpeedOfSound;
576 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
577 break;
578 }
579
580 case 105: {
584
585 // Reset velocities
586 for(MInt i = 0; i < nDim; i++) {
587 u[CV::UU[i]] = 0.0;
588 }
589 // u[ CV::U ] = 0.04 * x[1] * exp( -0.04 * log(2.0) *
590 // (x[0] * x[0] + x[1]* x[1]) );
591 // u[ CV::V ] = -0.04 * ( x[0] - 67.0) *
592 // exp( -0.04 * log(2.0) * (x[0] * x[0] + x[1] * x[1] ) );
593
594 u[CV::P] = 2.0 - exp(-1.0 / 2.0 * (x[0] * x[0] + x[1] * x[1] + x[2] * x[2]));
595
596 // Initialize mean velocities, mean density, mean speed of sound and its
597 // derivatives
598 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
599 nodeVars[CV::RHO0] = m_meanDensity;
600 nodeVars[CV::C0] = m_meanSpeedOfSound;
601 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
602 break;
603 }
604
605 case 106: {
609
610 // Reset velocities
611 for(MInt i = 0; i < nDim; i++) {
612 u[CV::UU[i]] = 0.0;
613 }
614
615 u[CV::P] = exp(-1.0 / 2.0 * (x[0] * x[0] + x[1] * x[1] + x[2] * x[2]));
616
617 // Initialize mean velocities, mean density, mean speed of sound and its
618 // derivatives
619 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
620 nodeVars[CV::RHO0] = m_meanDensity;
621 nodeVars[CV::C0] = m_meanSpeedOfSound;
622 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
623 break;
624 }
625
626 case 107: {
630
631 // Reset velocities
632 for(MInt i = 0; i < nDim; i++) {
633 u[CV::UU[i]] = 0.0;
634 }
635
636 u[CV::P] = 2.0 - exp(-1.0 / 2.0 * ((x[0] + 7.5) * (x[0] + 7.5) + x[1] * x[1] + x[2] * x[2]));
637
638 // Initialize mean velocities, mean density, mean speed of sound and its
639 // derivatives
640 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
641 nodeVars[CV::RHO0] = m_meanDensity;
642 nodeVars[CV::C0] = m_meanSpeedOfSound;
643 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
644 break;
645 }
646
647 case 115: {
652
653 // Reset velocities
654 for(MInt i = 0; i < nDim; i++) {
655 u[CV::UU[i]] = 0.0;
656 }
657
658 // Set pressure distribution
659 const MFloat r = x[0] * x[0] + x[1] * x[1] + (x[2] + 45.0) * (x[2] + 45.0);
660 u[CV::P] = exp(1 + 0.28 * exp(-r / (2.5))) - exp(1.0);
661
662 // Initialize mean velocities, mean density, mean speed of sound and its
663 // derivatives
664 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
665 nodeVars[CV::RHO0] = m_meanDensity;
666 nodeVars[CV::C0] = m_meanSpeedOfSound;
667 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
668 break;
669 }
670
671 case 201: {
675
676 // Reset velocities
677 for(MInt i = 0; i < nDim; i++) {
678 u[CV::UU[i]] = 0.0;
679 }
680
681 // Set pressure distribution
682 u[CV::P] = exp(-1.0 * log(2.0) / 25 * (x[0] * x[0] + (x[1] - 25.0) * (x[1] - 25.0)));
683
684 // Initialize mean velocities, mean density, mean speed of sound and its
685 // derivatives
686 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
687 nodeVars[CV::RHO0] = m_meanDensity;
688 nodeVars[CV::C0] = m_meanSpeedOfSound;
689 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
690 break;
691 }
692
693 case 202: {
698
699 u[CV::U] = sin(PI * x[0]);
700 u[CV::V] = 0.0;
701 u[CV::P] = sin(PI * x[0]);
702
703 // Initialize mean velocities, mean density, mean speed of sound and its
704 // derivatives
705 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
706 nodeVars[CV::RHO0] = m_meanDensity;
707 nodeVars[CV::C0] = m_meanSpeedOfSound;
708 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
709 break;
710 }
711
712
713 case 790: {
717
718 // Initial condition for setting everything to zero
719 // Averaged quantities are read from the given mean file
720 std::fill_n(u, s_noVariables, 0.0);
721
722 // Set constant speed of sound and reset its derivatives
723 if(m_constantSpeedOfSound) {
724 nodeVars[CV::C0] = m_meanSpeedOfSound;
725 for(MInt i = 0; i < nDim; i++) {
726 nodeVars[CV::DC0[i]] = 0.0;
727 }
728 }
729 break;
730 }
731
732 case 1661: {
736 std::fill_n(u, s_noVariables, 0.0); // zero initial condition
737
738 // Set default mean velocitites everywhere
739 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
740
741 // Initialize mean velocities, mean density, mean speed of sound and its derivatives inside a
742 // box, use zero/rho0/c0 values outside
743 MBool pointInBox = true;
744 for(MInt dim = 0; dim < nDim; dim++) {
745 if(fabs(x[dim]) > 0.75) {
746 pointInBox = false;
747 }
748 }
749 if(pointInBox) {
750 const MFloat uj = 0.9 * m_meanSpeedOfSound;
751 const MFloat r0 = 0.1;
752 const MFloat delOmega = 0.05 * r0;
753 MFloat rTemp = 0.0;
754 for(MInt dim = 1; dim < nDim; dim++) {
755 rTemp += x[dim] * x[dim];
756 }
757 const MFloat r = sqrt(rTemp);
758 nodeVars[CV::U0] = uj * (0.5 + 0.5 * tanh((r0 - r) / (2.0 * delOmega)));
759 nodeVars[CV::V0] = 0.3 * (1.0 + 0.1 * x[1]);
760 nodeVars[CV::RHO0] = m_meanDensity * (1.0 + 0.1 * x[0]);
761 nodeVars[CV::C0] = m_meanSpeedOfSound * (1.0 + 0.1 * x[1]);
762
763 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.1);
764 nodeVars[CV::DC0[0]] = 0.01;
765 } else {
766 std::fill_n(&nodeVars[CV::UU0[0]], nDim, 0.0);
767 nodeVars[CV::RHO0] = m_meanDensity;
768 nodeVars[CV::C0] = m_meanSpeedOfSound;
769 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
770 }
771 break;
772 }
773
774 case 1662: {
778 std::fill_n(u, s_noVariables, 0.0); // zero initial condition
779
780 const MFloat uj = 0.9 * m_meanSpeedOfSound;
781 const MFloat r0 = 0.1;
782 const MFloat delOmega = 0.05 * r0;
783 MFloat rTemp = 0.0;
784 for(MInt dim = 1; dim < nDim; dim++) {
785 rTemp += x[dim] * x[dim];
786 }
787 const MFloat r = sqrt(rTemp);
788
789 // Set default node variables
790 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
791 nodeVars[CV::RHO0] = m_meanDensity;
792 nodeVars[CV::C0] = m_meanSpeedOfSound;
793 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
794
795 // Initialize mean velocities, mean density, mean speed of sound and its derivatives for
796 // x > 0.5
797 if(x[0] >= 0.5) {
798 nodeVars[CV::U0] = uj * (0.5 + 0.5 * tanh((r0 - r) / (2.0 * delOmega)));
799 nodeVars[CV::C0] = m_meanSpeedOfSound * (1.0 + 0.1 * x[1]);
800 nodeVars[CV::DC0[0]] = 0.01;
801 // Set density only for y > 0.5
802 if(x[1] >= 0.5) {
803 nodeVars[CV::RHO0] = m_meanDensity * (1.0 + 0.1 * x[0]);
804 }
805 }
806
807 break;
808 }
809 default:
810 mTerm(1, AT_, "The specified initial condition (" + std::to_string(this->m_initialCondition) + ") is not valid!");
811 }
812}
813
814
826template <MInt nDim>
827void DgSysEqnAcousticPerturb<nDim>::calcFlux(const MFloat* const nodeVars,
828 const MFloat* const q,
829 const MInt noNodes1D,
830 MFloat* const flux) const {
831 // TRACE();
832
833 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
834 const MFloatTensor U(const_cast<MFloat*>(q), noNodes1D, noNodes1D, noNodes1D3, s_noVariables);
835 const MFloatTensor U0(const_cast<MFloat*>(nodeVars), noNodes1D, noNodes1D, noNodes1D3, s_noNodeVars);
836 MFloatTensor f(flux, noNodes1D, noNodes1D, noNodes1D3, nDim, s_noVariables);
837
838 // The following flux equations are calculated here (if 2D, consider all
839 // z-components to be zero; for an incompressible case: cBar=1, rhoBar=1):
840
841 // Flux in x-direction:
842 // f_x(CV[u]) = u*uBar + v*vBar + w*wBar + p/rhoBar
843 // f_x(CV[v]) = 0
844 // f_x(CV[w]) = 0
845 // f_x(CV[p]) = rhoBar*cBar^2*u + uBar*p
846
847 // Flux in y-direction:
848 // f_y(CV[u]) = 0
849 // f_y(CV[v]) = u*uBar + v*vBar + w*wBar + p/rhoBar
850 // f_y(CV[w]) = 0
851 // f_y(CV[p]) = rhoBar*cBar^2*v + vBar*p
852
853 // Flux in z-direction:
854 // f_z(CV[u]) = 0
855 // f_z(CV[v]) = 0
856 // f_z(CV[w]) = u*uBar + v*vBar + w*wBar + p/rhoBar
857 // f_z(CV[p]) = rhoBar*cBar^2*w + wBar*p
858
859 for(MInt i = 0; i < noNodes1D; i++) {
860 for(MInt j = 0; j < noNodes1D; j++) {
861 for(MInt k = 0; k < noNodes1D3; k++) {
862 for(MInt d = 0; d < nDim; d++) {
863 // Momentum fluxes
864 for(MInt fd = 0; fd < nDim; fd++) {
865 f(i, j, k, d, CV::UU[fd]) = 0.0;
866 }
867 f(i, j, k, d, CV::UU[d]) = std::inner_product(&U0(i, j, k, CV::UU0[0]),
868 &U0(i, j, k, CV::UU0[0]) + nDim,
869 &U(i, j, k, CV::UU[0]),
870 U(i, j, k, CV::P) / U0(i, j, k, CV::RHO0));
871
872 // Energy flux
873 f(i, j, k, d, CV::P) = U0(i, j, k, CV::RHO0) * POW2(U0(i, j, k, CV::C0)) * U(i, j, k, CV::UU[d])
874 + U0(i, j, k, CV::UU0[d]) * U(i, j, k, CV::P);
875 }
876 }
877 }
878 }
879}
880
881
895template <MInt nDim>
896void DgSysEqnAcousticPerturb<nDim>::calcFlux1D(const MFloat* const nodeVars,
897 const MFloat* const q,
898 const MInt noNodes1D,
899 const MInt dirId,
900 MFloat* const flux) const {
901 // TRACE();
902
903 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
904 const MFloatTensor U(const_cast<MFloat*>(q), noNodes1D, noNodes1D3, s_noVariables);
905 const MFloatTensor U0(const_cast<MFloat*>(nodeVars), noNodes1D, noNodes1D3, s_noNodeVars);
906 MFloatTensor f(flux, noNodes1D, noNodes1D3, s_noVariables);
907
908 for(MInt i = 0; i < noNodes1D; i++) {
909 for(MInt j = 0; j < noNodes1D3; j++) {
910 // Momentum fluxes
911 for(MInt d = 0; d < nDim; d++) {
912 if(d == dirId) {
913 f(i, j, CV::UU[d]) = std::inner_product(&U0(i, j, CV::UU0[0]),
914 &U0(i, j, CV::UU0[0]) + nDim,
915 &U(i, j, CV::UU[0]),
916 U(i, j, CV::P) / U0(i, j, CV::RHO0));
917 } else {
918 f(i, j, CV::UU[d]) = 0.0;
919 }
920 }
921
922 // Energy flux
923 f(i, j, CV::P) = U0(i, j, CV::RHO0) * POW2(U0(i, j, CV::C0)) * U(i, j, CV::UU[dirId])
924 + U0(i, j, CV::UU0[dirId]) * U(i, j, CV::P);
925 }
926 }
927}
928
929
941template <MInt nDim>
942MFloat DgSysEqnAcousticPerturb<nDim>::cflFactor(const MInt polyDeg) const {
943 if(polyDeg > s_maxPolyDeg) {
944 TERMM(1, "Polynomial degree exceeds maximum supported");
945 }
946 // Note: The factor "0.95" is used to avoid numerical instabilities when using the full
947 // cfl factor value.
948
949 const MFloat factor = 0.95 * m_cflFactor[nDim - 2][polyDeg];
950
951 return factor;
952}
953
954
968template <MInt nDim>
969void DgSysEqnAcousticPerturb<nDim>::calcSource(const MFloat* const nodeVars, const MFloat* const u,
970 const MInt noNodes1D, const MFloat t, const MFloat* const x,
971 MFloat* const src) const {
972 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
973 MFloatTensor s(src, noNodes1D, noNodes1D, noNodes1D3, s_noVariables);
974 const MFloatTensor X(const_cast<MFloat*>(x), noNodes1D, noNodes1D, noNodes1D3, nDim);
975 const MFloatTensor U(const_cast<MFloat*>(u), noNodes1D, noNodes1D, noNodes1D3, s_noVariables);
976 const MFloatTensor U0(const_cast<MFloat*>(nodeVars), noNodes1D, noNodes1D, noNodes1D3, s_noNodeVars);
977
981 switch(this->m_sourceTerm) {
982 case 0: {
984 s.set(F0);
985 break;
986 }
987
988 case 1: {
990 mTerm(1, AT_, "The APE-1 sources are not yet implemented!");
991 break;
992 }
993
994 case 2: {
996 mTerm(1, AT_, "The APE-2 sources are not yet implemented!");
997 break;
998 }
999
1000 case 3: {
1002 mTerm(1, AT_, "The APE-3 sources are not yet implemented!");
1003 break;
1004 }
1005
1006 case 4: {
1008 mTerm(1, AT_, "The APE-4 sources are not yet implemented!");
1009 break;
1010 }
1011
1012 case 5: {
1025
1026 // Set domain-specific values and user-defined properties
1027 const MFloat frequency = this->m_initialNumberWaves; // = rate of change / oscillations
1028 const MFloat A = F1B5; // = max. amplitude of the i.c.
1029 const MFloat length = F2; // = charactersitic length
1030 const MFloat a = F1B2; // = advection velocity
1031 const MFloat c = F2; // = constant offset
1032
1033 // Calculate initialization value
1034 const MFloat omega = F2 * PI * frequency / length;
1035
1036 for(MInt i = 0; i < noNodes1D; i++) {
1037 for(MInt j = 0; j < noNodes1D; j++) {
1038 for(MInt k = 0; k < noNodes1D3; k++) {
1039 const MFloat sumX = std::accumulate(&X(i, j, k, 0), &X(i, j, k, 0) + nDim, F0);
1040 const MFloat inner = omega * (sumX - a * t);
1041 const MFloat sumVelBarMa = std::accumulate(&U0(i, j, k, CV::UU0[0]), &U0(i, j, k, CV::UU0[0]) + nDim, -a);
1042 const MFloat F2sinA = F2 * sin(inner) * A;
1043 const MFloat cosA = cos(inner) * A;
1044
1045 // Set velocity sources
1046 for(MInt n = 0; n < nDim; n++) {
1047 s(i, j, k, CV::UU[n]) = omega * cosA * (F2 * c + sumVelBarMa + F2sinA);
1048 }
1049 // Set pressure source
1050 s(i, j, k, CV::P) = omega * cosA * (nDim + F2 * c * sumVelBarMa + sumVelBarMa * F2sinA);
1051 }
1052 }
1053 }
1054 break;
1055 }
1056
1057 case 6:
1061 {
1062 const MFloat pos[3] = {0.0, 0.0, 0.0}; // Source center
1063 const MFloat r = 0.1; // Radius
1064 const MFloat A = 1.0; // Amplitude
1065 const MFloat f = this->m_initialNumberWaves; // Frequency
1066 // Calculate constant factors
1067 const MFloat factor = A * sin(2 * PI * f * t);
1068 const MFloat expFactor = -1.0 / (2 * r * r);
1069
1070 for(MInt i = 0; i < noNodes1D; i++) {
1071 for(MInt j = 0; j < noNodes1D; j++) {
1072 for(MInt k = 0; k < noNodes1D3; k++) {
1073 // Set velocity sources to zero
1074 std::fill_n(&s(i, j, k, CV::UU[0]), nDim, F0);
1075
1076 // Compute squared distance between point and source center
1077 MFloat distSquared = F0;
1078 for(MInt n = 0; n < nDim; n++) {
1079 distSquared += POW2(X(i, j, k, n) - pos[n]);
1080 }
1081
1082 // Set pressure source
1083 s(i, j, k, CV::P) = factor * exp(expFactor * distSquared);
1084 }
1085 }
1086 }
1087 break;
1088 }
1089
1090 case 60:
1094 {
1095 const MFloat pos[3] = {25.0, 0.0, 0.0};
1096 const MFloat r = 0.25; // Radius
1097 const MFloat A = 0.1; // Amplitude
1098 const MFloat f = 0.5; // Frequency
1099 // Calculate constant factors
1100 const MFloat factor = A * sin(2 * PI * f * t);
1101 const MFloat expFactor = -1.0 / (2 * r * r);
1102
1103 for(MInt i = 0; i < noNodes1D; i++) {
1104 for(MInt j = 0; j < noNodes1D; j++) {
1105 for(MInt k = 0; k < noNodes1D3; k++) {
1106 // Set velocity sources to zero
1107 std::fill_n(&s(i, j, k, CV::UU[0]), nDim, F0);
1108
1109 // Compute squared distance between point and source center
1110 MFloat distSquared = F0;
1111 for(MInt n = 0; n < nDim; n++) {
1112 distSquared += POW2(X(i, j, k, n) - pos[n]);
1113 }
1114
1115 // Set pressure source
1116 s(i, j, k, CV::P) = factor * exp(expFactor * distSquared);
1117 }
1118 }
1119 }
1120 break;
1121 }
1122
1123 case 7: {
1124 TERMM(1, "dipole source needs testing");
1130 const MFloat d = 0.01; // Separation distance of the two monopole sources
1131 const MFloat pos1[3] = {-0.5 * d, 0.0, 0.0}; // Source center monopole #1
1132 const MFloat pos2[3] = {0.5 * d, 0.0, 0.0}; // Source center monopole #2
1133
1134 const MFloat r = 0.1; // Radius
1135 const MFloat A = 1.0; // Amplitude
1136 const MFloat f = 1.0; // Frequency
1137
1138 // Calculate constant factors
1139 const MFloat factor1 = A * sin(2 * PI * f * t);
1140 const MFloat factor2 = A * sin(2 * PI * (f * t + 0.5)); // Opposite phase for monopole #2
1141
1142 const MFloat expFactor = -1.0 / (2 * r * r);
1143
1144 for(MInt i = 0; i < noNodes1D; i++) {
1145 for(MInt j = 0; j < noNodes1D; j++) {
1146 for(MInt k = 0; k < noNodes1D3; k++) {
1147 // Set velocity sources to zero
1148 std::fill_n(&s(i, j, k, CV::UU[0]), nDim, F0);
1149
1150 // Compute squared distances between point and source centers
1151 MFloat distSquared1 = F0;
1152 MFloat distSquared2 = F0;
1153 for(MInt n = 0; n < nDim; n++) {
1154 distSquared1 += POW2(X(i, j, k, n) - pos1[n]);
1155 distSquared2 += POW2(X(i, j, k, n) - pos2[n]);
1156 }
1157
1158 // Set pressure source
1159 s(i, j, k, CV::P) = factor1 * exp(expFactor * distSquared1) + factor2 * exp(expFactor * distSquared2);
1160 }
1161 }
1162 }
1163 break;
1164 }
1165
1166 case 70: {
1168 const MFloat freq = 20.0; // 1.0/30.0;
1169 const MFloat omega = 2.0 * PI * freq;
1170 const MFloat eps = 0.5;
1171 const MFloat alpha = std::log(2.0) / 2.0;
1172 const MFloat x_s = 125.0;
1173 const MFloat y_s = 0.0;
1174
1175 for(MInt i = 0; i < noNodes1D; i++) {
1176 for(MInt j = 0; j < noNodes1D; j++) {
1177 for(MInt k = 0; k < noNodes1D3; k++) {
1178 // Calculate source term
1179 const MFloat x_norm = X(i, j, k, 0) - x_s;
1180 const MFloat y_norm = X(i, j, k, 1) - y_s;
1181 const MFloat f = eps * std::exp(-alpha * (x_norm * x_norm + y_norm * y_norm));
1182
1183 // Apply source term
1184 const MFloat source = f * std::sin(omega * t);
1185 s(i, j, k, CV::U) = source;
1186 s(i, j, k, CV::V) = source;
1187 s(i, j, k, CV::P) = 0.0;
1188 }
1189 }
1190 }
1191 break;
1192 }
1193
1194 case 71: {
1196 const MFloat freq = 1.0 / 5.0; // 20.0;//1.0/30.0;
1197 const MFloat omega = 2.0 * PI * freq;
1198 const MFloat eps = 1e-4; // 0.5;
1199 const MFloat alpha = std::log(2.0) / 2.0;
1200 const MFloat x_s = 125.0;
1201 const MFloat y_s = 0.0;
1202
1203 for(MInt i = 0; i < noNodes1D; i++) {
1204 for(MInt j = 0; j < noNodes1D; j++) {
1205 for(MInt k = 0; k < noNodes1D3; k++) {
1206 // Calculate source term
1207 const MFloat x_norm = X(i, j, k, 0) - x_s;
1208 const MFloat y_norm = X(i, j, k, 1) - y_s;
1209 const MFloat f = eps * std::exp(-alpha * (x_norm * x_norm + y_norm * y_norm));
1210
1211 // Apply source term
1212 const MFloat source = f * std::sin(omega * t);
1213 s(i, j, k, CV::U) = source;
1214 s(i, j, k, CV::V) = 0.0;
1215 s(i, j, k, CV::P) = 0.0;
1216 }
1217 }
1218 }
1219 break;
1220 }
1221
1222 case 800: {
1223 TERMM_IF_COND(nDim == 3, "only useful in 2D");
1229
1230 // Determine position-independent quantities
1231 const MFloat gamma = 1.0;
1232 const MFloat sigma = 1.0;
1233 const MFloat omega = gamma / 4. / PI;
1234 const MFloat alpha = gamma * gamma / 8 / PI / PI / sigma / sigma;
1235 const MFloat theta = omega * t;
1236 const MFloat bx = std::cos(theta);
1237 const MFloat by = std::sin(theta);
1238
1251 // Read filter properties
1252 const MFloat filterSlopeWidth = Context::getSolverProperty<MFloat>("filterSlopeWidth", this->m_solverId, AT_);
1253
1254 std::array<MFloat, nDim> filterRegionMin{};
1255 std::array<MFloat, nDim> filterRegionMax{};
1256 for(MInt i = 0; i < 2; i++) {
1257 filterRegionMin[i] = Context::getSolverProperty<MFloat>("filterRegionMin", this->m_solverId, AT_, i);
1258 filterRegionMax[i] = Context::getSolverProperty<MFloat>("filterRegionMax", this->m_solverId, AT_, i);
1259 }
1260
1261 // Determine position-dependent quantities and apply sources
1262 for(MInt i = 0; i < noNodes1D; i++) {
1263 for(MInt j = 0; j < noNodes1D; j++) {
1264 // Calculate source term
1265 const MFloat rxPos = X(i, j, 0, 0) - bx;
1266 const MFloat ryPos = X(i, j, 0, 1) - by;
1267 const MFloat rPos = std::sqrt(rxPos * rxPos + ryPos * ryPos);
1268 const MFloat rxNeg = X(i, j, 0, 0) + bx;
1269 const MFloat ryNeg = X(i, j, 0, 1) + by;
1270 const MFloat rNeg = std::sqrt(rxNeg * rxNeg + ryNeg * ryNeg);
1271 const MFloat qmEwert =
1272 alpha * (std::exp(-rPos * rPos / 2. / sigma / sigma) - std::exp(-rNeg * rNeg / 2. / sigma / sigma));
1273 const MFloat qmEwertX = qmEwert * std::cos(theta);
1274 const MFloat qmEwertY = qmEwert * std::sin(theta);
1275
1276 // Determine filter
1277 const MFloat filter = maia::filter::slope::cosbox<nDim>(&filterRegionMin[0], &filterRegionMax[0],
1278 filterSlopeWidth, &X(i, j, 0, 0));
1279
1280 // Apply source term
1281 s(i, j, 0, CV::U) = filter * qmEwertX;
1282 s(i, j, 0, CV::V) = filter * qmEwertY;
1283 s(i, j, 0, CV::P) = 0.0;
1284 }
1285 }
1286 break;
1287 }
1288
1289 default:
1290 mTerm(1, AT_, "The specified source term is not valid!");
1291 }
1292
1293 // Return if compresible source term should not be considered
1294 if(!m_compressibleSourceTerm) {
1295 return;
1296 }
1297
1298 // Add compressible source term to p'
1299 // q_cons = 2 * (rhoBar * cBar * u' + uBar * 1/cBar * p') * grad(cBar)
1300 for(MInt i = 0; i < noNodes1D; i++) {
1301 for(MInt j = 0; j < noNodes1D; j++) {
1302 for(MInt k = 0; k < noNodes1D3; k++) {
1303 for(MInt dim = 0; dim < nDim; dim++) {
1304 s(i, j, k, CV::P) += 2
1305 * (U0(i, j, k, CV::RHO0) * U0(i, j, k, CV::C0) * U(i, j, k, CV::UU[dim])
1306 + U0(i, j, k, CV::UU0[dim]) / U0(i, j, k, CV::C0) * U(i, j, k, CV::P))
1307 * U0(i, j, k, CV::DC0[dim]);
1308 }
1309 }
1310 }
1311 }
1312}
1313
1327template <MInt nDim>
1328void DgSysEqnAcousticPerturb<nDim>::calcSpongeSource(const MFloat* NotUsed(nodeVars), const MFloat* u,
1329 const MInt noNodes1D, const MFloat* eta, MFloat* src) const {
1330 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
1331
1332 // gamma = heat capacity ratio for monoatomic gas
1333 MFloat gammaMinusOne = 5.0 / 3.0 - 1.0;
1334 MFloat FgammaMinusOne = 1.0 / gammaMinusOne;
1335
1336 const MFloatTensor spongeEta(const_cast<MFloat*>(eta), noNodes1D, noNodes1D, noNodes1D3);
1337 MFloatTensor uState(const_cast<MFloat*>(u), noNodes1D, noNodes1D, noNodes1D3, s_noVariables);
1338 MFloatTensor spongeSource(src, noNodes1D, noNodes1D, noNodes1D3, s_noVariables);
1339
1340 for(MInt i = 0; i < noNodes1D; i++) {
1341 for(MInt j = 0; j < noNodes1D; j++) {
1342 for(MInt k = 0; k < noNodes1D3; k++) {
1343 MFloat deltaP = (m_spongePressureInfy - uState(i, j, k, CV::P)) * FgammaMinusOne;
1344
1345 spongeSource(i, j, k, CV::P) = m_spongeSigma * spongeEta(i, j, k) * deltaP;
1346 for(MInt l = 0; l < nDim; l++) {
1347 spongeSource(i, j, k, CV::UU[l]) = 0.0;
1348 }
1349 }
1350 }
1351 }
1352}
1353
1354
1367template <MInt nDim>
1368MFloat DgSysEqnAcousticPerturb<nDim>::getTimeStep(const MFloat* const nodeVars,
1369 const MFloat* const NotUsed(u),
1370 const MInt noNodes1D,
1371 const MFloat invJacobian,
1372 const MInt sbpMode) const {
1373 // TRACE();
1374
1375 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
1376 const MFloat inf = std::numeric_limits<MFloat>::infinity();
1377 const MFloatTensor U0(const_cast<MFloat*>(nodeVars), noNodes1D, noNodes1D, noNodes1D3, s_noNodeVars);
1378
1379 MFloat maxLambda[] = {-inf, -inf, -inf};
1380 for(MInt i = 0; i < noNodes1D; i++) {
1381 for(MInt j = 0; j < noNodes1D; j++) {
1382 for(MInt k = 0; k < noNodes1D3; k++) {
1383 for(MInt d = 0; d < nDim; d++) {
1384 maxLambda[d] = std::max(maxLambda[d], fabs(U0(i, j, k, CV::UU0[d])) + U0(i, j, k, CV::C0));
1385 }
1386 }
1387 }
1388 }
1389
1390 MFloat dt;
1391 if(sbpMode) {
1392 dt = this->cfl() * F2 / (invJacobian * std::accumulate(&maxLambda[0], &maxLambda[0] + nDim, F0) * (noNodes1D - 1));
1393 } else {
1394 const MInt polyDeg = noNodes1D - 1;
1395 dt = this->cflScaled(polyDeg) * cflFactor(polyDeg) * F2
1396 / (invJacobian * std::accumulate(&maxLambda[0], &maxLambda[0] + nDim, F0));
1397 }
1398 return dt;
1399}
1400
1401
1422template <MInt nDim>
1423void DgSysEqnAcousticPerturb<nDim>::calcRiemann(const MFloat* nodeVarsL,
1424 const MFloat* nodeVarsR,
1425 const MFloat* stateL,
1426 const MFloat* stateR,
1427 const MInt noNodes1D,
1428 const MInt dirId,
1429 MFloat* flux) const {
1430 // TRACE();
1431
1432 // Solve Riemann problem
1433 switch(this->m_riemannSolver) {
1434 case 0: // local Lax-Friedrichs flux
1435 {
1436 calcRiemannLaxFriedich(nodeVarsL, nodeVarsR, stateL, stateR, noNodes1D, dirId, flux);
1437 break;
1438 }
1439
1440 case 1: // Roe's flux
1441 {
1442 TERMM(1, "Implementation does not work for internal fluxes, see comments");
1443 calcRiemannRoe(nodeVarsL, nodeVarsR, stateL, stateR, noNodes1D, dirId, flux);
1444 break;
1445 }
1446
1447 default:
1448 mTerm(1, AT_, "The specified Riemann solver is not valid!");
1449 break;
1450 }
1451}
1452
1453
1460template <MInt nDim>
1461void DgSysEqnAcousticPerturb<nDim>::calcRiemannLaxFriedich(const MFloat* nodeVarsL,
1462 const MFloat* nodeVarsR,
1463 const MFloat* stateL,
1464 const MFloat* stateR,
1465 const MInt noNodes1D,
1466 const MInt dirId,
1467 MFloat* flux) const {
1468 // TRACE();
1469
1470 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
1471 const MFloatTensor uL(const_cast<MFloat*>(stateL), noNodes1D, noNodes1D3, s_noVariables);
1472 const MFloatTensor uR(const_cast<MFloat*>(stateR), noNodes1D, noNodes1D3, s_noVariables);
1473
1474 const MFloatTensor nL(const_cast<MFloat*>(nodeVarsL), noNodes1D, noNodes1D3, s_noNodeVars);
1475 const MFloatTensor nR(const_cast<MFloat*>(nodeVarsR), noNodes1D, noNodes1D3, s_noNodeVars);
1476
1477#ifndef _OPENMP
1478 MFloatScratchSpace maxLambda(noNodes1D, noNodes1D3, AT_, "maxLambda");
1479#else
1480 MFloatTensor maxLambda(noNodes1D, noNodes1D3);
1481#endif
1482
1483 // Calculate maximum eigenvalue
1484 for(MInt i = 0; i < noNodes1D; i++) {
1485 for(MInt j = 0; j < noNodes1D3; j++) {
1486 maxLambda(i, j) = std::max(fabs(nL(i, j, CV::UU0[dirId])) + nL(i, j, CV::C0),
1487 fabs(nR(i, j, CV::UU0[dirId])) + nR(i, j, CV::C0));
1488 }
1489 }
1490
1491#ifndef _OPENMP
1492 MFloatScratchSpace fluxL(noNodes1D, noNodes1D3, s_noVariables, AT_, "fluxL");
1493 MFloatScratchSpace fluxR(noNodes1D, noNodes1D3, s_noVariables, AT_, "fluxR");
1494#else
1495 MFloatTensor fluxL(noNodes1D, noNodes1D3, s_noVariables);
1496 MFloatTensor fluxR(noNodes1D, noNodes1D3, s_noVariables);
1497#endif
1498
1499 // Calculate flux from left and right state
1500 calcFlux1D(nodeVarsL, stateL, noNodes1D, dirId, &fluxL[0]);
1501 calcFlux1D(nodeVarsR, stateR, noNodes1D, dirId, &fluxR[0]);
1502
1503 // Solve Riemann problem
1504 MFloatTensor riemann(flux, noNodes1D, noNodes1D3, s_noVariables);
1505 for(MInt i = 0; i < noNodes1D; i++) {
1506 for(MInt j = 0; j < noNodes1D3; j++) {
1507 for(MInt n = 0; n < s_noVariables; n++) {
1508 riemann(i, j, n) = 0.5 * ((fluxL(i, j, n) + fluxR(i, j, n)) - maxLambda(i, j) * (uR(i, j, n) - uL(i, j, n)));
1509 }
1510 }
1511 }
1512}
1513
1514
1515/*---2D Fluxes in Primitive form---
1516 * The fluxes computed at surface, we could derive
1517 * inward and outward waves compositions via characteristic analysis,
1518 * and enforce inward outwaves to zero.
1519 *
1520 * Ref. Bauer's paper:
1521 * "-Application of a Discontinuous Galerkin Method to Discretize
1522 * Acoustic Perturbation Equations, AIAA, May, 2011
1523 *
1524 * \author H.-J. Cheng <h-j.cheng@aia.rwth-aachen.de>
1525 * \date 2014-12-25
1526 *
1527 * labels:DG IMPORTANT: Please note that the mean flow matrix calculation is just a hack
1528 * right now and only works if this method is used for a boundary condition, not
1529 * for internal fluxes!
1530 */
1531template <>
1532inline void DgSysEqnAcousticPerturb<2>::calcRiemannRoe(const MFloat* nodeVarsL,
1533 const MFloat* nodeVarsR,
1534 const MFloat* stateL,
1535 const MFloat* stateR,
1536 const MInt noNodes1D,
1537 const MInt dirId,
1538 MFloat* flux) const {
1539 // TRACE();
1540
1541 const MInt nDim = 2;
1542 const MFloatTensor uL(const_cast<MFloat*>(stateL), noNodes1D, noVars());
1543 const MFloatTensor uR(const_cast<MFloat*>(stateR), noNodes1D, noVars());
1544 const MFloatTensor nL(const_cast<MFloat*>(nodeVarsL), noNodes1D, noNodeVars());
1545 const MFloatTensor nR(const_cast<MFloat*>(nodeVarsR), noNodes1D, noNodeVars());
1546 MFloatTensor f(flux, noNodes1D, s_noVariables);
1547
1548 // TODO labels:DG use node variables instead
1549 // Define mean flow constants
1550 const MFloat c = 1.0;
1551 const MFloat rho = 1.0;
1552
1553 // Normal direction
1554 const std::array<MInt, nDim> n = {{1 - dirId, dirId}};
1555
1556 // Characteristic boundary conditions
1557 for(MInt i = 0; i < noNodes1D; i++) {
1558 // FIXME labels:DG,totest The following is just a hack and probably does not work when the
1559 // Roe Riemann solver is used for internal fluxes
1560 // Determine mean flow
1561 const std::array<MFloat, nDim> meanVelocity = {
1562 {0.5 * (nL(i, CV::U0) + nR(i, CV::U0)), 0.5 * (nL(i, CV::V0) + nR(i, CV::V0))}};
1563
1564 // Mean flow matrix with simplified transformation matrix
1565 const std::array<MFloat, nDim> um = {
1566 {meanVelocity[0] * n[0] + meanVelocity[1] * n[1], -meanVelocity[0] * n[1] + meanVelocity[1] * n[0]}};
1567
1568 // stateL(vecL) and stateR(vecR) variables multiplied by transformation
1569 // matrix (see reference)
1570 const std::array<MFloat, nDim> vecL = {
1571 {uL(i, CV::UU[0]) * n[0] + uL(i, CV::UU[1]) * n[1], -uL(i, CV::UU[0]) * n[1] + uL(i, CV::UU[1]) * n[0]}};
1572 const std::array<MFloat, nDim> vecR = {
1573 {uR(i, CV::UU[0]) * n[0] + uR(i, CV::UU[1]) * n[1], -uR(i, CV::UU[0]) * n[1] + uR(i, CV::UU[1]) * n[0]}};
1574
1575 // Calculate characteristic waves
1576 // L1 and L2 composite from vecL and vecR
1577 const MFloat L1 = 0.5 * (uR(i, CV::P) + rho * c * (-vecR[0] + um[1] / (um[0] - c) * (-vecR[1])));
1578 const MFloat L2 = 0.5 * (uL(i, CV::P) + rho * c * (vecL[0] + um[1] / (um[0] + c) * vecL[1]));
1579
1580 // Calculate Riemann flux
1581 f(i, CV::P) = (um[0] - c) * (L1) + (um[0] + c) * (L2);
1582 f(i, CV::UU[0]) = 1 / rho / c * (-1 * ((um[0] - c) * L1) + ((um[0] + c) * L2)) * n[0];
1583 f(i, CV::UU[1]) = 1 / rho / c * (-1 * ((um[0] - c) * L1) + ((um[0] + c) * L2)) * n[1];
1584 }
1585}
1586
1587
1588/*---3D Fluxes in Primitive form---
1589 * The fluxes computed at surface, we could derive
1590 * inward and outward waves compositions via characteristic analysis,
1591 * and enforce inward outwaves to zero.
1592 * Ref. Bauer Paper.
1593 *
1594 * Ref. Bauer's paper:
1595 * "-Application of a Discontinuous Galerkin Method to Discretize
1596 * Acoustic Perturbation Equations, AIAA, May, 2011
1597 *
1598 * \author H.-J. Cheng <h-j.cheng@aia.rwth-aachen.de>
1599 * \date 2014-12-25
1600 *
1601 * labels:DG IMPORTANT: Please note that the mean flow matrix calculation is just a hack
1602 * right now and only works if this method is used for a boundary condition, not
1603 * for internal fluxes!
1604 */
1605template <>
1606inline void DgSysEqnAcousticPerturb<3>::calcRiemannRoe(const MFloat* nodeVarsL,
1607 const MFloat* nodeVarsR,
1608 const MFloat* stateL,
1609 const MFloat* stateR,
1610 const MInt noNodes1D,
1611 const MInt dirId,
1612 MFloat* flux) const {
1613 // TRACE();
1614
1615 const MInt nDim = 3;
1616 const MFloatTensor uL(const_cast<MFloat*>(stateL), noNodes1D, noNodes1D, noVars());
1617 const MFloatTensor uR(const_cast<MFloat*>(stateR), noNodes1D, noNodes1D, noVars());
1618 const MFloatTensor nL(const_cast<MFloat*>(nodeVarsL), noNodes1D, noNodes1D, noNodeVars());
1619 const MFloatTensor nR(const_cast<MFloat*>(nodeVarsR), noNodes1D, noNodes1D, noNodeVars());
1620 MFloatTensor f(flux, noNodes1D, noNodes1D, s_noVariables);
1621
1622 // TODO labels:DG use node variables instead
1623 // Define mean flow constants
1624 const MFloat c = 1.0;
1625 const MFloat rho = 1.0;
1626
1627 // Normal direction
1628 const std::array<MInt, nDim> n = {{dirId == 0 ? 1 : 0, dirId == 1 ? 1 : 0, dirId == 2 ? 1 : 0}};
1629
1630 // Characteristic boundary conditions
1631 // stateL and stateR multiplied by transformation matrix (see reference)
1632 for(MInt i = 0; i < noNodes1D; i++) {
1633 for(MInt j = 0; j < noNodes1D; j++) {
1634 // FIXME labels:DG,totest The following is just a hack and probably does not work when the
1635 // Roe Riemann solver is used for internal fluxes
1636 // Determine mean flow
1637 const std::array<MFloat, nDim> meanVelocity = {{0.5 * (nL(i, j, CV::U0) + nR(i, j, CV::U0)),
1638 0.5 * (nL(i, j, CV::V0) + nR(i, j, CV::V0)),
1639 0.5 * (nL(i, j, CV::W0) + nR(i, j, CV::W0))}};
1640
1641 // Mean flow matrix with transformation matrix (see ref.)
1642 const std::array<MFloat, nDim> um = {{meanVelocity[0] * n[0] + meanVelocity[1] * n[1] + meanVelocity[2] * n[2],
1643 meanVelocity[0] * n[1] + meanVelocity[1] * n[2] + meanVelocity[2] * n[0],
1644 meanVelocity[0] * n[2] + meanVelocity[1] * n[0] + meanVelocity[2] * n[1]}};
1645
1646 const std::array<MFloat, nDim> vecL = {
1647 {uL(i, j, CV::UU[0]) * n[0] + uL(i, j, CV::UU[1]) * n[1] + uL(i, j, CV::UU[2]) * n[2],
1648 uL(i, j, CV::UU[0]) * n[1] + uL(i, j, CV::UU[1]) * n[2] + uL(i, j, CV::UU[2]) * n[0],
1649 uL(i, j, CV::UU[0]) * n[2] + uL(i, j, CV::UU[1]) * n[0] + uL(i, j, CV::UU[2]) * n[1]}};
1650
1651 const std::array<MFloat, nDim> vecR = {
1652 {uR(i, j, CV::UU[0]) * n[0] + uR(i, j, CV::UU[1]) * n[1] + uR(i, j, CV::UU[2]) * n[2],
1653 uR(i, j, CV::UU[0]) * n[1] + uR(i, j, CV::UU[1]) * n[2] + uR(i, j, CV::UU[2]) * n[0],
1654 uR(i, j, CV::UU[0]) * n[2] + uR(i, j, CV::UU[1]) * n[0] + uR(i, j, CV::UU[2]) * n[1]}};
1655
1656 // Calculate characteristic waves
1657 // L1 and L2 composite from vecL and vecR
1658 const MFloat L1 =
1659 0.5
1660 * (uR(i, j, CV::P)
1661 + rho * c * (-vecL[0] + um[1] / (um[0] - c) * (-vecL[1]) + um[2] / (um[0] - c) * (-vecL[2])));
1662 const MFloat L2 =
1663 0.5
1664 * (uL(i, j, CV::P) + rho * c * (vecR[0] + um[1] / (um[0] + c) * vecR[1] + um[2] / (um[0] + c) * (vecR[2])));
1665 // Calculate Riemann flux
1666 f(i, j, CV::P) = (um[0] - c) * (L1) + (um[0] + c) * (L2);
1667 f(i, j, CV::UU[0]) = 1 / rho / c * (-1 * ((um[0] - c) * L1) + ((um[0] + c) * L2)) * n[0];
1668 f(i, j, CV::UU[1]) = 1 / rho / c * (-1 * ((um[0] - c) * L1) + ((um[0] + c) * L2)) * n[1];
1669 f(i, j, CV::UU[2]) = 1 / rho / c * (-1 * ((um[0] - c) * L1) + ((um[0] + c) * L2)) * n[2];
1670 }
1671 }
1672}
1673
1674
1685template <MInt nDim>
1686void DgSysEqnAcousticPerturb<nDim>::consToPrim(const MFloat* cons, MFloat* prim) const {
1687 // Very easy...
1688 std::copy(cons, cons + s_noVariables, prim);
1689}
1690
1691
1702template <MInt nDim>
1703void DgSysEqnAcousticPerturb<nDim>::primToCons(const MFloat* prim, MFloat* cons) const {
1704 // Very easy...
1705 std::copy(prim, prim + s_noVariables, cons);
1706}
1707
1708
1715template <MInt nDim>
1716void DgSysEqnAcousticPerturb<nDim>::getDefaultNodeVars(MFloat* const nodeVars) const {
1717 // Mean velocities
1718 std::copy_n(std::begin(m_meanVelocity), nDim, &nodeVars[CV::UU0[0]]);
1719 // Mean density
1720 nodeVars[CV::RHO0] = m_meanDensity;
1721 // Mean speed of sound
1722 nodeVars[CV::C0] = m_meanSpeedOfSound;
1723 // Derivatives of mean speed of sound, always zero
1724 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
1725}
1726
1727
1729template <MInt nDim>
1730void DgSysEqnAcousticPerturb<nDim>::getDefaultNodeVarsBody(MFloat* const nodeVars) const {
1731 // Mean velocities
1732 std::fill_n(&nodeVars[CV::UU0[0]], nDim, 0.0);
1733 // Mean density
1734 nodeVars[CV::RHO0] = m_meanDensity;
1735 // Mean speed of sound
1736 nodeVars[CV::C0] = m_meanSpeedOfSound;
1737 // Derivatives of mean speed of sound, always zero
1738 std::fill_n(&nodeVars[CV::DC0[0]], nDim, 0.0);
1739}
1740
1741
1743template <MInt nDim>
1744MBool DgSysEqnAcousticPerturb<nDim>::extendNodeVar(const MInt varId) const {
1745 // Do not extend...
1746 if(varId == CV::C0) return false; // mean speed of sound ...
1747 for(MInt i = 0; i < nDim; i++) {
1748 if(varId == CV::DC0[i]) { // and its derivatives
1749 return false;
1750 }
1751 }
1752 return true;
1753}
1754
1755
1756#endif // DGSYSEQNACOUSTICPERTURB_H_
DgSysEqnAcousticPerturb::calcRiemann
void calcRiemann(const MFloat *nodeVarsL, const MFloat *nodeVarsR, const MFloat *stateL, const MFloat *stateR, const MInt noNodes1D, const MInt dirId, MFloat *flux) const
Calculates the numerical flux at a surface given two states (left and right).
Definition: dgcartesiansyseqnacousticperturb.h:1423
DgSysEqnAcousticPerturb::getDefaultNodeVars
void getDefaultNodeVars(MFloat *const nodeVars) const
Return the default node variables.
Definition: dgcartesiansyseqnacousticperturb.h:1716
DgSysEqnAcousticPerturb::DgSysEqnAcousticPerturb
DgSysEqnAcousticPerturb(MInt solverId)
Constructor calls parent constructor & loads all necessary properties for this equation.
Definition: dgcartesiansyseqnacousticperturb.h:155
DgSysEqnAcousticPerturb::primToCons
void primToCons(const MFloat *prim, MFloat *cons) const
Calculates a set of conservative variables from a set of primitive variables.
Definition: dgcartesiansyseqnacousticperturb.h:1703
DgSysEqnAcousticPerturb::calcInitialCondition
void calcInitialCondition(const MFloat t, const MFloat *x, MFloat *const nodeVars, MFloat *u) const
Calculate the initial condition for a certain point in space.
Definition: dgcartesiansyseqnacousticperturb.h:393
DgSysEqnAcousticPerturb::calcFlux
void calcFlux(const MFloat *const nodeVars, const MFloat *const q, const MInt noNodes1D, MFloat *const flux) const
Calculates the physical fluxes in all dimensions for all integrations points within a cell.
Definition: dgcartesiansyseqnacousticperturb.h:827
DgSysEqnAcousticPerturb::s_primVarNames
static const MString s_primVarNames[s_noVariables]
Definition: dgcartesiansyseqnacousticperturb.h:68
DgSysEqnAcousticPerturb::cflFactor
MFloat cflFactor(const MInt polyDeg) const
Definition: dgcartesiansyseqnacousticperturb.h:942
DgSysEqnAcousticPerturb::calcRiemannLaxFriedich
void calcRiemannLaxFriedich(const MFloat *nodeVarsL, const MFloat *nodeVarsR, const MFloat *stateL, const MFloat *stateR, const MInt noNodes1D, const MInt dirId, MFloat *flux) const
Calculate Riemann flux using the local Lax-Friedrichs flux scheme.
Definition: dgcartesiansyseqnacousticperturb.h:1461
DgSysEqnAcousticPerturb::calcFlux1D
void calcFlux1D(const MFloat *const nodeVars, const MFloat *const q, const MInt noNodes1D, const MInt dirId, MFloat *const flux) const
Calculates the physical fluxes in dimension dirId for all integrations points on an element face.
Definition: dgcartesiansyseqnacousticperturb.h:896
DgSysEqnAcousticPerturb::m_cflFactor
std::array< std::array< MFloat, s_maxPolyDeg+1 >, 2 > m_cflFactor
Definition: dgcartesiansyseqnacousticperturb.h:85
DgSysEqnAcousticPerturb::getTimeStep
MFloat getTimeStep(const MFloat *nodeVars, const MFloat *const u, const MInt noNodes1D, const MFloat invJacobian, const MInt sbpMode) const
Calculate the time step for an explicit time stepping scheme for a given element.
Definition: dgcartesiansyseqnacousticperturb.h:1368
DgSysEqnAcousticPerturb::extendNodeVar
MBool extendNodeVar(const MInt varId) const
Return if the given node variable should be extended.
Definition: dgcartesiansyseqnacousticperturb.h:1744
DgSysEqnAcousticPerturb::getDefaultNodeVarsBody
void getDefaultNodeVarsBody(MFloat *const nodeVars) const
Return the default node variables inside a body.
Definition: dgcartesiansyseqnacousticperturb.h:1730
DgSysEqnAcousticPerturb::calcSource
void calcSource(const MFloat *const nodeVars, const MFloat *const u, const MInt noNodes1D, const MFloat t, const MFloat *const x, MFloat *const src) const
Calculates the source terms for all integration points within a cell.
Definition: dgcartesiansyseqnacousticperturb.h:969
DgSysEqnAcousticPerturb::calcRiemannRoe
void calcRiemannRoe(const MFloat *nodeVarsL, const MFloat *nodeVarsR, const MFloat *stateL, const MFloat *stateR, const MInt noNodes1D, const MInt dirId, MFloat *flux) const
DgSysEqnAcousticPerturb::consToPrim
void consToPrim(const MFloat *cons, MFloat *prim) const
Calculates a set of primitive variables from a set of conservative variables.
Definition: dgcartesiansyseqnacousticperturb.h:1686
DgSysEqnAcousticPerturb::calcSpongeSource
void calcSpongeSource(const MFloat *nodeVars, const MFloat *u, const MInt noNodes1D, const MFloat *eta, MFloat *src) const
Calculates the sponge source terms for all integration points within a cell.
Definition: dgcartesiansyseqnacousticperturb.h:1328
DgSysEqnAcousticPerturb::s_nodeVarNames
static const MString s_nodeVarNames[s_noNodeVars]
Definition: dgcartesiansyseqnacousticperturb.h:73
DgSysEqnAcousticPerturb::m_meanVelocity
std::array< MFloat, nDim > m_meanVelocity
Definition: dgcartesiansyseqnacousticperturb.h:75
DgSysEqnAcousticPerturb::s_consVarNames
static const MString s_consVarNames[s_noVariables]
Definition: dgcartesiansyseqnacousticperturb.h:67
DgSysEqn
Base class for concrete system-of-equations classes (CRTP).
Definition: dgcartesiansyseqn.h:28
ScratchSpace
This class is a ScratchSpace.
Definition: scratch.h:758
string2enum
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
Definition: enums.cpp:20
DG_TIMEINTEGRATION_CARPENTER_4_5
@ DG_TIMEINTEGRATION_CARPENTER_4_5
Definition: enums.h:321
DG_TIMEINTEGRATION_TOULORGEC_3_7
@ DG_TIMEINTEGRATION_TOULORGEC_3_7
Definition: enums.h:325
DG_TIMEINTEGRATION_TOULORGEC_4_8
@ DG_TIMEINTEGRATION_TOULORGEC_4_8
Definition: enums.h:322
DG_TIMEINTEGRATION_TOULORGEF_4_8
@ DG_TIMEINTEGRATION_TOULORGEF_4_8
Definition: enums.h:326
DG_TIMEINTEGRATION_NIEGEMANN_4_13
@ DG_TIMEINTEGRATION_NIEGEMANN_4_13
Definition: enums.h:324
DG_TIMEINTEGRATION_NIEGEMANN_4_14
@ DG_TIMEINTEGRATION_NIEGEMANN_4_14
Definition: enums.h:323
mTerm
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
POW2
constexpr Real POW2(const Real x)
Definition: functions.h:119
m_log
InfoOutFile m_log
Definition: globalvariables.cpp:57
MInt
int32_t MInt
Definition: maiatypes.h:62
MString
std::basic_string< char > MString
Definition: maiatypes.h:55
MFloat
double MFloat
Definition: maiatypes.h:52
MBool
bool MBool
Definition: maiatypes.h:58
maia::filter::slope::cos
T cos(const T a, const T b, const T x)
Cosine slope filter.
Definition: filter.h:125
DgSysEqnAcousticPerturb::CV::UU0
static constexpr const MInt UU0[nDim]
Definition: dgcartesiansyseqnacousticperturb.h:112
DgSysEqnAcousticPerturb::CV::DC0
static constexpr const MInt DC0[nDim]
Definition: dgcartesiansyseqnacousticperturb.h:117
DgSysEqnAcousticPerturb::CV::UU
static constexpr const MInt UU[nDim]
Definition: dgcartesiansyseqnacousticperturb.h:106
a
Definition: contexttypes.h:19