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