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