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