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 #ifndef shocktube_context_struct 65 #define shocktube_context_struct 66 typedef struct ShockTubeContext_ *ShockTubeContext; 67 struct ShockTubeContext_ { 68 CeedScalar Cyzb; 69 CeedScalar Byzb; 70 CeedScalar c_tau; 71 bool implicit; 72 bool yzb; 73 int stabilization; 74 }; 75 #endif 76 77 // ***************************************************************************** 78 // This function sets the initial conditions 79 // 80 // Temperature: 81 // T = P / (rho * R) 82 // Density: 83 // rho = 1.0 if x <= mid_point 84 // = 0.125 if x > mid_point 85 // Pressure: 86 // P = 1.0 if x <= mid_point 87 // = 0.1 if x > mid_point 88 // Velocity: 89 // u = 0 90 // Velocity/Momentum Density: 91 // Ui = rho ui 92 // Total Energy: 93 // E = P / (gamma - 1) + rho (u u)/2 94 // 95 // Constants: 96 // cv , Specific heat, constant volume 97 // cp , Specific heat, constant pressure 98 // mid_point , Location of initial domain mid_point 99 // gamma = cp / cv, Specific heat ratio 100 // 101 // ***************************************************************************** 102 103 // ***************************************************************************** 104 // This helper function provides support for the exact, time-dependent solution 105 // (currently not implemented) and IC formulation for Euler traveling vortex 106 // ***************************************************************************** 107 CEED_QFUNCTION_HELPER int Exact_ShockTube(CeedInt dim, CeedScalar time, 108 const CeedScalar X[], CeedInt Nf, CeedScalar q[], 109 void *ctx) { 110 111 // Context 112 const SetupContext context = (SetupContext)ctx; 113 const CeedScalar mid_point = context->mid_point; // Midpoint of the domain 114 const CeedScalar P_high = context->P_high; // Driver section pressure 115 const CeedScalar rho_high = context->rho_high; // Driver section density 116 const CeedScalar P_low = context->P_low; // Driven section pressure 117 const CeedScalar rho_low = context->rho_low; // Driven section density 118 119 // Setup 120 const CeedScalar gamma = 1.4; // ratio of specific heats 121 const CeedScalar x = X[0]; // Coordinates 122 123 CeedScalar rho, P, u[3] = {0.}; 124 125 // Initial Conditions 126 if (x <= mid_point) { 127 rho = rho_high; 128 P = P_high; 129 } else { 130 rho = rho_low; 131 P = P_low; 132 } 133 134 // Assign exact solution 135 q[0] = rho; 136 q[1] = rho * u[0]; 137 q[2] = rho * u[1]; 138 q[3] = rho * u[2]; 139 q[4] = P / (gamma-1.0) + rho * (u[0]*u[0]) / 2.; 140 141 // Return 142 return 0; 143 } 144 145 // ***************************************************************************** 146 // Helper function for computing flux Jacobian 147 // ***************************************************************************** 148 CEED_QFUNCTION_HELPER void ConvectiveFluxJacobian_Euler(CeedScalar dF[3][5][5], 149 const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 150 const CeedScalar gamma) { 151 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; // Velocity square 152 for (CeedInt i=0; i<3; i++) { // Jacobian matrices for 3 directions 153 for (CeedInt j=0; j<3; j++) { // Rows of each Jacobian matrix 154 dF[i][j+1][0] = ((i==j) ? ((gamma-1.)*(u_sq/2.)) : 0.) - u[i]*u[j]; 155 for (CeedInt k=0; k<3; k++) { // Columns of each Jacobian matrix 156 dF[i][0][k+1] = ((i==k) ? 1. : 0.); 157 dF[i][j+1][k+1] = ((j==k) ? u[i] : 0.) + 158 ((i==k) ? u[j] : 0.) - 159 ((i==j) ? u[k] : 0.) * (gamma-1.); 160 dF[i][4][k+1] = ((i==k) ? (E*gamma/rho - (gamma-1.)*u_sq/2.) : 0.) - 161 (gamma-1.)*u[i]*u[k]; 162 } 163 dF[i][j+1][4] = ((i==j) ? (gamma-1.) : 0.); 164 } 165 dF[i][4][0] = u[i] * ((gamma-1.)*u_sq - E*gamma/rho); 166 dF[i][4][4] = u[i] * gamma; 167 } 168 } 169 170 // ***************************************************************************** 171 // Helper function for calculating the covariant length scale in the direction 172 // of some 3 element input vector 173 // 174 // Where 175 // vec = vector that length is measured in the direction of 176 // h = covariant element length along vec 177 // ***************************************************************************** 178 CEED_QFUNCTION_HELPER CeedScalar Covariant_length_along_vector( 179 CeedScalar vec[3], const CeedScalar dXdx[3][3]) { 180 181 CeedScalar vec_norm = sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]); 182 CeedScalar vec_dot_jacobian[3] = {0.0}; 183 for (CeedInt i=0; i<3; i++) { 184 for (CeedInt j=0; j<3; j++) { 185 vec_dot_jacobian[i] += dXdx[j][i]*vec[i]; 186 } 187 } 188 CeedScalar norm_vec_dot_jacobian = sqrt(vec_dot_jacobian[0]*vec_dot_jacobian[0]+ 189 vec_dot_jacobian[1]*vec_dot_jacobian[1]+ 190 vec_dot_jacobian[2]*vec_dot_jacobian[2]); 191 CeedScalar h = 2.0 * vec_norm / norm_vec_dot_jacobian; 192 return h; 193 } 194 195 196 // ***************************************************************************** 197 // Helper function for computing Tau elements (stabilization constant) 198 // Model from: 199 // Stabilized Methods for Compressible Flows, Hughes et al 2010 200 // 201 // Spatial criterion #2 - Tau is a 3x3 diagonal matrix 202 // Tau[i] = c_tau h[i] Xi(Pe) / rho(A[i]) (no sum) 203 // 204 // Where 205 // c_tau = stabilization constant (0.5 is reported as "optimal") 206 // h[i] = 2 length(dxdX[i]) 207 // Pe = Peclet number ( Pe = sqrt(u u) / dot(dXdx,u) diffusivity ) 208 // Xi(Pe) = coth Pe - 1. / Pe (1. at large local Peclet number ) 209 // rho(A[i]) = spectral radius of the convective flux Jacobian i, 210 // wave speed in direction i 211 // ***************************************************************************** 212 CEED_QFUNCTION_HELPER void Tau_spatial(CeedScalar Tau_x[3], 213 const CeedScalar dXdx[3][3], const CeedScalar u[3], 214 const CeedScalar sound_speed, const CeedScalar c_tau) { 215 for (int i=0; i<3; i++) { 216 // length of element in direction i 217 CeedScalar h = 2 / sqrt(dXdx[0][i]*dXdx[0][i] + dXdx[1][i]*dXdx[1][i] + 218 dXdx[2][i]*dXdx[2][i]); 219 // fastest wave in direction i 220 CeedScalar fastest_wave = fabs(u[i]) + sound_speed; 221 Tau_x[i] = c_tau * h / fastest_wave; 222 } 223 } 224 225 // ***************************************************************************** 226 // This QFunction sets the initial conditions for shock tube 227 // ***************************************************************************** 228 CEED_QFUNCTION(ICsShockTube)(void *ctx, CeedInt Q, 229 const CeedScalar *const *in, CeedScalar *const *out) { 230 // Inputs 231 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 232 233 // Outputs 234 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 235 236 CeedPragmaSIMD 237 // Quadrature Point Loop 238 for (CeedInt i=0; i<Q; i++) { 239 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 240 CeedScalar q[5]; 241 242 Exact_ShockTube(3, 0., x, 5, q, ctx); 243 244 for (CeedInt j=0; j<5; j++) 245 q0[j][i] = q[j]; 246 } // End of Quadrature Point Loop 247 248 // Return 249 return 0; 250 } 251 252 // ***************************************************************************** 253 // This QFunction implements the following formulation of Euler equations 254 // with explicit time stepping method 255 // 256 // This is 3D Euler for compressible gas dynamics in conservation 257 // form with state variables of density, momentum density, and total 258 // energy density. 259 // 260 // State Variables: q = ( rho, U1, U2, U3, E ) 261 // rho - Mass Density 262 // Ui - Momentum Density, Ui = rho ui 263 // E - Total Energy Density, E = P / (gamma - 1) + rho (u u)/2 264 // 265 // Euler Equations: 266 // drho/dt + div( U ) = 0 267 // dU/dt + div( rho (u x u) + P I3 ) = 0 268 // dE/dt + div( (E + P) u ) = 0 269 // 270 // Equation of State: 271 // P = (gamma - 1) (E - rho (u u) / 2) 272 // 273 // Constants: 274 // cv , Specific heat, constant volume 275 // cp , Specific heat, constant pressure 276 // g , Gravity 277 // gamma = cp / cv, Specific heat ratio 278 // ***************************************************************************** 279 CEED_QFUNCTION(EulerShockTube)(void *ctx, CeedInt Q, 280 const CeedScalar *const *in, CeedScalar *const *out) { 281 // *INDENT-OFF* 282 // Inputs 283 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 284 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 285 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 286 // Outputs 287 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 288 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 289 290 const CeedScalar gamma = 1.4; 291 292 ShockTubeContext context = (ShockTubeContext)ctx; 293 const CeedScalar Cyzb = context->Cyzb; 294 const CeedScalar Byzb = context->Byzb; 295 const CeedScalar c_tau = context->c_tau; 296 297 CeedPragmaSIMD 298 // Quadrature Point Loop 299 for (CeedInt i=0; i<Q; i++) { 300 // *INDENT-OFF* 301 // Setup 302 // -- Interp in 303 const CeedScalar rho = q[0][i]; 304 const CeedScalar u[3] = {q[1][i] / rho, 305 q[2][i] / rho, 306 q[3][i] / rho 307 }; 308 const CeedScalar E = q[4][i]; 309 const CeedScalar drho[3] = {dq[0][0][i], 310 dq[1][0][i], 311 dq[2][0][i] 312 }; 313 const CeedScalar dU[3][3] = {{dq[0][1][i], 314 dq[1][1][i], 315 dq[2][1][i]}, 316 {dq[0][2][i], 317 dq[1][2][i], 318 dq[2][2][i]}, 319 {dq[0][3][i], 320 dq[1][3][i], 321 dq[2][3][i]} 322 }; 323 const CeedScalar dE[3] = {dq[0][4][i], 324 dq[1][4][i], 325 dq[2][4][i] 326 }; 327 // -- Interp-to-Interp q_data 328 const CeedScalar wdetJ = q_data[0][i]; 329 // -- Interp-to-Grad q_data 330 // ---- Inverse of change of coordinate matrix: X_i,j 331 // *INDENT-OFF* 332 const CeedScalar dXdx[3][3] = {{q_data[1][i], 333 q_data[2][i], 334 q_data[3][i]}, 335 {q_data[4][i], 336 q_data[5][i], 337 q_data[6][i]}, 338 {q_data[7][i], 339 q_data[8][i], 340 q_data[9][i]} 341 }; 342 // dU/dx 343 CeedScalar du[3][3] = {{0}}; 344 CeedScalar drhodx[3] = {0}; 345 CeedScalar dEdx[3] = {0}; 346 CeedScalar dUdx[3][3] = {{0}}; 347 CeedScalar dXdxdXdxT[3][3] = {{0}}; 348 for (int j=0; j<3; j++) { 349 for (int k=0; k<3; k++) { 350 du[j][k] = (dU[j][k] - drho[k]*u[j]) / rho; 351 drhodx[j] += drho[k] * dXdx[k][j]; 352 dEdx[j] += dE[k] * dXdx[k][j]; 353 for (int l=0; l<3; l++) { 354 dUdx[j][k] += dU[j][l] * dXdx[l][k]; 355 dXdxdXdxT[j][k] += dXdx[j][l]*dXdx[k][l]; //dXdx_j,k * dXdx_k,j 356 } 357 } 358 } 359 360 // *INDENT-ON* 361 const CeedScalar 362 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 363 E_internal = E - E_kinetic, 364 P = E_internal * (gamma - 1); // P = pressure 365 366 // The Physics 367 // Zero v and dv so all future terms can safely sum into it 368 for (int j=0; j<5; j++) { 369 v[j][i] = 0; 370 for (int k=0; k<3; k++) 371 dv[k][j][i] = 0; 372 } 373 374 // -- Density 375 // ---- u rho 376 for (int j=0; j<3; j++) 377 dv[j][0][i] += wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 378 rho*u[2]*dXdx[j][2]); 379 // -- Momentum 380 // ---- rho (u x u) + P I3 381 for (int j=0; j<3; j++) 382 for (int k=0; k<3; k++) 383 dv[k][j+1][i] += wdetJ*((rho*u[j]*u[0] + (j==0?P:0))*dXdx[k][0] + 384 (rho*u[j]*u[1] + (j==1?P:0))*dXdx[k][1] + 385 (rho*u[j]*u[2] + (j==2?P:0))*dXdx[k][2]); 386 // -- Total Energy Density 387 // ---- (E + P) u 388 for (int j=0; j<3; j++) 389 dv[j][4][i] += wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 390 u[2]*dXdx[j][2]); 391 392 // -- YZB stabilization 393 if (context->yzb) { 394 CeedScalar drho_norm = 0.0; // magnitude of the density gradient 395 CeedScalar j_vec[3] = {0.0}; // unit vector aligned with the density gradient 396 CeedScalar h_shock = 0.0; // element lengthscale 397 CeedScalar acoustic_vel = 0.0; // characteristic velocity, acoustic speed 398 CeedScalar tau_shock = 0.0; // timescale 399 CeedScalar nu_shock = 0.0; // artificial diffusion 400 401 // Unit vector aligned with the density gradient 402 drho_norm = sqrt(drhodx[0]*drhodx[0] + drhodx[1]*drhodx[1] + 403 drhodx[2]*drhodx[2]); 404 for (int j=0; j<3; j++) 405 j_vec[j] = drhodx[j] / (drho_norm + 1e-20); 406 407 if (drho_norm == 0.0) { 408 nu_shock = 0.0; 409 } else { 410 h_shock = Covariant_length_along_vector(j_vec, dXdx); 411 h_shock /= Cyzb; 412 acoustic_vel = sqrt(gamma*P/rho); 413 tau_shock = h_shock / (2*acoustic_vel) * pow(drho_norm * h_shock / rho, Byzb); 414 nu_shock = fabs(tau_shock * acoustic_vel * acoustic_vel); 415 } 416 417 for (int j=0; j<3; j++) 418 dv[j][0][i] -= wdetJ * nu_shock * drhodx[j]; 419 420 for (int k=0; k<3; k++) 421 for (int j=0; j<3; j++) 422 dv[j][k][i] -= wdetJ * nu_shock * du[k][j]; 423 424 for (int j=0; j<3; j++) 425 dv[j][4][i] -= wdetJ * nu_shock * dEdx[j]; 426 } 427 428 // Stabilization 429 // Need the Jacobian for the advective fluxes for stabilization 430 // indexed as: jacob_F_conv[direction][flux component][solution component] 431 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}}; 432 ConvectiveFluxJacobian_Euler(jacob_F_conv, rho, u, E, gamma); 433 434 435 // dqdx collects drhodx, dUdx and dEdx in one vector 436 CeedScalar dqdx[5][3]; 437 for (int j=0; j<3; j++) { 438 dqdx[0][j] = drhodx[j]; 439 dqdx[4][j] = dEdx[j]; 440 for (int k=0; k<3; k++) 441 dqdx[k+1][j] = dUdx[k][j]; 442 } 443 444 // strong_conv = dF/dq * dq/dx (Strong convection) 445 CeedScalar strong_conv[5] = {0}; 446 for (int j=0; j<3; j++) 447 for (int k=0; k<5; k++) 448 for (int l=0; l<5; l++) 449 strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; 450 451 // Stabilization 452 // -- Tau elements 453 const CeedScalar sound_speed = sqrt(gamma * P / rho); 454 CeedScalar Tau_x[3] = {0.}; 455 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau); 456 457 CeedScalar stab[5][3] = {0}; 458 switch (context->stabilization) { 459 case 0: // Galerkin 460 break; 461 case 1: // SU 462 for (int j=0; j<3; j++) 463 for (int k=0; k<5; k++) 464 for (int l=0; l<5; l++) { 465 stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l]; 466 } 467 for (int j=0; j<5; j++) 468 for (int k=0; k<3; k++) 469 dv[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 470 stab[j][1] * dXdx[k][1] + 471 stab[j][2] * dXdx[k][2]); 472 break; 473 } 474 475 } // End Quadrature Point Loop 476 477 // Return 478 return 0; 479 } 480 481 #endif // shocktube_h 482