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