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