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