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