1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 /// @file 5 /// Euler traveling vortex initial condition and operator for HONEE 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/types.h> 10 #ifndef CEED_RUNNING_JIT_PASS 11 #include <stdbool.h> 12 #endif 13 14 #include "utils.h" 15 16 typedef struct EulerContext_ *EulerContext; 17 struct EulerContext_ { 18 CeedScalar center[3]; 19 CeedScalar curr_time; 20 CeedScalar vortex_strength; 21 CeedScalar c_tau; 22 CeedScalar mean_velocity[3]; 23 bool implicit; 24 int euler_test; 25 int stabilization; // See StabilizationType: 0=none, 1=SU, 2=SUPG 26 }; 27 28 // ***************************************************************************** 29 // This function sets the initial conditions 30 // 31 // Temperature: 32 // T = 1 - (gamma - 1) vortex_strength**2 exp(1 - r**2) / (8 gamma pi**2) 33 // Density: 34 // rho = (T/S_vortex)^(1 / (gamma - 1)) 35 // Pressure: 36 // P = rho * T 37 // Velocity: 38 // ui = 1 + vortex_strength exp((1 - r**2)/2.) [yc - y, x - xc] / (2 pi) 39 // r = sqrt( (x - xc)**2 + (y - yc)**2 ) 40 // Velocity/Momentum Density: 41 // Ui = rho ui 42 // Total Energy: 43 // E = P / (gamma - 1) + rho (u u)/2 44 // 45 // Constants: 46 // cv , Specific heat, constant volume 47 // cp , Specific heat, constant pressure 48 // vortex_strength , Strength of vortex 49 // center , Location of bubble center 50 // gamma = cp / cv, Specific heat ratio 51 // 52 // ***************************************************************************** 53 54 // ***************************************************************************** 55 // This helper function provides support for the exact, time-dependent solution (currently not implemented) and IC formulation for Euler traveling 56 // vortex 57 // ***************************************************************************** 58 CEED_QFUNCTION_HELPER int Exact_Euler(CeedInt dim, CeedScalar time, const CeedScalar X[], CeedInt Nf, CeedScalar q[], void *ctx) { 59 // Context 60 const EulerContext context = (EulerContext)ctx; 61 const CeedScalar vortex_strength = context->vortex_strength; 62 const CeedScalar *center = context->center; // Center of the domain 63 const CeedScalar *mean_velocity = context->mean_velocity; 64 65 // Setup 66 const CeedScalar gamma = 1.4; 67 const CeedScalar cv = 2.5; 68 const CeedScalar R = 1.; 69 const CeedScalar x = X[0], y = X[1]; // Coordinates 70 // Vortex center 71 const CeedScalar xc = center[0] + mean_velocity[0] * time; 72 const CeedScalar yc = center[1] + mean_velocity[1] * time; 73 74 const CeedScalar x0 = x - xc; 75 const CeedScalar y0 = y - yc; 76 const CeedScalar r = sqrt(Square(x0) + Square(y0)); 77 const CeedScalar C = vortex_strength * exp((1. - r * r) / 2.) / (2. * M_PI); 78 const CeedScalar delta_T = -(gamma - 1.) * vortex_strength * vortex_strength * exp(1 - r * r) / (8. * gamma * M_PI * M_PI); 79 const CeedScalar S_vortex = 1; // no perturbation in the entropy P / rho^gamma 80 const CeedScalar S_bubble = (gamma - 1.) * vortex_strength * vortex_strength / (8. * gamma * M_PI * M_PI); 81 CeedScalar rho, P, T, E, u[3] = {0.}; 82 83 // Initial Conditions 84 switch (context->euler_test) { 85 case 0: // Traveling vortex 86 T = 1 + delta_T; 87 // P = rho * T 88 // P = S * rho^gamma 89 // Solve for rho, then substitute for P 90 rho = pow(T / S_vortex, 1 / (gamma - 1.)); 91 P = rho * T; 92 u[0] = mean_velocity[0] - C * y0; 93 u[1] = mean_velocity[1] + C * x0; 94 95 // Assign exact solution 96 q[0] = rho; 97 q[1] = rho * u[0]; 98 q[2] = rho * u[1]; 99 q[3] = rho * u[2]; 100 q[4] = P / (gamma - 1.) + rho * (u[0] * u[0] + u[1] * u[1]) / 2.; 101 break; 102 case 1: // Constant zero velocity, density constant, total energy constant 103 rho = 1.; 104 E = 2.; 105 106 // Assign exact solution 107 q[0] = rho; 108 q[1] = rho * u[0]; 109 q[2] = rho * u[1]; 110 q[3] = rho * u[2]; 111 q[4] = E; 112 break; 113 case 2: // Constant nonzero velocity, density constant, total energy constant 114 rho = 1.; 115 E = 2.; 116 u[0] = mean_velocity[0]; 117 u[1] = mean_velocity[1]; 118 119 // Assign exact solution 120 q[0] = rho; 121 q[1] = rho * u[0]; 122 q[2] = rho * u[1]; 123 q[3] = rho * u[2]; 124 q[4] = E; 125 break; 126 case 3: // Velocity zero, pressure constant (so density and internal energy will be non-constant), but the velocity should stay zero and the 127 // bubble won't diffuse 128 // (for Euler, where there is no thermal conductivity) 129 P = 1.; 130 T = 1. - S_bubble * exp(1. - r * r); 131 rho = P / (R * T); 132 133 // Assign exact solution 134 q[0] = rho; 135 q[1] = rho * u[0]; 136 q[2] = rho * u[1]; 137 q[3] = rho * u[2]; 138 q[4] = rho * (cv * T + (u[0] * u[0] + u[1] * u[1]) / 2.); 139 break; 140 case 4: // Constant nonzero velocity, pressure constant (so density and internal energy will be non-constant), 141 // It should be transported across the domain, but velocity stays constant 142 P = 1.; 143 T = 1. - S_bubble * exp(1. - r * r); 144 rho = P / (R * T); 145 u[0] = mean_velocity[0]; 146 u[1] = mean_velocity[1]; 147 148 // Assign exact solution 149 q[0] = rho; 150 q[1] = rho * u[0]; 151 q[2] = rho * u[1]; 152 q[3] = rho * u[2]; 153 q[4] = rho * (cv * T + (u[0] * u[0] + u[1] * u[1]) / 2.); 154 break; 155 case 5: // non-smooth thermal bubble - cylinder 156 P = 1.; 157 T = 1. - (r < 1. ? S_bubble : 0.); 158 rho = P / (R * T); 159 u[0] = mean_velocity[0]; 160 u[1] = mean_velocity[1]; 161 162 // Assign exact solution 163 q[0] = rho; 164 q[1] = rho * u[0]; 165 q[2] = rho * u[1]; 166 q[3] = rho * u[2]; 167 q[4] = rho * (cv * T + (u[0] * u[0] + u[1] * u[1]) / 2.); 168 break; 169 } 170 return 0; 171 } 172 173 // ***************************************************************************** 174 // Helper function for computing flux Jacobian 175 // ***************************************************************************** 176 CEED_QFUNCTION_HELPER void ConvectiveFluxJacobian_Euler(CeedScalar dF[3][5][5], const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 177 const CeedScalar gamma) { 178 CeedScalar u_sq = u[0] * u[0] + u[1] * u[1] + u[2] * u[2]; // Velocity square 179 for (CeedInt i = 0; i < 3; i++) { // Jacobian matrices for 3 directions 180 for (CeedInt j = 0; j < 3; j++) { // Rows of each Jacobian matrix 181 dF[i][j + 1][0] = ((i == j) ? ((gamma - 1.) * (u_sq / 2.)) : 0.) - u[i] * u[j]; 182 for (CeedInt k = 0; k < 3; k++) { // Columns of each Jacobian matrix 183 dF[i][0][k + 1] = ((i == k) ? 1. : 0.); 184 dF[i][j + 1][k + 1] = ((j == k) ? u[i] : 0.) + ((i == k) ? u[j] : 0.) - ((i == j) ? u[k] : 0.) * (gamma - 1.); 185 dF[i][4][k + 1] = ((i == k) ? (E * gamma / rho - (gamma - 1.) * u_sq / 2.) : 0.) - (gamma - 1.) * u[i] * u[k]; 186 } 187 dF[i][j + 1][4] = ((i == j) ? (gamma - 1.) : 0.); 188 } 189 dF[i][4][0] = u[i] * ((gamma - 1.) * u_sq - E * gamma / rho); 190 dF[i][4][4] = u[i] * gamma; 191 } 192 } 193 194 // ***************************************************************************** 195 // Helper function for computing Tau elements (stabilization constant) 196 // Model from: 197 // Stabilized Methods for Compressible Flows, Hughes et al 2010 198 // 199 // Spatial criterion #2 - Tau is a 3x3 diagonal matrix 200 // Tau[i] = c_tau h[i] Xi(Pe) / rho(A[i]) (no sum) 201 // 202 // Where 203 // c_tau = stabilization constant (0.5 is reported as "optimal") 204 // h[i] = 2 length(dxdX[i]) 205 // Pe = Peclet number ( Pe = sqrt(u u) / dot(dXdx,u) diffusivity ) 206 // Xi(Pe) = coth Pe - 1. / Pe (1. at large local Peclet number ) 207 // rho(A[i]) = spectral radius of the convective flux Jacobian i, wave speed in direction i 208 // ***************************************************************************** 209 CEED_QFUNCTION_HELPER void Tau_spatial(CeedScalar Tau_x[3], const CeedScalar dXdx[3][3], const CeedScalar u[3], const CeedScalar sound_speed, 210 const CeedScalar c_tau) { 211 for (CeedInt i = 0; i < 3; i++) { 212 // length of element in direction i 213 CeedScalar h = 2 / sqrt(Square(dXdx[0][i]) + Square(dXdx[1][i]) + Square(dXdx[2][i])); 214 // fastest wave in direction i 215 CeedScalar fastest_wave = fabs(u[i]) + sound_speed; 216 Tau_x[i] = c_tau * h / fastest_wave; 217 } 218 } 219 220 // ***************************************************************************** 221 // This QFunction sets the initial conditions for Euler traveling vortex 222 // ***************************************************************************** 223 CEED_QFUNCTION(ICsEuler)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 224 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 225 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 226 227 const EulerContext context = (EulerContext)ctx; 228 229 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 230 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 231 CeedScalar q[5] = {0.}; 232 233 Exact_Euler(3, context->curr_time, x, 5, q, ctx); 234 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 235 } 236 return 0; 237 } 238 239 // ***************************************************************************** 240 // This QFunction implements the following formulation of Euler equations with explicit time stepping method 241 // 242 // This is 3D Euler for compressible gas dynamics in conservation form with state variables of density, momentum density, and total energy density. 243 // 244 // State Variables: q = ( rho, U1, U2, U3, E ) 245 // rho - Mass Density 246 // Ui - Momentum Density, Ui = rho ui 247 // E - Total Energy Density, E = P / (gamma - 1) + rho (u u)/2 248 // 249 // Euler Equations: 250 // drho/dt + div( U ) = 0 251 // dU/dt + div( rho (u x u) + P I3 ) = 0 252 // dE/dt + div( (E + P) u ) = 0 253 // 254 // Equation of State: 255 // P = (gamma - 1) (E - rho (u u) / 2) 256 // 257 // Constants: 258 // cv , Specific heat, constant volume 259 // cp , Specific heat, constant pressure 260 // g , Gravity 261 // gamma = cp / cv, Specific heat ratio 262 // ***************************************************************************** 263 CEED_QFUNCTION(Euler)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 264 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 265 const CeedScalar(*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 266 const CeedScalar(*q_data) = in[2]; 267 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 268 CeedScalar(*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 269 270 EulerContext context = (EulerContext)ctx; 271 const CeedScalar c_tau = context->c_tau; 272 const CeedScalar gamma = 1.4; 273 274 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 275 // Setup 276 // -- Interp in 277 const CeedScalar rho = q[0][i]; 278 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 279 const CeedScalar E = q[4][i]; 280 const CeedScalar drho[3] = {dq[0][0][i], dq[1][0][i], dq[2][0][i]}; 281 const CeedScalar dU[3][3] = { 282 {dq[0][1][i], dq[1][1][i], dq[2][1][i]}, 283 {dq[0][2][i], dq[1][2][i], dq[2][2][i]}, 284 {dq[0][3][i], dq[1][3][i], dq[2][3][i]} 285 }; 286 const CeedScalar dE[3] = {dq[0][4][i], dq[1][4][i], dq[2][4][i]}; 287 CeedScalar wdetJ, dXdx[3][3]; 288 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 289 // dU/dx 290 CeedScalar drhodx[3] = {0.}; 291 CeedScalar dEdx[3] = {0.}; 292 CeedScalar dUdx[3][3] = {{0.}}; 293 CeedScalar dXdxdXdxT[3][3] = {{0.}}; 294 for (CeedInt j = 0; j < 3; j++) { 295 for (CeedInt k = 0; k < 3; k++) { 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 // Pressure 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 // --Stabilization terms 331 // ---- jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 332 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}}; 333 ConvectiveFluxJacobian_Euler(jacob_F_conv, rho, u, E, gamma); 334 335 // ---- dqdx collects drhodx, dUdx and dEdx in one vector 336 CeedScalar dqdx[5][3]; 337 for (CeedInt j = 0; j < 3; j++) { 338 dqdx[0][j] = drhodx[j]; 339 dqdx[4][j] = dEdx[j]; 340 for (CeedInt k = 0; k < 3; k++) dqdx[k + 1][j] = dUdx[k][j]; 341 } 342 343 // ---- strong_conv = dF/dq * dq/dx (Strong convection) 344 CeedScalar strong_conv[5] = {0.}; 345 for (CeedInt j = 0; j < 3; j++) { 346 for (CeedInt k = 0; k < 5; k++) { 347 for (CeedInt l = 0; l < 5; l++) strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; 348 } 349 } 350 351 // Stabilization 352 // -- Tau elements 353 const CeedScalar sound_speed = sqrt(gamma * P / rho); 354 CeedScalar Tau_x[3] = {0.}; 355 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau); 356 357 // -- Stabilization method: none or SU 358 CeedScalar stab[5][3] = {{0.}}; 359 switch (context->stabilization) { 360 case 0: // Galerkin 361 break; 362 case 1: // SU 363 for (CeedInt j = 0; j < 3; j++) { 364 for (CeedInt k = 0; k < 5; k++) { 365 for (CeedInt l = 0; l < 5; l++) stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l]; 366 } 367 } 368 369 for (CeedInt j = 0; j < 5; j++) { 370 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]); 371 } 372 break; 373 case 2: // SUPG is not implemented for explicit scheme 374 break; 375 } 376 } 377 return 0; 378 } 379 380 // ***************************************************************************** 381 // This QFunction implements the Euler equations with (mentioned above) with implicit time stepping method 382 // ***************************************************************************** 383 CEED_QFUNCTION(IFunction_Euler)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 384 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 385 const CeedScalar(*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 386 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 387 const CeedScalar(*q_data) = in[3]; 388 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 389 CeedScalar(*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 390 391 EulerContext context = (EulerContext)ctx; 392 const CeedScalar c_tau = context->c_tau; 393 const CeedScalar gamma = 1.4; 394 395 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 396 // Setup 397 // -- Interp in 398 const CeedScalar rho = q[0][i]; 399 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 400 const CeedScalar E = q[4][i]; 401 const CeedScalar drho[3] = {dq[0][0][i], dq[1][0][i], dq[2][0][i]}; 402 const CeedScalar dU[3][3] = { 403 {dq[0][1][i], dq[1][1][i], dq[2][1][i]}, 404 {dq[0][2][i], dq[1][2][i], dq[2][2][i]}, 405 {dq[0][3][i], dq[1][3][i], dq[2][3][i]} 406 }; 407 const CeedScalar dE[3] = {dq[0][4][i], dq[1][4][i], dq[2][4][i]}; 408 CeedScalar wdetJ, dXdx[3][3]; 409 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 410 // dU/dx 411 CeedScalar drhodx[3] = {0.}; 412 CeedScalar dEdx[3] = {0.}; 413 CeedScalar dUdx[3][3] = {{0.}}; 414 CeedScalar dXdxdXdxT[3][3] = {{0.}}; 415 for (CeedInt j = 0; j < 3; j++) { 416 for (CeedInt k = 0; k < 3; k++) { 417 drhodx[j] += drho[k] * dXdx[k][j]; 418 dEdx[j] += dE[k] * dXdx[k][j]; 419 for (CeedInt l = 0; l < 3; l++) { 420 dUdx[j][k] += dU[j][l] * dXdx[l][k]; 421 dXdxdXdxT[j][k] += dXdx[j][l] * dXdx[k][l]; // dXdx_j,k * dXdx_k,j 422 } 423 } 424 } 425 const CeedScalar E_kinetic = 0.5 * rho * (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]), E_internal = E - E_kinetic, 426 P = E_internal * (gamma - 1.); // P = pressure 427 428 // The Physics 429 // Zero v and dv so all future terms can safely sum into it 430 for (CeedInt j = 0; j < 5; j++) { 431 v[j][i] = 0.; 432 for (CeedInt k = 0; k < 3; k++) dv[k][j][i] = 0.; 433 } 434 //-----mass matrix 435 for (CeedInt j = 0; j < 5; j++) v[j][i] += wdetJ * q_dot[j][i]; 436 437 // -- Density 438 // ---- u rho 439 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]); 440 // -- Momentum 441 // ---- rho (u x u) + P I3 442 for (CeedInt j = 0; j < 3; j++) { 443 for (CeedInt k = 0; k < 3; k++) { 444 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] + 445 (rho * u[j] * u[2] + (j == 2 ? P : 0.)) * dXdx[k][2]); 446 } 447 } 448 // -- Total Energy Density 449 // ---- (E + P) u 450 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]); 451 452 // -- Stabilization terms 453 // ---- jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 454 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}}; 455 ConvectiveFluxJacobian_Euler(jacob_F_conv, rho, u, E, gamma); 456 457 // ---- dqdx collects drhodx, dUdx and dEdx in one vector 458 CeedScalar dqdx[5][3]; 459 for (CeedInt j = 0; j < 3; j++) { 460 dqdx[0][j] = drhodx[j]; 461 dqdx[4][j] = dEdx[j]; 462 for (CeedInt k = 0; k < 3; k++) dqdx[k + 1][j] = dUdx[k][j]; 463 } 464 465 // ---- strong_conv = dF/dq * dq/dx (Strong convection) 466 CeedScalar strong_conv[5] = {0.}; 467 for (CeedInt j = 0; j < 3; j++) { 468 for (CeedInt k = 0; k < 5; k++) { 469 for (CeedInt l = 0; l < 5; l++) strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; 470 } 471 } 472 473 // ---- Strong residual 474 CeedScalar strong_res[5]; 475 for (CeedInt j = 0; j < 5; j++) strong_res[j] = q_dot[j][i] + strong_conv[j]; 476 477 // Stabilization 478 // -- Tau elements 479 const CeedScalar sound_speed = sqrt(gamma * P / rho); 480 CeedScalar Tau_x[3] = {0.}; 481 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau); 482 483 // -- Stabilization method: none, SU, or SUPG 484 CeedScalar stab[5][3] = {{0.}}; 485 switch (context->stabilization) { 486 case 0: // Galerkin 487 break; 488 case 1: // SU 489 for (CeedInt j = 0; j < 3; j++) { 490 for (CeedInt k = 0; k < 5; k++) { 491 for (CeedInt l = 0; l < 5; l++) stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l]; 492 } 493 } 494 495 for (CeedInt j = 0; j < 5; j++) { 496 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]); 497 } 498 break; 499 case 2: // SUPG 500 for (CeedInt j = 0; j < 3; j++) { 501 for (CeedInt k = 0; k < 5; k++) { 502 for (CeedInt l = 0; l < 5; l++) stab[k][j] = jacob_F_conv[j][k][l] * Tau_x[j] * strong_res[l]; 503 } 504 } 505 506 for (CeedInt j = 0; j < 5; j++) { 507 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]); 508 } 509 break; 510 } 511 } 512 return 0; 513 } 514 // ***************************************************************************** 515 // This QFunction sets the inflow boundary conditions for the traveling vortex problem. 516 // 517 // Prescribed T_inlet and P_inlet are converted to conservative variables and applied weakly. 518 // ***************************************************************************** 519 CEED_QFUNCTION(TravelingVortex_Inflow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 520 const CeedScalar(*q_data_sur) = in[2]; 521 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 522 523 EulerContext context = (EulerContext)ctx; 524 const int euler_test = context->euler_test; 525 const bool is_implicit = context->implicit; 526 CeedScalar *mean_velocity = context->mean_velocity; 527 const CeedScalar cv = 2.5; 528 const CeedScalar R = 1.; 529 CeedScalar T_inlet; 530 CeedScalar P_inlet; 531 532 // For test cases 1 and 3 the background velocity is zero 533 if (euler_test == 1 || euler_test == 3) { 534 for (CeedInt i = 0; i < 3; i++) mean_velocity[i] = 0.; 535 } 536 537 // For test cases 1 and 2, T_inlet = T_inlet = 0.4 538 if (euler_test == 1 || euler_test == 2) T_inlet = P_inlet = .4; 539 else T_inlet = P_inlet = 1.; 540 541 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 542 CeedScalar wdetJb, normal[3]; 543 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal); 544 wdetJb *= is_implicit ? -1. : 1.; 545 546 // face_normal = Normal vector of the face 547 const CeedScalar face_normal = Dot3(normal, mean_velocity); 548 // The Physics 549 // Zero v so all future terms can safely sum into it 550 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0.; 551 552 // Implementing in/outflow BCs 553 if (face_normal > 0) { 554 } else { // inflow 555 const CeedScalar rho_inlet = P_inlet / (R * T_inlet); 556 const CeedScalar E_kinetic_inlet = (mean_velocity[0] * mean_velocity[0] + mean_velocity[1] * mean_velocity[1]) / 2.; 557 // incoming total energy 558 const CeedScalar E_inlet = rho_inlet * (cv * T_inlet + E_kinetic_inlet); 559 560 // The Physics 561 // -- Density 562 v[0][i] -= wdetJb * rho_inlet * face_normal; 563 564 // -- Momentum 565 for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho_inlet * face_normal * mean_velocity[j] + normal[j] * P_inlet); 566 567 // -- Total Energy Density 568 v[4][i] -= wdetJb * face_normal * (E_inlet + P_inlet); 569 } 570 } 571 return 0; 572 } 573 574 // ***************************************************************************** 575 // This QFunction sets the outflow boundary conditions for the Euler solver. 576 // 577 // Outflow BCs: 578 // The validity of the weak form of the governing equations is extended to the outflow. 579 // ***************************************************************************** 580 CEED_QFUNCTION(Euler_Outflow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 581 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 582 const CeedScalar(*q_data_sur) = in[2]; 583 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 584 585 EulerContext context = (EulerContext)ctx; 586 const bool is_implicit = context->implicit; 587 CeedScalar *mean_velocity = context->mean_velocity; 588 589 const CeedScalar gamma = 1.4; 590 591 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 592 // Setup 593 // -- Interp in 594 const CeedScalar rho = q[0][i]; 595 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 596 const CeedScalar E = q[4][i]; 597 598 CeedScalar wdetJb, normal[3]; 599 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal); 600 wdetJb *= is_implicit ? -1. : 1.; 601 602 // face_normal = Normal vector of the face 603 const CeedScalar face_normal = Dot3(normal, mean_velocity); 604 // The Physics 605 // Zero v so all future terms can safely sum into it 606 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0; 607 608 // Implementing in/outflow BCs 609 if (face_normal > 0) { // outflow 610 const CeedScalar E_kinetic = (u[0] * u[0] + u[1] * u[1]) / 2.; 611 const CeedScalar P = (E - E_kinetic * rho) * (gamma - 1.); // pressure 612 const CeedScalar u_normal = Dot3(normal, u); // Normal velocity 613 // The Physics 614 // -- Density 615 v[0][i] -= wdetJb * rho * u_normal; 616 617 // -- Momentum 618 for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho * u_normal * u[j] + normal[j] * P); 619 620 // -- Total Energy Density 621 v[4][i] -= wdetJb * u_normal * (E + P); 622 } 623 } 624 return 0; 625 } 626