1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Implementation of the Synthetic Turbulence Generation (STG) algorithm
6 /// presented in Shur et al. 2014
7 //
8 /// SetupSTG_Rand reads in the input files and fills in STGShur14Context.
9 /// Then STGShur14_CalcQF is run over quadrature points.
10 /// Before the program exits, TearDownSTG is run to free the memory of the allocated arrays.
11 #include <ceed/types.h>
12 #ifndef CEED_RUNNING_JIT_PASS
13 #include <math.h>
14 #endif
15
16 #include "newtonian_state.h"
17 #include "setupgeo_helpers.h"
18 #include "stg_shur14_type.h"
19 #include "utils.h"
20
21 #define STG_NMODES_MAX 1024
22
23 /*
24 * @brief Interpolate quantities from input profile to given location
25 *
26 * Assumed that prof_wd[i+1] > prof_wd[i] and prof_wd[0] = 0
27 * If wall_dist > prof_wd[-1], then the interpolation takes the values at prof_wd[-1]
28 *
29 * @param[in] wall_dist Distance to the nearest wall
30 * @param[out] ubar Mean velocity at wall_dist
31 * @param[out] cij Cholesky decomposition at wall_dist
32 * @param[out] eps Turbulent dissipation at wall_dist
33 * @param[out] lt Turbulent length scale at wall_dist
34 * @param[in] stg_ctx STGShur14Context for the problem
35 */
InterpolateProfile(const CeedScalar wall_dist,CeedScalar ubar[3],CeedScalar cij[6],CeedScalar * eps,CeedScalar * lt,const StgShur14Context stg_ctx)36 CEED_QFUNCTION_HELPER void InterpolateProfile(const CeedScalar wall_dist, CeedScalar ubar[3], CeedScalar cij[6], CeedScalar *eps, CeedScalar *lt,
37 const StgShur14Context stg_ctx) {
38 const CeedInt nprofs = stg_ctx->nprofs;
39 const CeedScalar *prof_wd = &stg_ctx->data[stg_ctx->offsets.wall_dist];
40 const CeedScalar *prof_eps = &stg_ctx->data[stg_ctx->offsets.eps];
41 const CeedScalar *prof_lt = &stg_ctx->data[stg_ctx->offsets.lt];
42 const CeedScalar *prof_ubar = &stg_ctx->data[stg_ctx->offsets.ubar];
43 const CeedScalar *prof_cij = &stg_ctx->data[stg_ctx->offsets.cij];
44 CeedInt idx = -1;
45
46 for (CeedInt i = 0; i < nprofs; i++) {
47 if (wall_dist < prof_wd[i]) {
48 idx = i;
49 break;
50 }
51 }
52
53 if (idx > 0) { // y within the bounds of prof_wd
54 CeedScalar coeff = (wall_dist - prof_wd[idx - 1]) / (prof_wd[idx] - prof_wd[idx - 1]);
55
56 ubar[0] = prof_ubar[0 * nprofs + idx - 1] + coeff * (prof_ubar[0 * nprofs + idx] - prof_ubar[0 * nprofs + idx - 1]);
57 ubar[1] = prof_ubar[1 * nprofs + idx - 1] + coeff * (prof_ubar[1 * nprofs + idx] - prof_ubar[1 * nprofs + idx - 1]);
58 ubar[2] = prof_ubar[2 * nprofs + idx - 1] + coeff * (prof_ubar[2 * nprofs + idx] - prof_ubar[2 * nprofs + idx - 1]);
59 cij[0] = prof_cij[0 * nprofs + idx - 1] + coeff * (prof_cij[0 * nprofs + idx] - prof_cij[0 * nprofs + idx - 1]);
60 cij[1] = prof_cij[1 * nprofs + idx - 1] + coeff * (prof_cij[1 * nprofs + idx] - prof_cij[1 * nprofs + idx - 1]);
61 cij[2] = prof_cij[2 * nprofs + idx - 1] + coeff * (prof_cij[2 * nprofs + idx] - prof_cij[2 * nprofs + idx - 1]);
62 cij[3] = prof_cij[3 * nprofs + idx - 1] + coeff * (prof_cij[3 * nprofs + idx] - prof_cij[3 * nprofs + idx - 1]);
63 cij[4] = prof_cij[4 * nprofs + idx - 1] + coeff * (prof_cij[4 * nprofs + idx] - prof_cij[4 * nprofs + idx - 1]);
64 cij[5] = prof_cij[5 * nprofs + idx - 1] + coeff * (prof_cij[5 * nprofs + idx] - prof_cij[5 * nprofs + idx - 1]);
65 *eps = prof_eps[idx - 1] + coeff * (prof_eps[idx] - prof_eps[idx - 1]);
66 *lt = prof_lt[idx - 1] + coeff * (prof_lt[idx] - prof_lt[idx - 1]);
67 } else { // y outside bounds of prof_wd
68 ubar[0] = prof_ubar[1 * nprofs - 1];
69 ubar[1] = prof_ubar[2 * nprofs - 1];
70 ubar[2] = prof_ubar[3 * nprofs - 1];
71 cij[0] = prof_cij[1 * nprofs - 1];
72 cij[1] = prof_cij[2 * nprofs - 1];
73 cij[2] = prof_cij[3 * nprofs - 1];
74 cij[3] = prof_cij[4 * nprofs - 1];
75 cij[4] = prof_cij[5 * nprofs - 1];
76 cij[5] = prof_cij[6 * nprofs - 1];
77 *eps = prof_eps[nprofs - 1];
78 *lt = prof_lt[nprofs - 1];
79 }
80 }
81
82 /*
83 * @brief Calculate spectrum coefficient, qn
84 *
85 * Calculates q_n at a given distance to the wall
86 *
87 * @param[in] kappa nth wavenumber
88 * @param[in] dkappa Difference between wavenumbers
89 * @param[in] keta Dissipation wavenumber
90 * @param[in] kcut Mesh-induced cutoff wavenumber
91 * @param[in] ke Energy-containing wavenumber
92 * @param[in] Ektot_inv Inverse of total turbulent kinetic energy of spectrum
93 * @returns qn Spectrum coefficient
94 */
Calc_qn(const CeedScalar kappa,const CeedScalar dkappa,const CeedScalar keta,const CeedScalar kcut,const CeedScalar ke,const CeedScalar Ektot_inv)95 CEED_QFUNCTION_HELPER CeedScalar Calc_qn(const CeedScalar kappa, const CeedScalar dkappa, const CeedScalar keta, const CeedScalar kcut,
96 const CeedScalar ke, const CeedScalar Ektot_inv) {
97 const CeedScalar feta_x_fcut = exp(-Square(12 * kappa / keta) - Cube(4 * Max(kappa - 0.9 * kcut, 0) / kcut));
98 return pow(kappa / ke, 4.) * pow(1 + 2.4 * Square(kappa / ke), -17. / 6) * feta_x_fcut * dkappa * Ektot_inv;
99 }
100
101 // Calculate hmax, ke, keta, and kcut
SpectrumConstants(const CeedScalar wall_dist,const CeedScalar eps,const CeedScalar lt,const CeedScalar hNodSep[3],const CeedScalar nu,CeedScalar * hmax,CeedScalar * ke,CeedScalar * keta,CeedScalar * kcut)102 CEED_QFUNCTION_HELPER void SpectrumConstants(const CeedScalar wall_dist, const CeedScalar eps, const CeedScalar lt, const CeedScalar hNodSep[3],
103 const CeedScalar nu, CeedScalar *hmax, CeedScalar *ke, CeedScalar *keta, CeedScalar *kcut) {
104 *hmax = Max(Max(hNodSep[0], hNodSep[1]), hNodSep[2]);
105 *ke = wall_dist == 0 ? 1e16 : 2 * M_PI / Min(2 * wall_dist, 3 * lt);
106 *keta = 2 * M_PI * pow(Cube(nu) / eps, -0.25);
107 *kcut = M_PI / Min(Max(Max(hNodSep[1], hNodSep[2]), 0.3 * (*hmax)) + 0.1 * wall_dist, *hmax);
108 }
109
110 /*
111 * @brief Calculate spectrum coefficients for STG
112 *
113 * Calculates q_n at a given distance to the wall
114 *
115 * @param[in] wall_dist Distance to the nearest wall
116 * @param[in] eps Turbulent dissipation w/rt wall_dist
117 * @param[in] lt Turbulent length scale w/rt wall_dist
118 * @param[in] h_node_sep Element lengths in coordinate directions
119 * @param[in] nu Dynamic Viscosity;
120 * @param[in] stg_ctx STGShur14Context for the problem
121 * @param[out] qn Spectrum coefficients, [nmodes]
122 */
CalcSpectrum(const CeedScalar wall_dist,const CeedScalar eps,const CeedScalar lt,const CeedScalar h_node_sep[3],const CeedScalar nu,CeedScalar qn[],const StgShur14Context stg_ctx)123 CEED_QFUNCTION_HELPER void CalcSpectrum(const CeedScalar wall_dist, const CeedScalar eps, const CeedScalar lt, const CeedScalar h_node_sep[3],
124 const CeedScalar nu, CeedScalar qn[], const StgShur14Context stg_ctx) {
125 const CeedInt nmodes = stg_ctx->nmodes;
126 const CeedScalar *kappa = &stg_ctx->data[stg_ctx->offsets.kappa];
127 CeedScalar hmax, ke, keta, kcut, Ektot = 0.0;
128
129 SpectrumConstants(wall_dist, eps, lt, h_node_sep, nu, &hmax, &ke, &keta, &kcut);
130
131 for (CeedInt n = 0; n < nmodes; n++) {
132 const CeedScalar dkappa = n == 0 ? kappa[0] : kappa[n] - kappa[n - 1];
133 qn[n] = Calc_qn(kappa[n], dkappa, keta, kcut, ke, 1.0);
134 Ektot += qn[n];
135 }
136
137 if (Ektot == 0) return;
138 for (CeedInt n = 0; n < nmodes; n++) qn[n] /= Ektot;
139 }
140
141 /******************************************************
142 * @brief Calculate u(x,t) for STG inflow condition
143 *
144 * @param[in] X Location to evaluate u(X,t)
145 * @param[in] t Time to evaluate u(X,t)
146 * @param[in] ubar Mean velocity at X
147 * @param[in] cij Cholesky decomposition at X
148 * @param[in] qn Wavemode amplitudes at X, [nmodes]
149 * @param[out] u Velocity at X and t
150 * @param[in] stg_ctx STGShur14Context for the problem
151 */
StgShur14Calc(const CeedScalar X[3],const CeedScalar t,const CeedScalar ubar[3],const CeedScalar cij[6],const CeedScalar qn[],CeedScalar u[3],const StgShur14Context stg_ctx)152 CEED_QFUNCTION_HELPER void StgShur14Calc(const CeedScalar X[3], const CeedScalar t, const CeedScalar ubar[3], const CeedScalar cij[6],
153 const CeedScalar qn[], CeedScalar u[3], const StgShur14Context stg_ctx) {
154 const CeedInt nmodes = stg_ctx->nmodes;
155 const CeedScalar *kappa = &stg_ctx->data[stg_ctx->offsets.kappa];
156 const CeedScalar *phi = &stg_ctx->data[stg_ctx->offsets.phi];
157 const CeedScalar *sigma = &stg_ctx->data[stg_ctx->offsets.sigma];
158 const CeedScalar *d = &stg_ctx->data[stg_ctx->offsets.d];
159 CeedScalar xdotd, vp[3] = {0.};
160 CeedScalar xhat[] = {0., X[1], X[2]};
161
162 CeedPragmaSIMD for (CeedInt n = 0; n < nmodes; n++) {
163 xhat[0] = (X[0] - stg_ctx->u0 * t) * Max(2 * kappa[0] / kappa[n], 0.1);
164 xdotd = 0.;
165 for (CeedInt i = 0; i < 3; i++) xdotd += d[i * nmodes + n] * xhat[i];
166 const CeedScalar cos_kxdp = cos(kappa[n] * xdotd + phi[n]);
167 vp[0] += sqrt(qn[n]) * sigma[0 * nmodes + n] * cos_kxdp;
168 vp[1] += sqrt(qn[n]) * sigma[1 * nmodes + n] * cos_kxdp;
169 vp[2] += sqrt(qn[n]) * sigma[2 * nmodes + n] * cos_kxdp;
170 }
171 for (CeedInt i = 0; i < 3; i++) vp[i] *= 2 * sqrt(1.5);
172
173 u[0] = ubar[0] + cij[0] * vp[0];
174 u[1] = ubar[1] + cij[3] * vp[0] + cij[1] * vp[1];
175 u[2] = ubar[2] + cij[4] * vp[0] + cij[5] * vp[1] + cij[2] * vp[2];
176 }
177
178 /******************************************************
179 * @brief Calculate u(x,t) for STG inflow condition
180 *
181 * @param[in] X Location to evaluate u(X,t)
182 * @param[in] t Time to evaluate u(X,t)
183 * @param[in] ubar Mean velocity at X
184 * @param[in] cij Cholesky decomposition at X
185 * @param[in] Ektot Total spectrum energy at this location
186 * @param[in] h_node_sep Element size in 3 directions
187 * @param[in] wall_dist Distance to closest wall
188 * @param[in] eps Turbulent dissipation
189 * @param[in] lt Turbulent length scale
190 * @param[out] u Velocity at X and t
191 * @param[in] stg_ctx STGShur14Context for the problem
192 */
StgShur14Calc_PrecompEktot(const CeedScalar X[3],const CeedScalar t,const CeedScalar ubar[3],const CeedScalar cij[6],const CeedScalar Ektot,const CeedScalar h_node_sep[3],const CeedScalar wall_dist,const CeedScalar eps,const CeedScalar lt,const CeedScalar nu,CeedScalar u[3],const StgShur14Context stg_ctx)193 CEED_QFUNCTION_HELPER void StgShur14Calc_PrecompEktot(const CeedScalar X[3], const CeedScalar t, const CeedScalar ubar[3], const CeedScalar cij[6],
194 const CeedScalar Ektot, const CeedScalar h_node_sep[3], const CeedScalar wall_dist,
195 const CeedScalar eps, const CeedScalar lt, const CeedScalar nu, CeedScalar u[3],
196 const StgShur14Context stg_ctx) {
197 const CeedInt nmodes = stg_ctx->nmodes;
198 const CeedScalar *kappa = &stg_ctx->data[stg_ctx->offsets.kappa];
199 const CeedScalar *phi = &stg_ctx->data[stg_ctx->offsets.phi];
200 const CeedScalar *sigma = &stg_ctx->data[stg_ctx->offsets.sigma];
201 const CeedScalar *d = &stg_ctx->data[stg_ctx->offsets.d];
202 CeedScalar hmax, ke, keta, kcut;
203 SpectrumConstants(wall_dist, eps, lt, h_node_sep, nu, &hmax, &ke, &keta, &kcut);
204 CeedScalar xdotd, vp[3] = {0.};
205 CeedScalar xhat[] = {0., X[1], X[2]};
206
207 CeedPragmaSIMD for (CeedInt n = 0; n < nmodes; n++) {
208 xhat[0] = (X[0] - stg_ctx->u0 * t) * Max(2 * kappa[0] / kappa[n], 0.1);
209 xdotd = 0.;
210 for (CeedInt i = 0; i < 3; i++) xdotd += d[i * nmodes + n] * xhat[i];
211 const CeedScalar cos_kxdp = cos(kappa[n] * xdotd + phi[n]);
212 const CeedScalar dkappa = n == 0 ? kappa[0] : kappa[n] - kappa[n - 1];
213 const CeedScalar qn = Calc_qn(kappa[n], dkappa, keta, kcut, ke, Ektot);
214 vp[0] += sqrt(qn) * sigma[0 * nmodes + n] * cos_kxdp;
215 vp[1] += sqrt(qn) * sigma[1 * nmodes + n] * cos_kxdp;
216 vp[2] += sqrt(qn) * sigma[2 * nmodes + n] * cos_kxdp;
217 }
218 for (CeedInt i = 0; i < 3; i++) vp[i] *= 2 * sqrt(1.5);
219
220 u[0] = ubar[0] + cij[0] * vp[0];
221 u[1] = ubar[1] + cij[3] * vp[0] + cij[1] * vp[1];
222 u[2] = ubar[2] + cij[4] * vp[0] + cij[5] * vp[1] + cij[2] * vp[2];
223 }
224
225 /**
226 @brief Calculate the element length scales based on dXdx
227
228 WARNING: This assumes the reference domain is [-1,1], which is not true for tetrahedral elements
229
230 @param[in] dXdx Inverse mapping Jacobian, d\xi/dx
231 @param[in] scale Scale factor for the element lengths
232 @param[out] lengths The element lengths in each cartesian direction
233 **/
CalculateElementLengths(CeedScalar dXdx[3][3],CeedScalar scale,CeedScalar lengths[3])234 CEED_QFUNCTION_HELPER void CalculateElementLengths(CeedScalar dXdx[3][3], CeedScalar scale, CeedScalar lengths[3]) {
235 for (CeedInt j = 0; j < 3; j++) lengths[j] = 2 / sqrt(Square(dXdx[0][j]) + Square(dXdx[1][j]) + Square(dXdx[2][j]));
236 ScaleN(lengths, scale, 3);
237 }
238
239 // Create preprocessed input for the stg calculation
240 //
241 // stg_data[0] = 1 / Ektot (inverse of total spectrum energy)
StgShur14Preprocess(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)242 CEED_QFUNCTION(StgShur14Preprocess)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
243 const CeedScalar *dXdx_q = in[0];
244 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
245
246 CeedScalar(*stg_data) = (CeedScalar(*))out[0];
247
248 CeedScalar ubar[3], cij[6], eps, lt;
249 const StgShur14Context stg_ctx = (StgShur14Context)ctx;
250 const CeedScalar mu = stg_ctx->newt_ctx.gas.mu;
251 const CeedScalar theta0 = stg_ctx->theta0;
252 const CeedScalar P0 = stg_ctx->P0;
253 const CeedScalar Rd = GasConstant(stg_ctx->newt_ctx.gas);
254 const CeedScalar rho = P0 / (Rd * theta0);
255 const CeedScalar nu = mu / rho;
256
257 const CeedInt nmodes = stg_ctx->nmodes;
258 const CeedScalar *kappa = &stg_ctx->data[stg_ctx->offsets.kappa];
259 CeedScalar hmax, ke, keta, kcut;
260
261 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
262 const CeedScalar wall_dist = x[1][i];
263 CeedScalar dXdx[3][3], h_node_sep[3];
264 StoredValuesUnpack(Q, i, 0, 9, dXdx_q, (CeedScalar *)dXdx);
265
266 CalculateElementLengths(dXdx, stg_ctx->h_scale_factor, h_node_sep);
267 InterpolateProfile(wall_dist, ubar, cij, &eps, <, stg_ctx);
268 SpectrumConstants(wall_dist, eps, lt, h_node_sep, nu, &hmax, &ke, &keta, &kcut);
269
270 // Calculate total TKE per spectrum
271 CeedScalar Ek_tot = 0;
272 CeedPragmaSIMD for (CeedInt n = 0; n < nmodes; n++) {
273 const CeedScalar dkappa = n == 0 ? kappa[0] : kappa[n] - kappa[n - 1];
274 Ek_tot += Calc_qn(kappa[n], dkappa, keta, kcut, ke, 1.0);
275 }
276 // avoid underflowed and poorly defined spectrum coefficients
277 stg_data[i] = Ek_tot != 0 ? 1 / Ek_tot : 0;
278 }
279 return 0;
280 }
281
282 // Extrude the STGInflow profile through out the domain for an initial condition
ICsStg(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)283 CEED_QFUNCTION(ICsStg)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
284 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
285 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1];
286 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
287
288 const StgShur14Context stg_ctx = (StgShur14Context)ctx;
289 const NewtonianIdealGasContext context = &stg_ctx->newt_ctx;
290 const NewtonianIGProperties gas = context->gas;
291 CeedScalar qn[STG_NMODES_MAX], u[3], ubar[3], cij[6], eps, lt;
292 const CeedScalar dx = stg_ctx->dx;
293 const CeedScalar time = stg_ctx->time;
294 const CeedScalar theta0 = stg_ctx->theta0;
295 const CeedScalar P0 = stg_ctx->P0;
296 const CeedScalar rho = P0 / (GasConstant(gas) * theta0);
297 const CeedScalar nu = gas.mu / rho;
298
299 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
300 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]};
301 CeedScalar dXdx[3][3];
302 CeedScalar h_node_sep[3];
303
304 InvertMappingJacobian_3D(Q, i, J, dXdx, NULL);
305 CalculateElementLengths(dXdx, stg_ctx->h_scale_factor, h_node_sep);
306 h_node_sep[0] = dx * stg_ctx->h_scale_factor;
307
308 InterpolateProfile(x_i[1], ubar, cij, &eps, <, stg_ctx);
309 if (stg_ctx->use_fluctuating_IC) {
310 CalcSpectrum(x_i[1], eps, lt, h_node_sep, nu, qn, stg_ctx);
311 StgShur14Calc(x_i, time, ubar, cij, qn, u, stg_ctx);
312 } else {
313 for (CeedInt j = 0; j < 3; j++) u[j] = ubar[j];
314 }
315
316 CeedScalar Y[5] = {P0, u[0], u[1], u[2], theta0}, q[5];
317 State s = StateFromY(gas, Y);
318 StateToQ(gas, s, q, context->state_var);
319 for (CeedInt j = 0; j < 5; j++) {
320 q0[j][i] = q[j];
321 }
322 }
323 return 0;
324 }
325
326 /********************************************************************
327 * @brief QFunction to calculate the inflow boundary condition
328 *
329 * This will loop through quadrature points, calculate the wavemode amplitudes
330 * at each location, then calculate the actual velocity.
331 */
StgShur14Inflow(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)332 CEED_QFUNCTION(StgShur14Inflow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
333 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
334 const CeedScalar(*q_data_sur) = in[2];
335 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3];
336
337 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
338 CeedScalar(*jac_data_sur) = out[1];
339
340 const StgShur14Context stg_ctx = (StgShur14Context)ctx;
341 const NewtonianIGProperties gas = stg_ctx->newt_ctx.gas;
342 CeedScalar qn[STG_NMODES_MAX], u[3], ubar[3], cij[6], eps, lt;
343 const bool is_implicit = stg_ctx->is_implicit;
344 const bool mean_only = stg_ctx->mean_only;
345 const bool prescribe_T = stg_ctx->prescribe_T;
346 const CeedScalar dx = stg_ctx->dx;
347 const CeedScalar mu = gas.mu;
348 const CeedScalar time = stg_ctx->time;
349 const CeedScalar theta0 = stg_ctx->theta0;
350 const CeedScalar P0 = stg_ctx->P0;
351 const CeedScalar cv = gas.cv;
352 const CeedScalar Rd = GasConstant(gas);
353 const CeedScalar gamma = HeatCapacityRatio(gas);
354
355 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
356 const CeedScalar rho = prescribe_T ? q[0][i] : P0 / (Rd * theta0);
357 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]};
358 CeedScalar wdetJb, dXdx[2][3], normal[3];
359 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal);
360 wdetJb *= is_implicit ? -1. : 1.;
361
362 CeedScalar h_node_sep[3];
363 h_node_sep[0] = dx;
364 for (CeedInt j = 1; j < 3; j++) h_node_sep[j] = 2 / sqrt(Square(dXdx[0][j]) + Square(dXdx[1][j]));
365 ScaleN(h_node_sep, stg_ctx->h_scale_factor, 3);
366
367 InterpolateProfile(X[1][i], ubar, cij, &eps, <, stg_ctx);
368 if (!mean_only) {
369 CalcSpectrum(X[1][i], eps, lt, h_node_sep, mu / rho, qn, stg_ctx);
370 StgShur14Calc(x, time, ubar, cij, qn, u, stg_ctx);
371 } else {
372 for (CeedInt j = 0; j < 3; j++) u[j] = ubar[j];
373 }
374
375 const CeedScalar E_kinetic = .5 * rho * Dot3(u, u);
376 CeedScalar E_internal, P;
377 if (prescribe_T) {
378 // Temperature is being set weakly (theta0) and for constant cv this sets E_internal
379 E_internal = rho * cv * theta0;
380 // Find pressure using
381 P = rho * Rd * theta0; // interior rho with exterior T
382 } else {
383 E_internal = q[4][i] - E_kinetic; // uses prescribed rho and u, E from solution
384 P = E_internal * (gamma - 1.);
385 }
386
387 const CeedScalar E = E_internal + E_kinetic;
388
389 // Velocity normal to the boundary
390 const CeedScalar u_normal = Dot3(normal, u);
391
392 // The Physics
393 // Zero v so all future terms can safely sum into it
394 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0.;
395
396 // The Physics
397 // -- Density
398 v[0][i] -= wdetJb * rho * u_normal;
399
400 // -- Momentum
401 for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho * u_normal * u[j] + normal[j] * P);
402
403 // -- Total Energy Density
404 v[4][i] -= wdetJb * u_normal * (E + P);
405
406 const CeedScalar U[] = {rho, u[0], u[1], u[2], E}, kmstress[6] = {0.};
407 StoredValuesPack(Q, i, 0, 5, U, jac_data_sur);
408 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur);
409 }
410 return 0;
411 }
412
StgShur14Inflow_Jacobian(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)413 CEED_QFUNCTION(StgShur14Inflow_Jacobian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
414 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
415 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
416 const CeedScalar(*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4];
417 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
418
419 const StgShur14Context stg_ctx = (StgShur14Context)ctx;
420 const NewtonianIGProperties gas = stg_ctx->newt_ctx.gas;
421 const bool implicit = stg_ctx->is_implicit;
422 const CeedScalar cv = gas.cv;
423 const CeedScalar Rd = GasConstant(gas);
424 const CeedScalar gamma = HeatCapacityRatio(gas);
425
426 const CeedScalar theta0 = stg_ctx->theta0;
427 const bool prescribe_T = stg_ctx->prescribe_T;
428
429 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
430 // Setup
431 // -- Interp-to-Interp q_data
432 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q).
433 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q).
434 // We can effect this by swapping the sign on this weight
435 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i];
436
437 // Calculate inflow values
438 CeedScalar velocity[3];
439 for (CeedInt j = 0; j < 3; j++) velocity[j] = jac_data_sur[5 + j][i];
440 // TODO This is almost certainly a bug. Velocity isn't stored here, only 0s.
441
442 // enabling user to choose between weak T and weak rho inflow
443 CeedScalar drho, dE, dP;
444 if (prescribe_T) {
445 // rho should be from the current solution
446 drho = dq[0][i];
447 CeedScalar dE_internal = drho * cv * theta0;
448 CeedScalar dE_kinetic = .5 * drho * Dot3(velocity, velocity);
449 dE = dE_internal + dE_kinetic;
450 dP = drho * Rd * theta0; // interior rho with exterior T
451 } else { // rho specified, E_internal from solution
452 drho = 0;
453 dE = dq[4][i];
454 dP = dE * (gamma - 1.);
455 }
456 const CeedScalar normal[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]};
457
458 const CeedScalar u_normal = Dot3(normal, velocity);
459
460 v[0][i] = -wdetJb * drho * u_normal;
461 for (int j = 0; j < 3; j++) v[j + 1][i] = -wdetJb * (drho * u_normal * velocity[j] + normal[j] * dP);
462 v[4][i] = -wdetJb * u_normal * (dE + dP);
463 }
464 return 0;
465 }
466
467 /********************************************************************
468 * @brief QFunction to calculate the strongly enforce inflow BC
469 *
470 * This QF is for the strong application of STG via libCEED (rather than
471 * through the native PETSc `DMAddBoundary` -> `bcFunc` method.
472 */
StgShur14InflowStrongQF(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)473 CEED_QFUNCTION(StgShur14InflowStrongQF)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
474 const CeedScalar *dXdx_q = in[0];
475 const CeedScalar(*coords)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
476 const CeedScalar(*scale) = (const CeedScalar(*))in[2];
477 const CeedScalar(*inv_Ektotal) = (const CeedScalar(*))in[3];
478 CeedScalar(*bcval)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
479
480 const StgShur14Context stg_ctx = (StgShur14Context)ctx;
481 const NewtonianIGProperties gas = stg_ctx->newt_ctx.gas;
482 CeedScalar u[3], ubar[3], cij[6], eps, lt;
483 const bool mean_only = stg_ctx->mean_only;
484 const CeedScalar time = stg_ctx->time;
485 const CeedScalar theta0 = stg_ctx->theta0;
486 const CeedScalar P0 = stg_ctx->P0;
487 const CeedScalar rho = P0 / (GasConstant(gas) * theta0);
488 const CeedScalar nu = gas.mu / rho;
489
490 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
491 const CeedScalar x[] = {coords[0][i], coords[1][i], coords[2][i]};
492 CeedScalar dXdx[3][3], h_node_sep[3];
493 StoredValuesUnpack(Q, i, 0, 9, dXdx_q, (CeedScalar *)dXdx);
494
495 CalculateElementLengths(dXdx, stg_ctx->h_scale_factor, h_node_sep);
496 InterpolateProfile(coords[1][i], ubar, cij, &eps, <, stg_ctx);
497 if (!mean_only) {
498 if (1) {
499 StgShur14Calc_PrecompEktot(x, time, ubar, cij, inv_Ektotal[i], h_node_sep, x[1], eps, lt, nu, u, stg_ctx);
500 } else { // Original way
501 CeedScalar qn[STG_NMODES_MAX];
502 CalcSpectrum(coords[1][i], eps, lt, h_node_sep, nu, qn, stg_ctx);
503 StgShur14Calc(x, time, ubar, cij, qn, u, stg_ctx);
504 }
505 } else {
506 for (CeedInt j = 0; j < 3; j++) u[j] = ubar[j];
507 }
508
509 CeedScalar Y[5] = {P0, u[0], u[1], u[2], theta0}, q[5];
510 State s = StateFromY(gas, Y);
511 StateToQ(gas, s, q, stg_ctx->newt_ctx.state_var);
512 switch (stg_ctx->newt_ctx.state_var) {
513 case STATEVAR_CONSERVATIVE:
514 q[4] = 0.; // Don't set energy
515 break;
516 case STATEVAR_PRIMITIVE:
517 q[0] = 0; // Don't set pressure
518 break;
519 case STATEVAR_ENTROPY:
520 q[0] = 0; // Don't set V_density
521 break;
522 }
523 for (CeedInt j = 0; j < 5; j++) {
524 bcval[j][i] = scale[i] * q[j];
525 }
526 }
527 return 0;
528 }
529