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