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