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