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