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