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