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 14 // Finite Volume WENO Schemes, Zhang, Zhang, and Shu (2011). 15 16 #ifndef eulervortex_h 17 #define eulervortex_h 18 19 #include <math.h> 20 #include <ceed.h> 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 63 // (currently not implemented) and IC formulation for Euler traveling vortex 64 // ***************************************************************************** 65 CEED_QFUNCTION_HELPER int Exact_Euler(CeedInt dim, CeedScalar time, 66 const CeedScalar X[], CeedInt Nf, CeedScalar q[], 67 void *ctx) { 68 // Context 69 const EulerContext context = (EulerContext)ctx; 70 const CeedScalar vortex_strength = context->vortex_strength; 71 const CeedScalar *center = context->center; // Center of the domain 72 const CeedScalar *mean_velocity = context->mean_velocity; 73 74 // Setup 75 const CeedScalar gamma = 1.4; 76 const CeedScalar cv = 2.5; 77 const CeedScalar R = 1.; 78 const CeedScalar x = X[0], y = X[1]; // Coordinates 79 // Vortex center 80 const CeedScalar xc = center[0] + mean_velocity[0] * time; 81 const CeedScalar yc = center[1] + mean_velocity[1] * time; 82 83 const CeedScalar x0 = x - xc; 84 const CeedScalar y0 = y - yc; 85 const CeedScalar r = sqrt( x0*x0 + y0*y0 ); 86 const CeedScalar C = vortex_strength * exp((1. - r*r)/2.) / (2. * M_PI); 87 const CeedScalar delta_T = - (gamma - 1.) * vortex_strength * vortex_strength * 88 exp(1 - r*r) / (8. * gamma * M_PI * M_PI); 89 const CeedScalar S_vortex = 1; // no perturbation in the entropy P / rho^gamma 90 const CeedScalar S_bubble = (gamma - 1.) * vortex_strength * vortex_strength / 91 (8.*gamma*M_PI*M_PI); 92 CeedScalar rho, P, T, E, u[3] = {0.}; 93 94 // Initial Conditions 95 switch (context->euler_test) { 96 case 0: // Traveling vortex 97 T = 1 + delta_T; 98 // P = rho * T 99 // P = S * rho^gamma 100 // Solve for rho, then substitute for P 101 rho = pow(T/S_vortex, 1 / (gamma - 1.)); 102 P = rho * T; 103 u[0] = mean_velocity[0] - C*y0; 104 u[1] = mean_velocity[1] + C*x0; 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] = P / (gamma - 1.) + rho * (u[0]*u[0] + u[1]*u[1]) / 2.; 112 break; 113 case 1: // Constant zero velocity, density constant, total energy constant 114 rho = 1.; 115 E = 2.; 116 117 // Assign exact solution 118 q[0] = rho; 119 q[1] = rho * u[0]; 120 q[2] = rho * u[1]; 121 q[3] = rho * u[2]; 122 q[4] = E; 123 break; 124 case 2: // Constant nonzero velocity, density constant, total energy constant 125 rho = 1.; 126 E = 2.; 127 u[0] = mean_velocity[0]; 128 u[1] = mean_velocity[1]; 129 130 // Assign exact solution 131 q[0] = rho; 132 q[1] = rho * u[0]; 133 q[2] = rho * u[1]; 134 q[3] = rho * u[2]; 135 q[4] = E; 136 break; 137 case 3: // Velocity zero, pressure constant 138 // (so density and internal energy will be non-constant), 139 // but the velocity should stay zero and the bubble won't diffuse 140 // (for Euler, where there is no thermal conductivity) 141 P = 1.; 142 T = 1. - S_bubble * exp(1. - r*r); 143 rho = P / (R*T); 144 145 // Assign exact solution 146 q[0] = rho; 147 q[1] = rho * u[0]; 148 q[2] = rho * u[1]; 149 q[3] = rho * u[2]; 150 q[4] = rho * (cv * T + (u[0]*u[0] + u[1]*u[1])/2.); 151 break; 152 case 4: // Constant nonzero velocity, pressure constant 153 // (so density and internal energy will be non-constant), 154 // it should be transported across the domain, but velocity stays constant 155 P = 1.; 156 T = 1. - S_bubble * exp(1. - r*r); 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 case 5: // non-smooth thermal bubble - cylinder 169 P = 1.; 170 T = 1. - (r < 1. ? S_bubble : 0.); 171 rho = P / (R*T); 172 u[0] = mean_velocity[0]; 173 u[1] = mean_velocity[1]; 174 175 // Assign exact solution 176 q[0] = rho; 177 q[1] = rho * u[0]; 178 q[2] = rho * u[1]; 179 q[3] = rho * u[2]; 180 q[4] = rho * (cv * T + (u[0]*u[0] + u[1]*u[1])/2.); 181 break; 182 } 183 // Return 184 return 0; 185 } 186 187 // ***************************************************************************** 188 // Helper function for computing flux Jacobian 189 // ***************************************************************************** 190 CEED_QFUNCTION_HELPER void ConvectiveFluxJacobian_Euler(CeedScalar dF[3][5][5], 191 const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 192 const CeedScalar gamma) { 193 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; // Velocity square 194 for (CeedInt i=0; i<3; i++) { // Jacobian matrices for 3 directions 195 for (CeedInt j=0; j<3; j++) { // Rows of each Jacobian matrix 196 dF[i][j+1][0] = ((i==j) ? ((gamma-1.)*(u_sq/2.)) : 0.) - u[i]*u[j]; 197 for (CeedInt k=0; k<3; k++) { // Columns of each Jacobian matrix 198 dF[i][0][k+1] = ((i==k) ? 1. : 0.); 199 dF[i][j+1][k+1] = ((j==k) ? u[i] : 0.) + 200 ((i==k) ? u[j] : 0.) - 201 ((i==j) ? u[k] : 0.) * (gamma-1.); 202 dF[i][4][k+1] = ((i==k) ? (E*gamma/rho - (gamma-1.)*u_sq/2.) : 0.) - 203 (gamma-1.)*u[i]*u[k]; 204 } 205 dF[i][j+1][4] = ((i==j) ? (gamma-1.) : 0.); 206 } 207 dF[i][4][0] = u[i] * ((gamma-1.)*u_sq - E*gamma/rho); 208 dF[i][4][4] = u[i] * gamma; 209 } 210 } 211 212 // ***************************************************************************** 213 // Helper function for computing Tau elements (stabilization constant) 214 // Model from: 215 // Stabilized Methods for Compressible Flows, Hughes et al 2010 216 // 217 // Spatial criterion #2 - Tau is a 3x3 diagonal matrix 218 // Tau[i] = c_tau h[i] Xi(Pe) / rho(A[i]) (no sum) 219 // 220 // Where 221 // c_tau = stabilization constant (0.5 is reported as "optimal") 222 // h[i] = 2 length(dxdX[i]) 223 // Pe = Peclet number ( Pe = sqrt(u u) / dot(dXdx,u) diffusivity ) 224 // Xi(Pe) = coth Pe - 1. / Pe (1. at large local Peclet number ) 225 // rho(A[i]) = spectral radius of the convective flux Jacobian i, 226 // wave speed in direction i 227 // ***************************************************************************** 228 CEED_QFUNCTION_HELPER void Tau_spatial(CeedScalar Tau_x[3], 229 const CeedScalar dXdx[3][3], const CeedScalar u[3], 230 const CeedScalar sound_speed, const CeedScalar c_tau) { 231 for (CeedInt i=0; i<3; i++) { 232 // length of element in direction i 233 CeedScalar h = 2 / sqrt(dXdx[0][i]*dXdx[0][i] + dXdx[1][i]*dXdx[1][i] + 234 dXdx[2][i]*dXdx[2][i]); 235 // fastest wave in direction i 236 CeedScalar fastest_wave = fabs(u[i]) + sound_speed; 237 Tau_x[i] = c_tau * h / fastest_wave; 238 } 239 } 240 241 // ***************************************************************************** 242 // This QFunction sets the initial conditions for Euler traveling vortex 243 // ***************************************************************************** 244 CEED_QFUNCTION(ICsEuler)(void *ctx, CeedInt Q, 245 const CeedScalar *const *in, CeedScalar *const *out) { 246 // Inputs 247 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 248 249 // Outputs 250 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 251 const EulerContext context = (EulerContext)ctx; 252 253 CeedPragmaSIMD 254 // Quadrature Point Loop 255 for (CeedInt i=0; i<Q; i++) { 256 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 257 CeedScalar q[5] = {0.}; 258 259 Exact_Euler(3, context->curr_time, x, 5, q, ctx); 260 261 for (CeedInt j=0; j<5; j++) 262 q0[j][i] = q[j]; 263 } // End of Quadrature Point Loop 264 265 // Return 266 return 0; 267 } 268 269 // ***************************************************************************** 270 // This QFunction implements the following formulation of Euler equations 271 // with explicit time stepping method 272 // 273 // This is 3D Euler for compressible gas dynamics in conservation 274 // form with state variables of density, momentum density, and total 275 // energy density. 276 // 277 // State Variables: q = ( rho, U1, U2, U3, E ) 278 // rho - Mass Density 279 // Ui - Momentum Density, Ui = rho ui 280 // E - Total Energy Density, E = P / (gamma - 1) + rho (u u)/2 281 // 282 // Euler Equations: 283 // drho/dt + div( U ) = 0 284 // dU/dt + div( rho (u x u) + P I3 ) = 0 285 // dE/dt + div( (E + P) u ) = 0 286 // 287 // Equation of State: 288 // P = (gamma - 1) (E - rho (u u) / 2) 289 // 290 // Constants: 291 // cv , Specific heat, constant volume 292 // cp , Specific heat, constant pressure 293 // g , Gravity 294 // gamma = cp / cv, Specific heat ratio 295 // ***************************************************************************** 296 CEED_QFUNCTION(Euler)(void *ctx, CeedInt Q, 297 const CeedScalar *const *in, CeedScalar *const *out) { 298 // *INDENT-OFF* 299 // Inputs 300 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 301 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 302 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 303 // Outputs 304 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 305 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 306 307 EulerContext context = (EulerContext)ctx; 308 const CeedScalar c_tau = context->c_tau; 309 const CeedScalar gamma = 1.4; 310 311 CeedPragmaSIMD 312 // Quadrature Point Loop 313 for (CeedInt i=0; i<Q; i++) { 314 // *INDENT-OFF* 315 // Setup 316 // -- Interp in 317 const CeedScalar rho = q[0][i]; 318 const CeedScalar u[3] = {q[1][i] / rho, 319 q[2][i] / rho, 320 q[3][i] / rho 321 }; 322 const CeedScalar E = q[4][i]; 323 const CeedScalar drho[3] = {dq[0][0][i], 324 dq[1][0][i], 325 dq[2][0][i] 326 }; 327 const CeedScalar dU[3][3] = {{dq[0][1][i], 328 dq[1][1][i], 329 dq[2][1][i]}, 330 {dq[0][2][i], 331 dq[1][2][i], 332 dq[2][2][i]}, 333 {dq[0][3][i], 334 dq[1][3][i], 335 dq[2][3][i]} 336 }; 337 const CeedScalar dE[3] = {dq[0][4][i], 338 dq[1][4][i], 339 dq[2][4][i] 340 }; 341 // -- Interp-to-Interp q_data 342 const CeedScalar wdetJ = q_data[0][i]; 343 // -- Interp-to-Grad q_data 344 // ---- Inverse of change of coordinate matrix: X_i,j 345 // *INDENT-OFF* 346 const CeedScalar dXdx[3][3] = {{q_data[1][i], 347 q_data[2][i], 348 q_data[3][i]}, 349 {q_data[4][i], 350 q_data[5][i], 351 q_data[6][i]}, 352 {q_data[7][i], 353 q_data[8][i], 354 q_data[9][i]} 355 }; 356 // *INDENT-ON* 357 // dU/dx 358 CeedScalar drhodx[3] = {0.}; 359 CeedScalar dEdx[3] = {0.}; 360 CeedScalar dUdx[3][3] = {{0.}}; 361 CeedScalar dXdxdXdxT[3][3] = {{0.}}; 362 for (CeedInt j=0; j<3; j++) { 363 for (CeedInt k=0; k<3; k++) { 364 drhodx[j] += drho[k] * dXdx[k][j]; 365 dEdx[j] += dE[k] * dXdx[k][j]; 366 for (CeedInt l=0; l<3; l++) { 367 dUdx[j][k] += dU[j][l] * dXdx[l][k]; 368 dXdxdXdxT[j][k] += dXdx[j][l]*dXdx[k][l]; //dXdx_j,k * dXdx_k,j 369 } 370 } 371 } 372 // Pressure 373 const CeedScalar 374 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 375 E_internal = E - E_kinetic, 376 P = E_internal * (gamma - 1.); // P = pressure 377 378 // The Physics 379 // Zero v and dv so all future terms can safely sum into it 380 for (CeedInt j=0; j<5; j++) { 381 v[j][i] = 0.; 382 for (CeedInt k=0; k<3; k++) 383 dv[k][j][i] = 0.; 384 } 385 386 // -- Density 387 // ---- u rho 388 for (CeedInt j=0; j<3; j++) 389 dv[j][0][i] += wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 390 rho*u[2]*dXdx[j][2]); 391 // -- Momentum 392 // ---- rho (u x u) + P I3 393 for (CeedInt j=0; j<3; j++) 394 for (CeedInt k=0; k<3; k++) 395 dv[k][j+1][i] += wdetJ*((rho*u[j]*u[0] + (j==0?P:0.))*dXdx[k][0] + 396 (rho*u[j]*u[1] + (j==1?P:0.))*dXdx[k][1] + 397 (rho*u[j]*u[2] + (j==2?P:0.))*dXdx[k][2]); 398 // -- Total Energy Density 399 // ---- (E + P) u 400 for (CeedInt j=0; j<3; j++) 401 dv[j][4][i] += wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 402 u[2]*dXdx[j][2]); 403 404 // --Stabilization terms 405 // ---- jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 406 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}}; 407 ConvectiveFluxJacobian_Euler(jacob_F_conv, rho, u, E, gamma); 408 409 // ---- dqdx collects drhodx, dUdx and dEdx in one vector 410 CeedScalar dqdx[5][3]; 411 for (CeedInt j=0; j<3; j++) { 412 dqdx[0][j] = drhodx[j]; 413 dqdx[4][j] = dEdx[j]; 414 for (CeedInt k=0; k<3; k++) 415 dqdx[k+1][j] = dUdx[k][j]; 416 } 417 418 // ---- strong_conv = dF/dq * dq/dx (Strong convection) 419 CeedScalar strong_conv[5] = {0.}; 420 for (CeedInt j=0; j<3; j++) 421 for (CeedInt k=0; k<5; k++) 422 for (CeedInt l=0; l<5; l++) 423 strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; 424 425 // Stabilization 426 // -- Tau elements 427 const CeedScalar sound_speed = sqrt(gamma * P / rho); 428 CeedScalar Tau_x[3] = {0.}; 429 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau); 430 431 // -- Stabilization method: none or SU 432 CeedScalar stab[5][3] = {{0.}}; 433 switch (context->stabilization) { 434 case 0: // Galerkin 435 break; 436 case 1: // SU 437 for (CeedInt j=0; j<3; j++) 438 for (CeedInt k=0; k<5; k++) 439 for (CeedInt l=0; l<5; l++) 440 stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l]; 441 442 for (CeedInt j=0; j<5; j++) 443 for (CeedInt k=0; k<3; k++) 444 dv[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 445 stab[j][1] * dXdx[k][1] + 446 stab[j][2] * dXdx[k][2]); 447 break; 448 case 2: // SUPG is not implemented for explicit scheme 449 break; 450 } 451 452 } // End Quadrature Point Loop 453 454 // Return 455 return 0; 456 } 457 458 // ***************************************************************************** 459 // This QFunction implements the Euler equations with (mentioned above) 460 // with implicit time stepping method 461 // 462 // ***************************************************************************** 463 CEED_QFUNCTION(IFunction_Euler)(void *ctx, CeedInt Q, 464 const CeedScalar *const *in, CeedScalar *const *out) { 465 // *INDENT-OFF* 466 // Inputs 467 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 468 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 469 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 470 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 471 // Outputs 472 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 473 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 474 475 EulerContext context = (EulerContext)ctx; 476 const CeedScalar c_tau = context->c_tau; 477 const CeedScalar gamma = 1.4; 478 479 CeedPragmaSIMD 480 // Quadrature Point Loop 481 for (CeedInt i=0; i<Q; i++) { 482 // *INDENT-OFF* 483 // Setup 484 // -- Interp in 485 const CeedScalar rho = q[0][i]; 486 const CeedScalar u[3] = {q[1][i] / rho, 487 q[2][i] / rho, 488 q[3][i] / rho 489 }; 490 const CeedScalar E = q[4][i]; 491 const CeedScalar drho[3] = {dq[0][0][i], 492 dq[1][0][i], 493 dq[2][0][i] 494 }; 495 const CeedScalar dU[3][3] = {{dq[0][1][i], 496 dq[1][1][i], 497 dq[2][1][i]}, 498 {dq[0][2][i], 499 dq[1][2][i], 500 dq[2][2][i]}, 501 {dq[0][3][i], 502 dq[1][3][i], 503 dq[2][3][i]} 504 }; 505 const CeedScalar dE[3] = {dq[0][4][i], 506 dq[1][4][i], 507 dq[2][4][i] 508 }; 509 // -- Interp-to-Interp q_data 510 const CeedScalar wdetJ = q_data[0][i]; 511 // -- Interp-to-Grad q_data 512 // ---- Inverse of change of coordinate matrix: X_i,j 513 // *INDENT-OFF* 514 const CeedScalar dXdx[3][3] = {{q_data[1][i], 515 q_data[2][i], 516 q_data[3][i]}, 517 {q_data[4][i], 518 q_data[5][i], 519 q_data[6][i]}, 520 {q_data[7][i], 521 q_data[8][i], 522 q_data[9][i]} 523 }; 524 // *INDENT-ON* 525 // dU/dx 526 CeedScalar drhodx[3] = {0.}; 527 CeedScalar dEdx[3] = {0.}; 528 CeedScalar dUdx[3][3] = {{0.}}; 529 CeedScalar dXdxdXdxT[3][3] = {{0.}}; 530 for (CeedInt j=0; j<3; j++) { 531 for (CeedInt k=0; k<3; k++) { 532 drhodx[j] += drho[k] * dXdx[k][j]; 533 dEdx[j] += dE[k] * dXdx[k][j]; 534 for (CeedInt l=0; l<3; l++) { 535 dUdx[j][k] += dU[j][l] * dXdx[l][k]; 536 dXdxdXdxT[j][k] += dXdx[j][l]*dXdx[k][l]; //dXdx_j,k * dXdx_k,j 537 } 538 } 539 } 540 const CeedScalar 541 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 542 E_internal = E - E_kinetic, 543 P = E_internal * (gamma - 1.); // P = pressure 544 545 // The Physics 546 // Zero v and dv so all future terms can safely sum into it 547 for (CeedInt j=0; j<5; j++) { 548 v[j][i] = 0.; 549 for (CeedInt k=0; k<3; k++) 550 dv[k][j][i] = 0.; 551 } 552 //-----mass matrix 553 for (CeedInt j=0; j<5; j++) 554 v[j][i] += wdetJ*q_dot[j][i]; 555 556 // -- Density 557 // ---- u rho 558 for (CeedInt j=0; j<3; j++) 559 dv[j][0][i] -= wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 560 rho*u[2]*dXdx[j][2]); 561 // -- Momentum 562 // ---- rho (u x u) + P I3 563 for (CeedInt j=0; j<3; j++) 564 for (CeedInt k=0; k<3; k++) 565 dv[k][j+1][i] -= wdetJ*((rho*u[j]*u[0] + (j==0?P:0.))*dXdx[k][0] + 566 (rho*u[j]*u[1] + (j==1?P:0.))*dXdx[k][1] + 567 (rho*u[j]*u[2] + (j==2?P:0.))*dXdx[k][2]); 568 // -- Total Energy Density 569 // ---- (E + P) u 570 for (CeedInt j=0; j<3; j++) 571 dv[j][4][i] -= wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 572 u[2]*dXdx[j][2]); 573 574 // -- Stabilization terms 575 // ---- jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 576 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}}; 577 ConvectiveFluxJacobian_Euler(jacob_F_conv, rho, u, E, gamma); 578 579 // ---- dqdx collects drhodx, dUdx and dEdx in one vector 580 CeedScalar dqdx[5][3]; 581 for (CeedInt j=0; j<3; j++) { 582 dqdx[0][j] = drhodx[j]; 583 dqdx[4][j] = dEdx[j]; 584 for (CeedInt k=0; k<3; k++) 585 dqdx[k+1][j] = dUdx[k][j]; 586 } 587 588 // ---- strong_conv = dF/dq * dq/dx (Strong convection) 589 CeedScalar strong_conv[5] = {0.}; 590 for (CeedInt j=0; j<3; j++) 591 for (CeedInt k=0; k<5; k++) 592 for (CeedInt l=0; l<5; l++) 593 strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; 594 595 // ---- Strong residual 596 CeedScalar strong_res[5]; 597 for (CeedInt j=0; j<5; j++) 598 strong_res[j] = q_dot[j][i] + strong_conv[j]; 599 600 // Stabilization 601 // -- Tau elements 602 const CeedScalar sound_speed = sqrt(gamma * P / rho); 603 CeedScalar Tau_x[3] = {0.}; 604 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau); 605 606 // -- Stabilization method: none, SU, or SUPG 607 CeedScalar stab[5][3] = {{0.}}; 608 switch (context->stabilization) { 609 case 0: // Galerkin 610 break; 611 case 1: // SU 612 for (CeedInt j=0; j<3; j++) 613 for (CeedInt k=0; k<5; k++) 614 for (CeedInt l=0; l<5; l++) 615 stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l]; 616 617 for (CeedInt j=0; j<5; j++) 618 for (CeedInt k=0; k<3; k++) 619 dv[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 620 stab[j][1] * dXdx[k][1] + 621 stab[j][2] * dXdx[k][2]); 622 break; 623 case 2: // SUPG 624 for (CeedInt j=0; j<3; j++) 625 for (CeedInt k=0; k<5; k++) 626 for (CeedInt l=0; l<5; l++) 627 stab[k][j] = jacob_F_conv[j][k][l] * Tau_x[j] * strong_res[l]; 628 629 for (CeedInt j=0; j<5; j++) 630 for (CeedInt k=0; k<3; k++) 631 dv[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 632 stab[j][1] * dXdx[k][1] + 633 stab[j][2] * dXdx[k][2]); 634 break; 635 } 636 } // End Quadrature Point Loop 637 638 // Return 639 return 0; 640 } 641 // ***************************************************************************** 642 // This QFunction sets the inflow boundary conditions for 643 // the traveling vortex problem. 644 // 645 // Prescribed T_inlet and P_inlet are converted to conservative variables 646 // and applied weakly. 647 // 648 // ***************************************************************************** 649 CEED_QFUNCTION(TravelingVortex_Inflow)(void *ctx, CeedInt Q, 650 const CeedScalar *const *in, 651 CeedScalar *const *out) { 652 // *INDENT-OFF* 653 // Inputs 654 const CeedScalar (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 655 // Outputs 656 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 657 // *INDENT-ON* 658 EulerContext context = (EulerContext)ctx; 659 const int euler_test = context->euler_test; 660 const bool implicit = context->implicit; 661 CeedScalar *mean_velocity = context->mean_velocity; 662 const CeedScalar cv = 2.5; 663 const CeedScalar R = 1.; 664 CeedScalar T_inlet; 665 CeedScalar P_inlet; 666 667 // For test cases 1 and 3 the background velocity is zero 668 if (euler_test == 1 || euler_test == 3) 669 for (CeedInt i=0; i<3; i++) mean_velocity[i] = 0.; 670 671 // For test cases 1 and 2, T_inlet = T_inlet = 0.4 672 if (euler_test == 1 || euler_test == 2) T_inlet = P_inlet = .4; 673 else T_inlet = P_inlet = 1.; 674 675 CeedPragmaSIMD 676 // Quadrature Point Loop 677 for (CeedInt i=0; i<Q; i++) { 678 // Setup 679 // -- Interp-to-Interp q_data 680 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 681 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 682 // We can effect this by swapping the sign on this weight 683 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 684 // ---- Normal vect 685 const CeedScalar norm[3] = {q_data_sur[1][i], 686 q_data_sur[2][i], 687 q_data_sur[3][i] 688 }; 689 690 // face_normal = Normal vector of the face 691 const CeedScalar face_normal = norm[0]*mean_velocity[0] + 692 norm[1]*mean_velocity[1] + 693 norm[2]*mean_velocity[2]; 694 // The Physics 695 // Zero v so all future terms can safely sum into it 696 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 697 698 // Implementing in/outflow BCs 699 if (face_normal > 0) { 700 } else { // inflow 701 const CeedScalar rho_inlet = P_inlet/(R*T_inlet); 702 const CeedScalar E_kinetic_inlet = (mean_velocity[0]*mean_velocity[0] + 703 mean_velocity[1]*mean_velocity[1]) / 2.; 704 // incoming total energy 705 const CeedScalar E_inlet = rho_inlet * (cv * T_inlet + E_kinetic_inlet); 706 707 // The Physics 708 // -- Density 709 v[0][i] -= wdetJb * rho_inlet * face_normal; 710 711 // -- Momentum 712 for (CeedInt j=0; j<3; j++) 713 v[j+1][i] -= wdetJb *(rho_inlet * face_normal * mean_velocity[j] + 714 norm[j] * P_inlet); 715 716 // -- Total Energy Density 717 v[4][i] -= wdetJb * face_normal * (E_inlet + P_inlet); 718 } 719 720 } // End Quadrature Point Loop 721 return 0; 722 } 723 724 // ***************************************************************************** 725 // This QFunction sets the outflow boundary conditions for 726 // the Euler solver. 727 // 728 // Outflow BCs: 729 // The validity of the weak form of the governing equations is 730 // extended to the outflow. 731 // 732 // ***************************************************************************** 733 CEED_QFUNCTION(Euler_Outflow)(void *ctx, CeedInt Q, 734 const CeedScalar *const *in, 735 CeedScalar *const *out) { 736 // *INDENT-OFF* 737 // Inputs 738 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 739 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 740 // Outputs 741 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 742 // *INDENT-ON* 743 EulerContext context = (EulerContext)ctx; 744 const bool implicit = context->implicit; 745 CeedScalar *mean_velocity = context->mean_velocity; 746 747 const CeedScalar gamma = 1.4; 748 749 CeedPragmaSIMD 750 // Quadrature Point Loop 751 for (CeedInt i=0; i<Q; i++) { 752 // Setup 753 // -- Interp in 754 const CeedScalar rho = q[0][i]; 755 const CeedScalar u[3] = {q[1][i] / rho, 756 q[2][i] / rho, 757 q[3][i] / rho 758 }; 759 const CeedScalar E = q[4][i]; 760 761 // -- Interp-to-Interp q_data 762 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 763 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 764 // We can effect this by swapping the sign on this weight 765 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 766 // ---- Normal vectors 767 const CeedScalar norm[3] = {q_data_sur[1][i], 768 q_data_sur[2][i], 769 q_data_sur[3][i] 770 }; 771 772 // face_normal = Normal vector of the face 773 const CeedScalar face_normal = norm[0]*mean_velocity[0] + 774 norm[1]*mean_velocity[1] + 775 norm[2]*mean_velocity[2]; 776 // The Physics 777 // Zero v so all future terms can safely sum into it 778 for (CeedInt j=0; j<5; j++) v[j][i] = 0; 779 780 // Implementing in/outflow BCs 781 if (face_normal > 0) { // outflow 782 const CeedScalar E_kinetic = (u[0]*u[0] + u[1]*u[1]) / 2.; 783 const CeedScalar P = (E - E_kinetic * rho) * (gamma - 1.); // pressure 784 const CeedScalar u_normal = norm[0]*u[0] + norm[1]*u[1] + 785 norm[2]*u[2]; // Normal velocity 786 // The Physics 787 // -- Density 788 v[0][i] -= wdetJb * rho * u_normal; 789 790 // -- Momentum 791 for (CeedInt j=0; j<3; j++) 792 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + norm[j] * P); 793 794 // -- Total Energy Density 795 v[4][i] -= wdetJb * u_normal * (E + P); 796 } 797 } // End Quadrature Point Loop 798 return 0; 799 } 800 801 // ***************************************************************************** 802 803 #endif // eulervortex_h 804