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