1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Shock tube initial condition and Euler equation operator for Navier-Stokes 19 /// example using PETSc - modified from eulervortex.h 20 21 // Model from: 22 // On the Order of Accuracy and Numerical Performance of Two Classes of 23 // Finite Volume WENO Schemes, Zhang, Zhang, and Shu (2011). 24 25 #ifndef shocktube_h 26 #define shocktube_h 27 28 #include <math.h> 29 30 #ifndef M_PI 31 #define M_PI 3.14159265358979323846 32 #endif 33 34 #ifndef setup_context_struct 35 #define setup_context_struct 36 typedef struct SetupContext_ *SetupContext; 37 struct SetupContext_ { 38 CeedScalar theta0; 39 CeedScalar thetaC; 40 CeedScalar P0; 41 CeedScalar N; 42 CeedScalar cv; 43 CeedScalar cp; 44 CeedScalar g[3]; 45 CeedScalar rc; 46 CeedScalar lx; 47 CeedScalar ly; 48 CeedScalar lz; 49 CeedScalar center[3]; 50 CeedScalar dc_axis[3]; 51 CeedScalar wind[3]; 52 CeedScalar time; 53 CeedScalar mid_point; 54 CeedScalar P_high; 55 CeedScalar rho_high; 56 CeedScalar P_low; 57 CeedScalar rho_low; 58 int wind_type; // See WindType: 0=ROTATION, 1=TRANSLATION 59 int bubble_type; // See BubbleType: 0=SPHERE, 1=CYLINDER 60 int bubble_continuity_type; // See BubbleContinuityType: 0=SMOOTH, 1=BACK_SHARP 2=THICK 61 }; 62 #endif 63 64 typedef struct ShockTubeContext_ *ShockTubeContext; 65 struct ShockTubeContext_ { 66 CeedScalar Cyzb; 67 CeedScalar Byzb; 68 CeedScalar c_tau; 69 bool implicit; 70 bool yzb; 71 int stabilization; 72 }; 73 74 // ***************************************************************************** 75 // This function sets the initial conditions 76 // 77 // Temperature: 78 // T = P / (rho * R) 79 // Density: 80 // rho = 1.0 if x <= mid_point 81 // = 0.125 if x > mid_point 82 // Pressure: 83 // P = 1.0 if x <= mid_point 84 // = 0.1 if x > mid_point 85 // Velocity: 86 // u = 0 87 // Velocity/Momentum Density: 88 // Ui = rho ui 89 // Total Energy: 90 // E = P / (gamma - 1) + rho (u u)/2 91 // 92 // Constants: 93 // cv , Specific heat, constant volume 94 // cp , Specific heat, constant pressure 95 // mid_point , Location of initial domain mid_point 96 // gamma = cp / cv, Specific heat ratio 97 // 98 // ***************************************************************************** 99 100 // ***************************************************************************** 101 // This helper function provides support for the exact, time-dependent solution 102 // (currently not implemented) and IC formulation for Euler traveling vortex 103 // ***************************************************************************** 104 CEED_QFUNCTION_HELPER int Exact_ShockTube(CeedInt dim, CeedScalar time, 105 const CeedScalar X[], CeedInt Nf, CeedScalar q[], 106 void *ctx) { 107 108 // Context 109 const SetupContext context = (SetupContext)ctx; 110 const CeedScalar mid_point = context->mid_point; // Midpoint of the domain 111 const CeedScalar P_high = context->P_high; // Driver section pressure 112 const CeedScalar rho_high = context->rho_high; // Driver section density 113 const CeedScalar P_low = context->P_low; // Driven section pressure 114 const CeedScalar rho_low = context->rho_low; // Driven section density 115 116 // Setup 117 const CeedScalar gamma = 1.4; // ratio of specific heats 118 const CeedScalar x = X[0]; // Coordinates 119 120 CeedScalar rho, P, u[3] = {0.}; 121 122 // Initial Conditions 123 if (x <= mid_point) { 124 rho = rho_high; 125 P = P_high; 126 } else { 127 rho = rho_low; 128 P = P_low; 129 } 130 131 // Assign exact solution 132 q[0] = rho; 133 q[1] = rho * u[0]; 134 q[2] = rho * u[1]; 135 q[3] = rho * u[2]; 136 q[4] = P / (gamma-1.0) + rho * (u[0]*u[0]) / 2.; 137 138 // Return 139 return 0; 140 } 141 142 // ***************************************************************************** 143 // Helper function for computing flux Jacobian 144 // ***************************************************************************** 145 CEED_QFUNCTION_HELPER void ConvectiveFluxJacobian_Euler(CeedScalar dF[3][5][5], 146 const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 147 const CeedScalar gamma) { 148 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; // Velocity square 149 for (CeedInt i=0; i<3; i++) { // Jacobian matrices for 3 directions 150 for (CeedInt j=0; j<3; j++) { // Rows of each Jacobian matrix 151 dF[i][j+1][0] = ((i==j) ? ((gamma-1.)*(u_sq/2.)) : 0.) - u[i]*u[j]; 152 for (CeedInt k=0; k<3; k++) { // Columns of each Jacobian matrix 153 dF[i][0][k+1] = ((i==k) ? 1. : 0.); 154 dF[i][j+1][k+1] = ((j==k) ? u[i] : 0.) + 155 ((i==k) ? u[j] : 0.) - 156 ((i==j) ? u[k] : 0.) * (gamma-1.); 157 dF[i][4][k+1] = ((i==k) ? (E*gamma/rho - (gamma-1.)*u_sq/2.) : 0.) - 158 (gamma-1.)*u[i]*u[k]; 159 } 160 dF[i][j+1][4] = ((i==j) ? (gamma-1.) : 0.); 161 } 162 dF[i][4][0] = u[i] * ((gamma-1.)*u_sq - E*gamma/rho); 163 dF[i][4][4] = u[i] * gamma; 164 } 165 } 166 167 // ***************************************************************************** 168 // Helper function for calculating the covariant length scale in the direction 169 // of some 3 element input vector 170 // 171 // Where 172 // vec = vector that length is measured in the direction of 173 // h = covariant element length along vec 174 // ***************************************************************************** 175 CEED_QFUNCTION_HELPER CeedScalar Covariant_length_along_vector( 176 CeedScalar vec[3], const CeedScalar dXdx[3][3]) { 177 178 CeedScalar vec_norm = sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]); 179 CeedScalar vec_dot_jacobian[3] = {0.0}; 180 for (CeedInt i=0; i<3; i++) { 181 for (CeedInt j=0; j<3; j++) { 182 vec_dot_jacobian[i] += dXdx[j][i]*vec[i]; 183 } 184 } 185 CeedScalar norm_vec_dot_jacobian = sqrt(vec_dot_jacobian[0]*vec_dot_jacobian[0]+ 186 vec_dot_jacobian[1]*vec_dot_jacobian[1]+ 187 vec_dot_jacobian[2]*vec_dot_jacobian[2]); 188 CeedScalar h = 2.0 * vec_norm / norm_vec_dot_jacobian; 189 return h; 190 } 191 192 193 // ***************************************************************************** 194 // Helper function for computing Tau elements (stabilization constant) 195 // Model from: 196 // Stabilized Methods for Compressible Flows, Hughes et al 2010 197 // 198 // Spatial criterion #2 - Tau is a 3x3 diagonal matrix 199 // Tau[i] = c_tau h[i] Xi(Pe) / rho(A[i]) (no sum) 200 // 201 // Where 202 // c_tau = stabilization constant (0.5 is reported as "optimal") 203 // h[i] = 2 length(dxdX[i]) 204 // Pe = Peclet number ( Pe = sqrt(u u) / dot(dXdx,u) diffusivity ) 205 // Xi(Pe) = coth Pe - 1. / Pe (1. at large local Peclet number ) 206 // rho(A[i]) = spectral radius of the convective flux Jacobian i, 207 // wave speed in direction i 208 // ***************************************************************************** 209 CEED_QFUNCTION_HELPER void Tau_spatial(CeedScalar Tau_x[3], 210 const CeedScalar dXdx[3][3], const CeedScalar u[3], 211 const CeedScalar sound_speed, const CeedScalar c_tau) { 212 for (int i=0; i<3; i++) { 213 // length of element in direction i 214 CeedScalar h = 2 / sqrt(dXdx[0][i]*dXdx[0][i] + dXdx[1][i]*dXdx[1][i] + 215 dXdx[2][i]*dXdx[2][i]); 216 // fastest wave in direction i 217 CeedScalar fastest_wave = fabs(u[i]) + sound_speed; 218 Tau_x[i] = c_tau * h / fastest_wave; 219 } 220 } 221 222 // ***************************************************************************** 223 // This QFunction sets the initial conditions for shock tube 224 // ***************************************************************************** 225 CEED_QFUNCTION(ICsShockTube)(void *ctx, CeedInt Q, 226 const CeedScalar *const *in, CeedScalar *const *out) { 227 // Inputs 228 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 229 230 // Outputs 231 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 232 233 CeedPragmaSIMD 234 // Quadrature Point Loop 235 for (CeedInt i=0; i<Q; i++) { 236 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 237 CeedScalar q[5]; 238 239 Exact_ShockTube(3, 0., x, 5, q, ctx); 240 241 for (CeedInt j=0; j<5; j++) 242 q0[j][i] = q[j]; 243 } // End of Quadrature Point Loop 244 245 // Return 246 return 0; 247 } 248 249 // ***************************************************************************** 250 // This QFunction implements the following formulation of Euler equations 251 // with explicit time stepping method 252 // 253 // This is 3D Euler for compressible gas dynamics in conservation 254 // form with state variables of density, momentum density, and total 255 // energy density. 256 // 257 // State Variables: q = ( rho, U1, U2, U3, E ) 258 // rho - Mass Density 259 // Ui - Momentum Density, Ui = rho ui 260 // E - Total Energy Density, E = P / (gamma - 1) + rho (u u)/2 261 // 262 // Euler Equations: 263 // drho/dt + div( U ) = 0 264 // dU/dt + div( rho (u x u) + P I3 ) = 0 265 // dE/dt + div( (E + P) u ) = 0 266 // 267 // Equation of State: 268 // P = (gamma - 1) (E - rho (u u) / 2) 269 // 270 // Constants: 271 // cv , Specific heat, constant volume 272 // cp , Specific heat, constant pressure 273 // g , Gravity 274 // gamma = cp / cv, Specific heat ratio 275 // ***************************************************************************** 276 CEED_QFUNCTION(EulerShockTube)(void *ctx, CeedInt Q, 277 const CeedScalar *const *in, CeedScalar *const *out) { 278 // *INDENT-OFF* 279 // Inputs 280 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 281 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 282 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 283 // Outputs 284 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 285 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 286 287 const CeedScalar gamma = 1.4; 288 289 ShockTubeContext context = (ShockTubeContext)ctx; 290 const CeedScalar Cyzb = context->Cyzb; 291 const CeedScalar Byzb = context->Byzb; 292 const CeedScalar c_tau = context->c_tau; 293 294 CeedPragmaSIMD 295 // Quadrature Point Loop 296 for (CeedInt i=0; i<Q; i++) { 297 // *INDENT-OFF* 298 // Setup 299 // -- Interp in 300 const CeedScalar rho = q[0][i]; 301 const CeedScalar u[3] = {q[1][i] / rho, 302 q[2][i] / rho, 303 q[3][i] / rho 304 }; 305 const CeedScalar E = q[4][i]; 306 const CeedScalar drho[3] = {dq[0][0][i], 307 dq[1][0][i], 308 dq[2][0][i] 309 }; 310 const CeedScalar dU[3][3] = {{dq[0][1][i], 311 dq[1][1][i], 312 dq[2][1][i]}, 313 {dq[0][2][i], 314 dq[1][2][i], 315 dq[2][2][i]}, 316 {dq[0][3][i], 317 dq[1][3][i], 318 dq[2][3][i]} 319 }; 320 const CeedScalar dE[3] = {dq[0][4][i], 321 dq[1][4][i], 322 dq[2][4][i] 323 }; 324 // -- Interp-to-Interp q_data 325 const CeedScalar wdetJ = q_data[0][i]; 326 // -- Interp-to-Grad q_data 327 // ---- Inverse of change of coordinate matrix: X_i,j 328 // *INDENT-OFF* 329 const CeedScalar dXdx[3][3] = {{q_data[1][i], 330 q_data[2][i], 331 q_data[3][i]}, 332 {q_data[4][i], 333 q_data[5][i], 334 q_data[6][i]}, 335 {q_data[7][i], 336 q_data[8][i], 337 q_data[9][i]} 338 }; 339 // dU/dx 340 CeedScalar du[3][3] = {{0}}; 341 CeedScalar drhodx[3] = {0}; 342 CeedScalar dEdx[3] = {0}; 343 CeedScalar dUdx[3][3] = {{0}}; 344 CeedScalar dXdxdXdxT[3][3] = {{0}}; 345 for (int j=0; j<3; j++) { 346 for (int k=0; k<3; k++) { 347 du[j][k] = (dU[j][k] - drho[k]*u[j]) / rho; 348 drhodx[j] += drho[k] * dXdx[k][j]; 349 dEdx[j] += dE[k] * dXdx[k][j]; 350 for (int l=0; l<3; l++) { 351 dUdx[j][k] += dU[j][l] * dXdx[l][k]; 352 dXdxdXdxT[j][k] += dXdx[j][l]*dXdx[k][l]; //dXdx_j,k * dXdx_k,j 353 } 354 } 355 } 356 357 // *INDENT-ON* 358 const CeedScalar 359 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 360 E_internal = E - E_kinetic, 361 P = E_internal * (gamma - 1); // P = pressure 362 363 // The Physics 364 // Zero v and dv so all future terms can safely sum into it 365 for (int j=0; j<5; j++) { 366 v[j][i] = 0; 367 for (int k=0; k<3; k++) 368 dv[k][j][i] = 0; 369 } 370 371 // -- Density 372 // ---- u rho 373 for (int j=0; j<3; j++) 374 dv[j][0][i] += wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 375 rho*u[2]*dXdx[j][2]); 376 // -- Momentum 377 // ---- rho (u x u) + P I3 378 for (int j=0; j<3; j++) 379 for (int k=0; k<3; k++) 380 dv[k][j+1][i] += wdetJ*((rho*u[j]*u[0] + (j==0?P:0))*dXdx[k][0] + 381 (rho*u[j]*u[1] + (j==1?P:0))*dXdx[k][1] + 382 (rho*u[j]*u[2] + (j==2?P:0))*dXdx[k][2]); 383 // -- Total Energy Density 384 // ---- (E + P) u 385 for (int j=0; j<3; j++) 386 dv[j][4][i] += wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 387 u[2]*dXdx[j][2]); 388 389 // -- YZB stabilization 390 if (context->yzb) { 391 CeedScalar drho_norm = 0.0; // magnitude of the density gradient 392 CeedScalar j_vec[3] = {0.0}; // unit vector aligned with the density gradient 393 CeedScalar h_shock = 0.0; // element lengthscale 394 CeedScalar acoustic_vel = 0.0; // characteristic velocity, acoustic speed 395 CeedScalar tau_shock = 0.0; // timescale 396 CeedScalar nu_shock = 0.0; // artificial diffusion 397 398 // Unit vector aligned with the density gradient 399 drho_norm = sqrt(drhodx[0]*drhodx[0] + drhodx[1]*drhodx[1] + 400 drhodx[2]*drhodx[2]); 401 for (int j=0; j<3; j++) 402 j_vec[j] = drhodx[j] / (drho_norm + 1e-20); 403 404 if (drho_norm == 0.0) { 405 nu_shock = 0.0; 406 } else { 407 h_shock = Covariant_length_along_vector(j_vec, dXdx); 408 h_shock /= Cyzb; 409 acoustic_vel = sqrt(gamma*P/rho); 410 tau_shock = h_shock / (2*acoustic_vel) * pow(drho_norm * h_shock / rho, Byzb); 411 nu_shock = fabs(tau_shock * acoustic_vel * acoustic_vel); 412 } 413 414 for (int j=0; j<3; j++) 415 dv[j][0][i] -= wdetJ * nu_shock * drhodx[j]; 416 417 for (int k=0; k<3; k++) 418 for (int j=0; j<3; j++) 419 dv[j][k][i] -= wdetJ * nu_shock * du[k][j]; 420 421 for (int j=0; j<3; j++) 422 dv[j][4][i] -= wdetJ * nu_shock * dEdx[j]; 423 } 424 425 // Stabilization 426 // Need the Jacobian for the advective fluxes for stabilization 427 // indexed as: jacob_F_conv[direction][flux component][solution component] 428 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}}; 429 ConvectiveFluxJacobian_Euler(jacob_F_conv, rho, u, E, gamma); 430 431 432 // dqdx collects drhodx, dUdx and dEdx in one vector 433 CeedScalar dqdx[5][3]; 434 for (int j=0; j<3; j++) { 435 dqdx[0][j] = drhodx[j]; 436 dqdx[4][j] = dEdx[j]; 437 for (int k=0; k<3; k++) 438 dqdx[k+1][j] = dUdx[k][j]; 439 } 440 441 // strong_conv = dF/dq * dq/dx (Strong convection) 442 CeedScalar strong_conv[5] = {0}; 443 for (int j=0; j<3; j++) 444 for (int k=0; k<5; k++) 445 for (int l=0; l<5; l++) 446 strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; 447 448 // Stabilization 449 // -- Tau elements 450 const CeedScalar sound_speed = sqrt(gamma * P / rho); 451 CeedScalar Tau_x[3] = {0.}; 452 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau); 453 454 CeedScalar stab[5][3] = {0}; 455 switch (context->stabilization) { 456 case 0: // Galerkin 457 break; 458 case 1: // SU 459 for (int j=0; j<3; j++) 460 for (int k=0; k<5; k++) 461 for (int l=0; l<5; l++) { 462 stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l]; 463 } 464 for (int j=0; j<5; j++) 465 for (int k=0; k<3; k++) 466 dv[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 467 stab[j][1] * dXdx[k][1] + 468 stab[j][2] * dXdx[k][2]); 469 break; 470 } 471 472 } // End Quadrature Point Loop 473 474 // Return 475 return 0; 476 } 477 478 #endif // shocktube_h 479