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