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