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 24 #ifndef M_PI 25 #define M_PI 3.14159265358979323846 26 #endif 27 28 #define STG_NMODES_MAX 1024 29 30 CEED_QFUNCTION_HELPER CeedScalar Max(CeedScalar a, CeedScalar b) { return a < b ? b : a; } 31 CEED_QFUNCTION_HELPER CeedScalar Min(CeedScalar a, CeedScalar b) { return a < b ? a : b; } 32 33 /* 34 * @brief Interpolate quantities from input profile to given location 35 * 36 * Assumed that prof_dw[i+1] > prof_dw[i] and prof_dw[0] = 0 37 * If dw > prof_dw[-1], then the interpolation takes the values at prof_dw[-1] 38 * 39 * @param[in] dw Distance to the nearest wall 40 * @param[out] ubar Mean velocity at dw 41 * @param[out] cij Cholesky decomposition at dw 42 * @param[out] eps Turbulent dissipation at dw 43 * @param[out] lt Turbulent length scale at dw 44 * @param[in] stg_ctx STGShur14Context for the problem 45 */ 46 CEED_QFUNCTION_HELPER void InterpolateProfile(const CeedScalar dw, 47 CeedScalar ubar[3], CeedScalar cij[6], CeedScalar *eps, CeedScalar *lt, 48 const STGShur14Context stg_ctx) { 49 50 const CeedInt nprofs = stg_ctx->nprofs; 51 const CeedScalar *prof_dw = &stg_ctx->data[stg_ctx->offsets.prof_dw]; 52 const CeedScalar *prof_eps = &stg_ctx->data[stg_ctx->offsets.eps]; 53 const CeedScalar *prof_lt = &stg_ctx->data[stg_ctx->offsets.lt]; 54 const CeedScalar *prof_ubar = &stg_ctx->data[stg_ctx->offsets.ubar]; 55 const CeedScalar *prof_cij = &stg_ctx->data[stg_ctx->offsets.cij]; 56 CeedInt idx=-1; 57 58 for(CeedInt i=0; i<nprofs; i++) { 59 if (dw < prof_dw[i]) { 60 idx = i; 61 break; 62 } 63 } 64 65 if (idx > 0) { // y within the bounds of prof_dw 66 CeedScalar coeff = (dw - prof_dw[idx-1]) / (prof_dw[idx] - prof_dw[idx-1]); 67 68 //*INDENT-OFF* 69 ubar[0] = prof_ubar[0*nprofs+idx-1] + coeff*( prof_ubar[0*nprofs+idx] - prof_ubar[0*nprofs+idx-1] ); 70 ubar[1] = prof_ubar[1*nprofs+idx-1] + coeff*( prof_ubar[1*nprofs+idx] - prof_ubar[1*nprofs+idx-1] ); 71 ubar[2] = prof_ubar[2*nprofs+idx-1] + coeff*( prof_ubar[2*nprofs+idx] - prof_ubar[2*nprofs+idx-1] ); 72 cij[0] = prof_cij[0*nprofs+idx-1] + coeff*( prof_cij[0*nprofs+idx] - prof_cij[0*nprofs+idx-1] ); 73 cij[1] = prof_cij[1*nprofs+idx-1] + coeff*( prof_cij[1*nprofs+idx] - prof_cij[1*nprofs+idx-1] ); 74 cij[2] = prof_cij[2*nprofs+idx-1] + coeff*( prof_cij[2*nprofs+idx] - prof_cij[2*nprofs+idx-1] ); 75 cij[3] = prof_cij[3*nprofs+idx-1] + coeff*( prof_cij[3*nprofs+idx] - prof_cij[3*nprofs+idx-1] ); 76 cij[4] = prof_cij[4*nprofs+idx-1] + coeff*( prof_cij[4*nprofs+idx] - prof_cij[4*nprofs+idx-1] ); 77 cij[5] = prof_cij[5*nprofs+idx-1] + coeff*( prof_cij[5*nprofs+idx] - prof_cij[5*nprofs+idx-1] ); 78 *eps = prof_eps[idx-1] + coeff*( prof_eps[idx] - prof_eps[idx-1] ); 79 *lt = prof_lt[idx-1] + coeff*( prof_lt[idx] - prof_lt[idx-1] ); 80 //*INDENT-ON* 81 } else { // y outside bounds of prof_dw 82 ubar[0] = prof_ubar[1*nprofs-1]; 83 ubar[1] = prof_ubar[2*nprofs-1]; 84 ubar[2] = prof_ubar[3*nprofs-1]; 85 cij[0] = prof_cij[1*nprofs-1]; 86 cij[1] = prof_cij[2*nprofs-1]; 87 cij[2] = prof_cij[3*nprofs-1]; 88 cij[3] = prof_cij[4*nprofs-1]; 89 cij[4] = prof_cij[5*nprofs-1]; 90 cij[5] = prof_cij[6*nprofs-1]; 91 *eps = prof_eps[nprofs-1]; 92 *lt = prof_lt[nprofs-1]; 93 } 94 } 95 96 /* 97 * @brief Calculate spectrum coefficients for STG 98 * 99 * Calculates q_n at a given distance to the wall 100 * 101 * @param[in] dw Distance to the nearest wall 102 * @param[in] eps Turbulent dissipation w/rt dw 103 * @param[in] lt Turbulent length scale w/rt dw 104 * @param[in] h Element lengths in coordinate directions 105 * @param[in] nu Dynamic Viscosity; 106 * @param[in] stg_ctx STGShur14Context for the problem 107 * @param[out] qn Spectrum coefficients, [nmodes] 108 */ 109 void CEED_QFUNCTION_HELPER(CalcSpectrum)(const CeedScalar dw, 110 const CeedScalar eps, const CeedScalar lt, const CeedScalar h[3], 111 const CeedScalar nu, CeedScalar qn[], const STGShur14Context stg_ctx) { 112 113 const CeedInt nmodes = stg_ctx->nmodes; 114 const CeedScalar *kappa = &stg_ctx->data[stg_ctx->offsets.kappa]; 115 116 const CeedScalar hmax = Max( Max(h[0], h[1]), h[2]); 117 const CeedScalar ke = 2*M_PI/Min(2*dw, 3*lt); 118 const CeedScalar keta = 2*M_PI*pow(pow(nu,3.0)/eps, -0.25); 119 const CeedScalar kcut = 120 M_PI/ Min( Max(Max(h[1], h[2]), 0.3*hmax) + 0.1*dw, hmax ); 121 CeedScalar fcut, feta, Ektot=0.0; 122 123 for(CeedInt n=0; n<nmodes; n++) { 124 feta = exp(-Square(12*kappa[n]/keta)); 125 fcut = exp( -pow(4*Max(kappa[n] - 0.9*kcut, 0)/kcut, 3.) ); 126 qn[n] = pow(kappa[n]/ke, 4.) 127 * pow(1 + 2.4*Square(kappa[n]/ke),-17./6)*feta*fcut; 128 qn[n] *= n==0 ? kappa[0] : kappa[n] - kappa[n-1]; 129 Ektot += qn[n]; 130 } 131 132 if (Ektot == 0) return; 133 for(CeedInt n=0; n<nmodes; n++) qn[n] /= Ektot; 134 } 135 136 /****************************************************** 137 * @brief Calculate u(x,t) for STG inflow condition 138 * 139 * @param[in] X Location to evaluate u(X,t) 140 * @param[in] t Time to evaluate u(X,t) 141 * @param[in] ubar Mean velocity at X 142 * @param[in] cij Cholesky decomposition at X 143 * @param[in] qn Wavemode amplitudes at X, [nmodes] 144 * @param[out] u Velocity at X and t 145 * @param[in] stg_ctx STGShur14Context for the problem 146 */ 147 void CEED_QFUNCTION_HELPER(STGShur14_Calc)(const CeedScalar X[3], 148 const CeedScalar t, const CeedScalar ubar[3], const CeedScalar cij[6], 149 const CeedScalar qn[], CeedScalar u[3], 150 const STGShur14Context stg_ctx) { 151 152 //*INDENT-OFF* 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 //*INDENT-ON* 159 CeedScalar xdotd, vp[3] = {0.}; 160 CeedScalar xhat[] = {0., X[1], X[2]}; 161 162 CeedPragmaSIMD 163 for(CeedInt n=0; n<nmodes; n++) { 164 xhat[0] = (X[0] - stg_ctx->u0*t)*Max(2*kappa[0]/kappa[n], 0.1); 165 xdotd = 0.; 166 for(CeedInt i=0; i<3; i++) xdotd += d[i*nmodes+n]*xhat[i]; 167 const CeedScalar cos_kxdp = cos(kappa[n]*xdotd + phi[n]); 168 vp[0] += sqrt(qn[n])*sigma[0*nmodes+n] * cos_kxdp; 169 vp[1] += sqrt(qn[n])*sigma[1*nmodes+n] * cos_kxdp; 170 vp[2] += sqrt(qn[n])*sigma[2*nmodes+n] * cos_kxdp; 171 } 172 for(CeedInt i=0; i<3; i++) vp[i] *= 2*sqrt(1.5); 173 174 u[0] = ubar[0] + cij[0]*vp[0]; 175 u[1] = ubar[1] + cij[3]*vp[0] + cij[1]*vp[1]; 176 u[2] = ubar[2] + cij[4]*vp[0] + cij[5]*vp[1] + cij[2]*vp[2]; 177 } 178 179 /******************************************************************** 180 * @brief QFunction to calculate the inflow boundary condition 181 * 182 * This will loop through quadrature points, calculate the wavemode amplitudes 183 * at each location, then calculate the actual velocity. 184 */ 185 CEED_QFUNCTION(STGShur14_Inflow)(void *ctx, CeedInt Q, 186 const CeedScalar *const *in, 187 CeedScalar *const *out) { 188 189 //*INDENT-OFF* 190 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA]) in[0], 191 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA]) in[1], 192 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA]) in[2]; 193 194 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[0]; 195 196 //*INDENT-ON* 197 198 const STGShur14Context stg_ctx = (STGShur14Context) ctx; 199 CeedScalar qn[STG_NMODES_MAX], u[3], ubar[3], cij[6], eps, lt; 200 const bool is_implicit = stg_ctx->is_implicit; 201 const bool mean_only = stg_ctx->mean_only; 202 const bool prescribe_T = stg_ctx->prescribe_T; 203 const CeedScalar dx = stg_ctx->dx; 204 const CeedScalar mu = stg_ctx->newtonian_ctx.mu; 205 const CeedScalar time = stg_ctx->time; 206 const CeedScalar theta0 = stg_ctx->theta0; 207 const CeedScalar P0 = stg_ctx->P0; 208 const CeedScalar cv = stg_ctx->newtonian_ctx.cv; 209 const CeedScalar cp = stg_ctx->newtonian_ctx.cp; 210 const CeedScalar Rd = cp - cv; 211 const CeedScalar gamma = cp/cv; 212 213 CeedPragmaSIMD 214 for(CeedInt i=0; i<Q; i++) { 215 const CeedScalar rho = prescribe_T ? q[0][i] : P0 / (Rd * theta0); 216 const CeedScalar x[] = { X[0][i], X[1][i], X[2][i] }; 217 const CeedScalar dXdx[2][3] = { 218 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 219 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 220 }; 221 222 CeedScalar h[3]; 223 for (CeedInt j=0; j<3; j++) 224 h[j] = 2/sqrt(dXdx[0][j]*dXdx[0][j] + dXdx[1][j]*dXdx[1][j]); 225 h[0] = dx; 226 227 InterpolateProfile(X[1][i], ubar, cij, &eps, <, stg_ctx); 228 if (!mean_only) { 229 CalcSpectrum(X[1][i], eps, lt, h, mu/rho, qn, stg_ctx); 230 STGShur14_Calc(x, time, ubar, cij, qn, u, stg_ctx); 231 } else { 232 for (CeedInt j=0; j<3; j++) u[j] = ubar[j]; 233 } 234 235 const CeedScalar E_kinetic = .5 * rho * (u[0]*u[0] + 236 u[1]*u[1] + 237 u[2]*u[2]); 238 CeedScalar E_internal, P; 239 if (prescribe_T) { 240 // Temperature is being set weakly (theta0) and for constant cv this sets E_internal 241 E_internal = rho * cv * theta0; 242 // Find pressure using 243 P = rho * Rd * theta0; // interior rho with exterior T 244 } else { 245 E_internal = q[4][i] - E_kinetic; // uses prescribed rho and u, E from solution 246 P = E_internal * (gamma - 1.); 247 } 248 249 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 250 // ---- Normal vect 251 const CeedScalar norm[3] = {q_data_sur[1][i], 252 q_data_sur[2][i], 253 q_data_sur[3][i] 254 }; 255 256 const CeedScalar E = E_internal + E_kinetic; 257 258 // Velocity normal to the boundary 259 const CeedScalar u_normal = norm[0]*u[0] + 260 norm[1]*u[1] + 261 norm[2]*u[2]; 262 // The Physics 263 // Zero v so all future terms can safely sum into it 264 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 265 266 // The Physics 267 // -- Density 268 v[0][i] -= wdetJb * rho * u_normal; 269 270 // -- Momentum 271 for (CeedInt j=0; j<3; j++) 272 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + 273 norm[j] * P); 274 275 // -- Total Energy Density 276 v[4][i] -= wdetJb * u_normal * (E + P); 277 } 278 return 0; 279 } 280 281 /* Compute boundary integral for strong STG enforcement 282 * 283 * This assumes that density is set strongly and temperature is allowed to 284 * float 285 */ 286 CEED_QFUNCTION(STGShur14_Inflow_Strong)(void *ctx, CeedInt Q, 287 const CeedScalar *const *in, 288 CeedScalar *const *out) { 289 290 //*INDENT-OFF* 291 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA]) in[0], 292 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA]) in[1]; 293 294 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[0]; 295 296 //*INDENT-ON* 297 298 const STGShur14Context stg_ctx = (STGShur14Context) ctx; 299 const bool is_implicit = stg_ctx->is_implicit; 300 const CeedScalar cv = stg_ctx->newtonian_ctx.cv; 301 const CeedScalar cp = stg_ctx->newtonian_ctx.cp; 302 const CeedScalar gamma = cp/cv; 303 304 CeedPragmaSIMD 305 for(CeedInt i=0; i<Q; i++) { 306 const CeedScalar rho = q[0][i]; 307 const CeedScalar u[] = {q[1][i]/rho, q[2][i]/rho, q[3][i]/rho}; 308 const CeedScalar E_kinetic = .5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]); 309 const CeedScalar E_internal = q[4][i] - E_kinetic; 310 const CeedScalar P = E_internal * (gamma - 1.); 311 312 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 313 // ---- Normal vect 314 const CeedScalar norm[3] = {q_data_sur[1][i], 315 q_data_sur[2][i], 316 q_data_sur[3][i] 317 }; 318 319 const CeedScalar E = E_internal + E_kinetic; 320 321 // Velocity normal to the boundary 322 const CeedScalar u_normal = norm[0]*u[0] + 323 norm[1]*u[1] + 324 norm[2]*u[2]; 325 // The Physics 326 // Zero v so all future terms can safely sum into it 327 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 328 329 // The Physics 330 // -- Density 331 v[0][i] -= wdetJb * rho * u_normal; 332 333 // -- Momentum 334 for (CeedInt j=0; j<3; j++) 335 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + 336 norm[j] * P); 337 338 // -- Total Energy Density 339 v[4][i] -= wdetJb * u_normal * (E + P); 340 } 341 return 0; 342 } 343 344 #endif // stg_shur14_h 345