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