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 <math.h> 20 #include <ceed.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_dw[i+1] > prof_dw[i] and prof_dw[0] = 0 31 * If dw > prof_dw[-1], then the interpolation takes the values at prof_dw[-1] 32 * 33 * @param[in] dw Distance to the nearest wall 34 * @param[out] ubar Mean velocity at dw 35 * @param[out] cij Cholesky decomposition at dw 36 * @param[out] eps Turbulent dissipation at dw 37 * @param[out] lt Turbulent length scale at dw 38 * @param[in] stg_ctx STGShur14Context for the problem 39 */ 40 CEED_QFUNCTION_HELPER void InterpolateProfile(const CeedScalar dw, 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_dw = &stg_ctx->data[stg_ctx->offsets.prof_dw]; 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 (dw < prof_dw[i]) { 54 idx = i; 55 break; 56 } 57 } 58 59 if (idx > 0) { // y within the bounds of prof_dw 60 CeedScalar coeff = (dw - prof_dw[idx-1]) / (prof_dw[idx] - prof_dw[idx-1]); 61 62 //*INDENT-OFF* 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_dw 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 coefficients for STG 92 * 93 * Calculates q_n at a given distance to the wall 94 * 95 * @param[in] dw Distance to the nearest wall 96 * @param[in] eps Turbulent dissipation w/rt dw 97 * @param[in] lt Turbulent length scale w/rt dw 98 * @param[in] h Element lengths in coordinate directions 99 * @param[in] nu Dynamic Viscosity; 100 * @param[in] stg_ctx STGShur14Context for the problem 101 * @param[out] qn Spectrum coefficients, [nmodes] 102 */ 103 void CEED_QFUNCTION_HELPER(CalcSpectrum)(const CeedScalar dw, 104 const CeedScalar eps, const CeedScalar lt, const CeedScalar h[3], 105 const CeedScalar nu, CeedScalar qn[], const STGShur14Context stg_ctx) { 106 107 const CeedInt nmodes = stg_ctx->nmodes; 108 const CeedScalar *kappa = &stg_ctx->data[stg_ctx->offsets.kappa]; 109 110 const CeedScalar hmax = Max( Max(h[0], h[1]), h[2]); 111 const CeedScalar ke = dw==0 ? 1e16 : 2*M_PI/Min(2*dw, 3*lt); 112 const CeedScalar keta = 2*M_PI*pow(Cube(nu)/eps, -0.25); 113 const CeedScalar kcut = 114 M_PI/ Min( Max(Max(h[1], h[2]), 0.3*hmax) + 0.1*dw, hmax ); 115 CeedScalar fcut, feta, Ektot=0.0; 116 117 for(CeedInt n=0; n<nmodes; n++) { 118 feta = exp(-Square(12*kappa[n]/keta)); 119 fcut = exp( -pow(4*Max(kappa[n] - 0.9*kcut, 0)/kcut, 3.) ); 120 qn[n] = pow(kappa[n]/ke, 4.) 121 * pow(1 + 2.4*Square(kappa[n]/ke),-17./6)*feta*fcut; 122 qn[n] *= n==0 ? kappa[0] : kappa[n] - kappa[n-1]; 123 Ektot += qn[n]; 124 } 125 126 if (Ektot == 0) return; 127 for(CeedInt n=0; n<nmodes; n++) qn[n] /= Ektot; 128 } 129 130 /****************************************************** 131 * @brief Calculate u(x,t) for STG inflow condition 132 * 133 * @param[in] X Location to evaluate u(X,t) 134 * @param[in] t Time to evaluate u(X,t) 135 * @param[in] ubar Mean velocity at X 136 * @param[in] cij Cholesky decomposition at X 137 * @param[in] qn Wavemode amplitudes at X, [nmodes] 138 * @param[out] u Velocity at X and t 139 * @param[in] stg_ctx STGShur14Context for the problem 140 */ 141 void CEED_QFUNCTION_HELPER(STGShur14_Calc)(const CeedScalar X[3], 142 const CeedScalar t, const CeedScalar ubar[3], const CeedScalar cij[6], 143 const CeedScalar qn[], CeedScalar u[3], 144 const STGShur14Context stg_ctx) { 145 146 //*INDENT-OFF* 147 const CeedInt nmodes = stg_ctx->nmodes; 148 const CeedScalar *kappa = &stg_ctx->data[stg_ctx->offsets.kappa]; 149 const CeedScalar *phi = &stg_ctx->data[stg_ctx->offsets.phi]; 150 const CeedScalar *sigma = &stg_ctx->data[stg_ctx->offsets.sigma]; 151 const CeedScalar *d = &stg_ctx->data[stg_ctx->offsets.d]; 152 //*INDENT-ON* 153 CeedScalar xdotd, vp[3] = {0.}; 154 CeedScalar xhat[] = {0., X[1], X[2]}; 155 156 CeedPragmaSIMD 157 for(CeedInt n=0; n<nmodes; n++) { 158 xhat[0] = (X[0] - stg_ctx->u0*t)*Max(2*kappa[0]/kappa[n], 0.1); 159 xdotd = 0.; 160 for(CeedInt i=0; i<3; i++) xdotd += d[i*nmodes+n]*xhat[i]; 161 const CeedScalar cos_kxdp = cos(kappa[n]*xdotd + phi[n]); 162 vp[0] += sqrt(qn[n])*sigma[0*nmodes+n] * cos_kxdp; 163 vp[1] += sqrt(qn[n])*sigma[1*nmodes+n] * cos_kxdp; 164 vp[2] += sqrt(qn[n])*sigma[2*nmodes+n] * cos_kxdp; 165 } 166 for(CeedInt i=0; i<3; i++) vp[i] *= 2*sqrt(1.5); 167 168 u[0] = ubar[0] + cij[0]*vp[0]; 169 u[1] = ubar[1] + cij[3]*vp[0] + cij[1]*vp[1]; 170 u[2] = ubar[2] + cij[4]*vp[0] + cij[5]*vp[1] + cij[2]*vp[2]; 171 } 172 173 // Extrude the STGInflow profile through out the domain for an initial condition 174 CEED_QFUNCTION(ICsSTG)(void *ctx, CeedInt Q, 175 const CeedScalar *const *in, CeedScalar *const *out) { 176 // Inputs 177 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 178 179 // Outputs 180 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 181 182 const STGShur14Context stg_ctx = (STGShur14Context) ctx; 183 CeedScalar u[3], cij[6], eps, lt; 184 const CeedScalar theta0 = stg_ctx->theta0; 185 const CeedScalar P0 = stg_ctx->P0; 186 const CeedScalar cv = stg_ctx->newtonian_ctx.cv; 187 const CeedScalar cp = stg_ctx->newtonian_ctx.cp; 188 const CeedScalar Rd = cp - cv; 189 const CeedScalar rho = P0 / (Rd * theta0); 190 191 CeedPragmaSIMD 192 for(CeedInt i=0; i<Q; i++) { 193 InterpolateProfile(X[1][i], u, cij, &eps, <, stg_ctx); 194 195 q0[0][i] = rho; 196 q0[1][i] = u[0] * rho; 197 q0[2][i] = u[1] * rho; 198 q0[3][i] = u[2] * rho; 199 q0[4][i] = rho * (0.5 * Dot3(u, u) + cv * theta0); 200 } // End of Quadrature Point Loop 201 return 0; 202 } 203 204 /******************************************************************** 205 * @brief QFunction to calculate the inflow boundary condition 206 * 207 * This will loop through quadrature points, calculate the wavemode amplitudes 208 * at each location, then calculate the actual velocity. 209 */ 210 CEED_QFUNCTION(STGShur14_Inflow)(void *ctx, CeedInt Q, 211 const CeedScalar *const *in, 212 CeedScalar *const *out) { 213 214 //*INDENT-OFF* 215 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA]) in[0], 216 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA]) in[2], 217 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA]) in[3]; 218 219 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[0], 220 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[1]; 221 222 //*INDENT-ON* 223 224 const STGShur14Context stg_ctx = (STGShur14Context) ctx; 225 CeedScalar qn[STG_NMODES_MAX], u[3], ubar[3], cij[6], eps, lt; 226 const bool is_implicit = stg_ctx->is_implicit; 227 const bool mean_only = stg_ctx->mean_only; 228 const bool prescribe_T = stg_ctx->prescribe_T; 229 const CeedScalar dx = stg_ctx->dx; 230 const CeedScalar mu = stg_ctx->newtonian_ctx.mu; 231 const CeedScalar time = stg_ctx->time; 232 const CeedScalar theta0 = stg_ctx->theta0; 233 const CeedScalar P0 = stg_ctx->P0; 234 const CeedScalar cv = stg_ctx->newtonian_ctx.cv; 235 const CeedScalar cp = stg_ctx->newtonian_ctx.cp; 236 const CeedScalar Rd = cp - cv; 237 const CeedScalar gamma = cp/cv; 238 239 CeedPragmaSIMD 240 for(CeedInt i=0; i<Q; i++) { 241 const CeedScalar rho = prescribe_T ? q[0][i] : P0 / (Rd * theta0); 242 const CeedScalar x[] = { X[0][i], X[1][i], X[2][i] }; 243 const CeedScalar dXdx[2][3] = { 244 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 245 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 246 }; 247 248 CeedScalar h[3]; 249 h[0] = dx; 250 for (CeedInt j=1; j<3; j++) 251 h[j] = 2/sqrt(Square(dXdx[0][j]) + Square(dXdx[1][j])); 252 253 InterpolateProfile(X[1][i], ubar, cij, &eps, <, stg_ctx); 254 if (!mean_only) { 255 CalcSpectrum(X[1][i], eps, lt, h, mu/rho, qn, stg_ctx); 256 STGShur14_Calc(x, time, ubar, cij, qn, u, stg_ctx); 257 } else { 258 for (CeedInt j=0; j<3; j++) u[j] = ubar[j]; 259 } 260 261 const CeedScalar E_kinetic = .5 * rho * Dot3(u, u); 262 CeedScalar E_internal, P; 263 if (prescribe_T) { 264 // Temperature is being set weakly (theta0) and for constant cv this sets E_internal 265 E_internal = rho * cv * theta0; 266 // Find pressure using 267 P = rho * Rd * theta0; // interior rho with exterior T 268 } else { 269 E_internal = q[4][i] - E_kinetic; // uses prescribed rho and u, E from solution 270 P = E_internal * (gamma - 1.); 271 } 272 273 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 274 // ---- Normal vect 275 const CeedScalar norm[3] = {q_data_sur[1][i], 276 q_data_sur[2][i], 277 q_data_sur[3][i] 278 }; 279 280 const CeedScalar E = E_internal + E_kinetic; 281 282 // Velocity normal to the boundary 283 const CeedScalar u_normal = Dot3(norm, u); 284 285 // The Physics 286 // Zero v so all future terms can safely sum into it 287 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 288 289 // The Physics 290 // -- Density 291 v[0][i] -= wdetJb * rho * u_normal; 292 293 // -- Momentum 294 for (CeedInt j=0; j<3; j++) 295 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + 296 norm[j] * P); 297 298 // -- Total Energy Density 299 v[4][i] -= wdetJb * u_normal * (E + P); 300 301 jac_data_sur[0][i] = rho; 302 jac_data_sur[1][i] = u[0]; 303 jac_data_sur[2][i] = u[1]; 304 jac_data_sur[3][i] = u[2]; 305 jac_data_sur[4][i] = E; 306 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = 0.; 307 } 308 return 0; 309 } 310 311 CEED_QFUNCTION(STGShur14_Inflow_Jacobian)(void *ctx, CeedInt Q, 312 const CeedScalar *const *in, 313 CeedScalar *const *out) { 314 // *INDENT-OFF* 315 // Inputs 316 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 317 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 318 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 319 // Outputs 320 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 321 // *INDENT-ON* 322 const STGShur14Context stg_ctx = (STGShur14Context)ctx; 323 const bool implicit = stg_ctx->is_implicit; 324 const CeedScalar cv = stg_ctx->newtonian_ctx.cv; 325 const CeedScalar cp = stg_ctx->newtonian_ctx.cp; 326 const CeedScalar Rd = cp - cv; 327 const CeedScalar gamma = cp/cv; 328 329 const CeedScalar theta0 = stg_ctx->theta0; 330 const bool prescribe_T = stg_ctx->prescribe_T; 331 332 CeedPragmaSIMD 333 // Quadrature Point Loop 334 for (CeedInt i=0; i<Q; i++) { 335 // Setup 336 // Setup 337 // -- Interp-to-Interp q_data 338 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 339 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 340 // We can effect this by swapping the sign on this weight 341 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 342 343 // Calculate inflow values 344 CeedScalar velocity[3]; 345 for (CeedInt j=0; j<3; j++) velocity[j] = jac_data_sur[5+j][i]; 346 347 // enabling user to choose between weak T and weak rho inflow 348 CeedScalar drho, dE, dP; 349 if (prescribe_T) { 350 // rho should be from the current solution 351 drho = dq[0][i]; 352 CeedScalar dE_internal = drho * cv * theta0; 353 CeedScalar dE_kinetic = .5 * drho * Dot3(velocity, velocity); 354 dE = dE_internal + dE_kinetic; 355 dP = drho * Rd * theta0; // interior rho with exterior T 356 } else { // rho specified, E_internal from solution 357 drho = 0; 358 dE = dq[4][i]; 359 dP = dE * (gamma - 1.); 360 } 361 const CeedScalar norm[3] = {q_data_sur[1][i], 362 q_data_sur[2][i], 363 q_data_sur[3][i] 364 }; 365 366 const CeedScalar u_normal = Dot3(norm, velocity); 367 368 v[0][i] = - wdetJb * drho * u_normal; 369 for (int j=0; j<3; j++) 370 v[j+1][i] = -wdetJb * (drho * u_normal * velocity[j] + norm[j] * dP); 371 v[4][i] = - wdetJb * u_normal * (dE + dP); 372 } // End Quadrature Point Loop 373 return 0; 374 } 375 376 /******************************************************************** 377 * @brief QFunction to calculate the strongly enforce inflow BC 378 * 379 * This QF is for the strong application of STG via libCEED (rather than 380 * through the native PETSc `DMAddBoundary` -> `bcFunc` method. 381 */ 382 CEED_QFUNCTION(STGShur14_Inflow_StrongQF)(void *ctx, CeedInt Q, 383 const CeedScalar *const *in, CeedScalar *const *out) { 384 385 //*INDENT-OFF* 386 const CeedScalar (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA]) in[0], 387 (*coords)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA]) in[1], 388 (*scale) = (const CeedScalar(*)) in[2]; 389 390 CeedScalar(*bcval)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[0]; 391 //*INDENT-ON* 392 393 const STGShur14Context stg_ctx = (STGShur14Context) ctx; 394 CeedScalar qn[STG_NMODES_MAX], u[3], ubar[3], cij[6], eps, lt; 395 const bool mean_only = stg_ctx->mean_only; 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 rho = P0 / (Rd * theta0); 405 406 CeedPragmaSIMD 407 for(CeedInt i=0; i<Q; i++) { 408 const CeedScalar x[] = { coords[0][i], coords[1][i], coords[2][i] }; 409 const CeedScalar dXdx[2][3] = { 410 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 411 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 412 }; 413 414 CeedScalar h[3]; 415 h[0] = dx; 416 for (CeedInt j=1; j<3; j++) 417 h[j] = 2/sqrt(Square(dXdx[0][j]) + Square(dXdx[1][j])); 418 419 InterpolateProfile(coords[1][i], ubar, cij, &eps, <, stg_ctx); 420 if (!mean_only) { 421 CalcSpectrum(coords[1][i], eps, lt, h, mu/rho, qn, stg_ctx); 422 STGShur14_Calc(x, time, ubar, cij, qn, u, stg_ctx); 423 } else { 424 for (CeedInt j=0; j<3; j++) u[j] = ubar[j]; 425 } 426 427 bcval[0][i] = scale[i] * rho; 428 bcval[1][i] = scale[i] * rho * u[0]; 429 bcval[2][i] = scale[i] * rho * u[1]; 430 bcval[3][i] = scale[i] * rho * u[2]; 431 bcval[4][i] = 0.; 432 } 433 return 0; 434 } 435 436 #endif // stg_shur14_h 437