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