1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Euler traveling vortex initial condition and operator for Navier-Stokes 19 /// example using PETSc 20 21 // Model from: 22 // On the Order of Accuracy and Numerical Performance of Two Classes of 23 // Finite Volume WENO Schemes, Zhang, Zhang, and Shu (2011). 24 25 #ifndef eulervortex_h 26 #define eulervortex_h 27 28 #include <math.h> 29 30 #ifndef M_PI 31 #define M_PI 3.14159265358979323846 32 #endif 33 34 #ifndef euler_context_struct 35 #define euler_context_struct 36 typedef struct EulerContext_ *EulerContext; 37 struct EulerContext_ { 38 CeedScalar center[3]; 39 CeedScalar curr_time; 40 CeedScalar vortex_strength; 41 CeedScalar c_tau; 42 CeedScalar mean_velocity[3]; 43 bool implicit; 44 int euler_test; 45 int stabilization; // See StabilizationType: 0=none, 1=SU, 2=SUPG 46 }; 47 #endif 48 49 // ***************************************************************************** 50 // This function sets the initial conditions 51 // 52 // Temperature: 53 // T = 1 - (gamma - 1) vortex_strength**2 exp(1 - r**2) / (8 gamma pi**2) 54 // Density: 55 // rho = (T/S_vortex)^(1 / (gamma - 1)) 56 // Pressure: 57 // P = rho * T 58 // Velocity: 59 // ui = 1 + vortex_strength exp((1 - r**2)/2.) [yc - y, x - xc] / (2 pi) 60 // r = sqrt( (x - xc)**2 + (y - yc)**2 ) 61 // Velocity/Momentum Density: 62 // Ui = rho ui 63 // Total Energy: 64 // E = P / (gamma - 1) + rho (u u)/2 65 // 66 // Constants: 67 // cv , Specific heat, constant volume 68 // cp , Specific heat, constant pressure 69 // vortex_strength , Strength of vortex 70 // center , Location of bubble center 71 // gamma = cp / cv, Specific heat ratio 72 // 73 // ***************************************************************************** 74 75 // ***************************************************************************** 76 // This helper function provides support for the exact, time-dependent solution 77 // (currently not implemented) and IC formulation for Euler traveling vortex 78 // ***************************************************************************** 79 CEED_QFUNCTION_HELPER int Exact_Euler(CeedInt dim, CeedScalar time, 80 const CeedScalar X[], CeedInt Nf, CeedScalar q[], 81 void *ctx) { 82 // Context 83 const EulerContext context = (EulerContext)ctx; 84 const CeedScalar vortex_strength = context->vortex_strength; 85 const CeedScalar *center = context->center; // Center of the domain 86 const CeedScalar *mean_velocity = context->mean_velocity; 87 88 // Setup 89 const CeedScalar gamma = 1.4; 90 const CeedScalar cv = 2.5; 91 const CeedScalar R = 1.; 92 const CeedScalar x = X[0], y = X[1]; // Coordinates 93 // Vortex center 94 const CeedScalar xc = center[0] + mean_velocity[0] * time; 95 const CeedScalar yc = center[1] + mean_velocity[1] * time; 96 97 const CeedScalar x0 = x - xc; 98 const CeedScalar y0 = y - yc; 99 const CeedScalar r = sqrt( x0*x0 + y0*y0 ); 100 const CeedScalar C = vortex_strength * exp((1. - r*r)/2.) / (2. * M_PI); 101 const CeedScalar delta_T = - (gamma - 1.) * vortex_strength * vortex_strength * 102 exp(1 - r*r) / (8. * gamma * M_PI * M_PI); 103 const CeedScalar S_vortex = 1; // no perturbation in the entropy P / rho^gamma 104 const CeedScalar S_bubble = (gamma - 1.) * vortex_strength * vortex_strength / 105 (8.*gamma*M_PI*M_PI); 106 CeedScalar rho, P, T, E, u[3] = {0.}; 107 108 // Initial Conditions 109 switch (context->euler_test) { 110 case 0: // Traveling vortex 111 T = 1 + delta_T; 112 // P = rho * T 113 // P = S * rho^gamma 114 // Solve for rho, then substitute for P 115 rho = pow(T/S_vortex, 1 / (gamma - 1.)); 116 P = rho * T; 117 u[0] = mean_velocity[0] - C*y0; 118 u[1] = mean_velocity[1] + C*x0; 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] = P / (gamma - 1.) + rho * (u[0]*u[0] + u[1]*u[1]) / 2.; 126 break; 127 case 1: // Constant zero velocity, density constant, total energy constant 128 rho = 1.; 129 E = 2.; 130 131 // Assign exact solution 132 q[0] = rho; 133 q[1] = rho * u[0]; 134 q[2] = rho * u[1]; 135 q[3] = rho * u[2]; 136 q[4] = E; 137 break; 138 case 2: // Constant nonzero velocity, density constant, total energy constant 139 rho = 1.; 140 E = 2.; 141 u[0] = mean_velocity[0]; 142 u[1] = mean_velocity[1]; 143 144 // Assign exact solution 145 q[0] = rho; 146 q[1] = rho * u[0]; 147 q[2] = rho * u[1]; 148 q[3] = rho * u[2]; 149 q[4] = E; 150 break; 151 case 3: // Velocity zero, pressure constant 152 // (so density and internal energy will be non-constant), 153 // but the velocity should stay zero and the bubble won't diffuse 154 // (for Euler, where there is no thermal conductivity) 155 P = 1.; 156 T = 1. - S_bubble * exp(1. - r*r); 157 rho = P / (R*T); 158 159 // Assign exact solution 160 q[0] = rho; 161 q[1] = rho * u[0]; 162 q[2] = rho * u[1]; 163 q[3] = rho * u[2]; 164 q[4] = rho * (cv * T + (u[0]*u[0] + u[1]*u[1])/2.); 165 break; 166 case 4: // Constant nonzero velocity, pressure constant 167 // (so density and internal energy will be non-constant), 168 // it should be transported across the domain, but velocity stays constant 169 P = 1.; 170 T = 1. - S_bubble * exp(1. - r*r); 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 (int 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 (int j=0; j<3; j++) { 363 for (int k=0; k<3; k++) { 364 drhodx[j] += drho[k] * dXdx[k][j]; 365 dEdx[j] += dE[k] * dXdx[k][j]; 366 for (int 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 (int j=0; j<5; j++) { 381 v[j][i] = 0.; 382 for (int k=0; k<3; k++) 383 dv[k][j][i] = 0.; 384 } 385 386 // -- Density 387 // ---- u rho 388 for (int 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 (int j=0; j<3; j++) 394 for (int 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 (int 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 // ---- Transpose of the Jacobian 410 CeedScalar jacob_F_conv_T[3][5][5]; 411 for (int j=0; j<3; j++) 412 for (int k=0; k<5; k++) 413 for (int l=0; l<5; l++) 414 jacob_F_conv_T[j][k][l] = jacob_F_conv[j][l][k]; 415 416 // ---- dqdx collects drhodx, dUdx and dEdx in one vector 417 CeedScalar dqdx[5][3]; 418 for (int j=0; j<3; j++) { 419 dqdx[0][j] = drhodx[j]; 420 dqdx[4][j] = dEdx[j]; 421 for (int k=0; k<3; k++) 422 dqdx[k+1][j] = dUdx[k][j]; 423 } 424 425 // ---- strong_conv = dF/dq * dq/dx (Strong convection) 426 CeedScalar strong_conv[5] = {0.}; 427 for (int j=0; j<3; j++) 428 for (int k=0; k<5; k++) 429 for (int l=0; l<5; l++) 430 strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; 431 432 // Stabilization 433 // -- Tau elements 434 const CeedScalar sound_speed = sqrt(gamma * P / rho); 435 CeedScalar Tau_x[3] = {0.}; 436 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau); 437 438 // -- Stabilization method: none or SU 439 CeedScalar stab[5][3]; 440 switch (context->stabilization) { 441 case 0: // Galerkin 442 break; 443 case 1: // SU 444 for (int j=0; j<3; j++) 445 for (int k=0; k<5; k++) 446 for (int l=0; l<5; l++) 447 stab[k][j] = jacob_F_conv_T[j][k][l] * Tau_x[j] * strong_conv[l]; 448 449 for (int j=0; j<5; j++) 450 for (int k=0; k<3; k++) 451 dv[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 452 stab[j][1] * dXdx[k][1] + 453 stab[j][2] * dXdx[k][2]); 454 break; 455 case 2: // SUPG is not implemented for explicit scheme 456 break; 457 } 458 459 } // End Quadrature Point Loop 460 461 // Return 462 return 0; 463 } 464 465 // ***************************************************************************** 466 // This QFunction implements the Euler equations with (mentioned above) 467 // with implicit time stepping method 468 // 469 // ***************************************************************************** 470 CEED_QFUNCTION(IFunction_Euler)(void *ctx, CeedInt Q, 471 const CeedScalar *const *in, CeedScalar *const *out) { 472 // *INDENT-OFF* 473 // Inputs 474 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 475 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 476 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 477 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 478 // Outputs 479 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 480 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 481 482 EulerContext context = (EulerContext)ctx; 483 const CeedScalar c_tau = context->c_tau; 484 const CeedScalar gamma = 1.4; 485 486 CeedPragmaSIMD 487 // Quadrature Point Loop 488 for (CeedInt i=0; i<Q; i++) { 489 // *INDENT-OFF* 490 // Setup 491 // -- Interp in 492 const CeedScalar rho = q[0][i]; 493 const CeedScalar u[3] = {q[1][i] / rho, 494 q[2][i] / rho, 495 q[3][i] / rho 496 }; 497 const CeedScalar E = q[4][i]; 498 const CeedScalar drho[3] = {dq[0][0][i], 499 dq[1][0][i], 500 dq[2][0][i] 501 }; 502 const CeedScalar dU[3][3] = {{dq[0][1][i], 503 dq[1][1][i], 504 dq[2][1][i]}, 505 {dq[0][2][i], 506 dq[1][2][i], 507 dq[2][2][i]}, 508 {dq[0][3][i], 509 dq[1][3][i], 510 dq[2][3][i]} 511 }; 512 const CeedScalar dE[3] = {dq[0][4][i], 513 dq[1][4][i], 514 dq[2][4][i] 515 }; 516 // -- Interp-to-Interp q_data 517 const CeedScalar wdetJ = q_data[0][i]; 518 // -- Interp-to-Grad q_data 519 // ---- Inverse of change of coordinate matrix: X_i,j 520 // *INDENT-OFF* 521 const CeedScalar dXdx[3][3] = {{q_data[1][i], 522 q_data[2][i], 523 q_data[3][i]}, 524 {q_data[4][i], 525 q_data[5][i], 526 q_data[6][i]}, 527 {q_data[7][i], 528 q_data[8][i], 529 q_data[9][i]} 530 }; 531 // *INDENT-ON* 532 // dU/dx 533 CeedScalar drhodx[3] = {0.}; 534 CeedScalar dEdx[3] = {0.}; 535 CeedScalar dUdx[3][3] = {{0.}}; 536 CeedScalar dXdxdXdxT[3][3] = {{0.}}; 537 for (int j=0; j<3; j++) { 538 for (int k=0; k<3; k++) { 539 drhodx[j] += drho[k] * dXdx[k][j]; 540 dEdx[j] += dE[k] * dXdx[k][j]; 541 for (int l=0; l<3; l++) { 542 dUdx[j][k] += dU[j][l] * dXdx[l][k]; 543 dXdxdXdxT[j][k] += dXdx[j][l]*dXdx[k][l]; //dXdx_j,k * dXdx_k,j 544 } 545 } 546 } 547 const CeedScalar 548 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 549 E_internal = E - E_kinetic, 550 P = E_internal * (gamma - 1.); // P = pressure 551 552 // The Physics 553 // Zero v and dv so all future terms can safely sum into it 554 for (int j=0; j<5; j++) { 555 v[j][i] = 0.; 556 for (int k=0; k<3; k++) 557 dv[k][j][i] = 0.; 558 } 559 //-----mass matrix 560 for (int j=0; j<5; j++) 561 v[j][i] += wdetJ*q_dot[j][i]; 562 563 // -- Density 564 // ---- u rho 565 for (int j=0; j<3; j++) 566 dv[j][0][i] -= wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 567 rho*u[2]*dXdx[j][2]); 568 // -- Momentum 569 // ---- rho (u x u) + P I3 570 for (int j=0; j<3; j++) 571 for (int k=0; k<3; k++) 572 dv[k][j+1][i] -= wdetJ*((rho*u[j]*u[0] + (j==0?P:0.))*dXdx[k][0] + 573 (rho*u[j]*u[1] + (j==1?P:0.))*dXdx[k][1] + 574 (rho*u[j]*u[2] + (j==2?P:0.))*dXdx[k][2]); 575 // -- Total Energy Density 576 // ---- (E + P) u 577 for (int j=0; j<3; j++) 578 dv[j][4][i] -= wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 579 u[2]*dXdx[j][2]); 580 581 // -- Stabilization terms 582 // ---- jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 583 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}}; 584 ConvectiveFluxJacobian_Euler(jacob_F_conv, rho, u, E, gamma); 585 586 // ---- Transpose of the Jacobian 587 CeedScalar jacob_F_conv_T[3][5][5]; 588 for (int j=0; j<3; j++) 589 for (int k=0; k<5; k++) 590 for (int l=0; l<5; l++) 591 jacob_F_conv_T[j][k][l] = jacob_F_conv[j][l][k]; 592 593 // ---- dqdx collects drhodx, dUdx and dEdx in one vector 594 CeedScalar dqdx[5][3]; 595 for (int j=0; j<3; j++) { 596 dqdx[0][j] = drhodx[j]; 597 dqdx[4][j] = dEdx[j]; 598 for (int k=0; k<3; k++) 599 dqdx[k+1][j] = dUdx[k][j]; 600 } 601 602 // ---- strong_conv = dF/dq * dq/dx (Strong convection) 603 CeedScalar strong_conv[5] = {0.}; 604 for (int j=0; j<3; j++) 605 for (int k=0; k<5; k++) 606 for (int l=0; l<5; l++) 607 strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; 608 609 // ---- Strong residual 610 CeedScalar strong_res[5]; 611 for (int j=0; j<5; j++) 612 strong_res[j] = q_dot[j][i] + strong_conv[j]; 613 614 // Stabilization 615 // -- Tau elements 616 const CeedScalar sound_speed = sqrt(gamma * P / rho); 617 CeedScalar Tau_x[3] = {0.}; 618 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau); 619 620 // -- Stabilization method: none, SU, or SUPG 621 CeedScalar stab[5][3]; 622 switch (context->stabilization) { 623 case 0: // Galerkin 624 break; 625 case 1: // SU 626 for (int j=0; j<3; j++) 627 for (int k=0; k<5; k++) 628 for (int l=0; l<5; l++) 629 stab[k][j] = jacob_F_conv_T[j][k][l] * Tau_x[j] * strong_conv[l]; 630 631 for (int j=0; j<5; j++) 632 for (int k=0; k<3; k++) 633 dv[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 634 stab[j][1] * dXdx[k][1] + 635 stab[j][2] * dXdx[k][2]); 636 break; 637 case 2: // SUPG 638 for (int j=0; j<3; j++) 639 for (int k=0; k<5; k++) 640 for (int l=0; l<5; l++) 641 stab[k][j] = jacob_F_conv_T[j][k][l] * Tau_x[j] * strong_res[l]; 642 643 for (int j=0; j<5; j++) 644 for (int k=0; k<3; k++) 645 dv[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 646 stab[j][1] * dXdx[k][1] + 647 stab[j][2] * dXdx[k][2]); 648 break; 649 } 650 } // End Quadrature Point Loop 651 652 // Return 653 return 0; 654 } 655 // ***************************************************************************** 656 // This QFunction sets the boundary conditions 657 // In this problem, only in/outflow BCs are implemented 658 // 659 // Inflow and outflow faces are determined based on 660 // sign(dot(mean_velocity, normal)): 661 // sign(dot(mean_velocity, normal)) > 0 : outflow BCs 662 // sign(dot(mean_velocity, normal)) < 0 : inflow BCs 663 // 664 // Outflow BCs: 665 // The validity of the weak form of the governing equations is 666 // extended to the outflow. 667 // 668 // Inflow BCs: 669 // Prescribed T_inlet and P_inlet are converted to conservative variables 670 // and applied weakly. 671 // 672 // ***************************************************************************** 673 CEED_QFUNCTION(Euler_Sur)(void *ctx, CeedInt Q, 674 const CeedScalar *const *in, 675 CeedScalar *const *out) { 676 // *INDENT-OFF* 677 // Inputs 678 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 679 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 680 // Outputs 681 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 682 // *INDENT-ON* 683 EulerContext context = (EulerContext)ctx; 684 const int euler_test = context->euler_test; 685 const bool implicit = context->implicit; 686 CeedScalar *mean_velocity = context->mean_velocity; 687 688 const CeedScalar gamma = 1.4; 689 const CeedScalar cv = 2.5; 690 const CeedScalar R = 1.; 691 CeedScalar T_inlet; 692 CeedScalar P_inlet; 693 694 // For test cases 1 and 3 the background velocity is zero 695 if (euler_test == 1 || euler_test == 3) 696 for (CeedInt i=0; i<3; i++) mean_velocity[i] = 0.; 697 698 // For test cases 1 and 2, T_inlet = T_inlet = 0.4 699 if (euler_test == 1 || euler_test == 2) T_inlet = P_inlet = .4; 700 else T_inlet = P_inlet = 1.; 701 702 CeedPragmaSIMD 703 // Quadrature Point Loop 704 for (CeedInt i=0; i<Q; i++) { 705 // Setup 706 // -- Interp in 707 const CeedScalar rho = q[0][i]; 708 const CeedScalar u[3] = {q[1][i] / rho, 709 q[2][i] / rho, 710 q[3][i] / rho 711 }; 712 const CeedScalar E = q[4][i]; 713 714 // -- Interp-to-Interp q_data 715 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 716 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 717 // We can effect this by swapping the sign on this weight 718 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 719 // ---- Normal vectors 720 const CeedScalar norm[3] = {q_data_sur[1][i], 721 q_data_sur[2][i], 722 q_data_sur[3][i] 723 }; 724 725 // face_normal = Normal vector of the face 726 const CeedScalar face_normal = norm[0]*mean_velocity[0] + 727 norm[1]*mean_velocity[1] + 728 norm[2]*mean_velocity[2]; 729 // The Physics 730 // Zero v so all future terms can safely sum into it 731 for (int j=0; j<5; j++) v[j][i] = 0.; 732 733 // Implementing in/outflow BCs 734 if (face_normal > 0) { // outflow 735 const CeedScalar E_kinetic = (u[0]*u[0] + u[1]*u[1]) / 2.; 736 const CeedScalar P = (E - E_kinetic * rho) * (gamma - 1.); // pressure 737 const CeedScalar u_normal = norm[0]*u[0] + norm[1]*u[1] + 738 norm[2]*u[2]; // Normal velocity 739 // The Physics 740 // -- Density 741 v[0][i] -= wdetJb * rho * u_normal; 742 743 // -- Momentum 744 for (int j=0; j<3; j++) 745 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + norm[j] * P); 746 747 // -- Total Energy Density 748 v[4][i] -= wdetJb * u_normal * (E + P); 749 750 } else { // inflow 751 const CeedScalar rho_inlet = P_inlet/(R*T_inlet); 752 const CeedScalar E_kinetic_inlet = (mean_velocity[0]*mean_velocity[0] + 753 mean_velocity[1]*mean_velocity[1]) / 2.; 754 // incoming total energy 755 const CeedScalar E_inlet = rho_inlet * (cv * T_inlet + E_kinetic_inlet); 756 757 // The Physics 758 // -- Density 759 v[0][i] -= wdetJb * rho_inlet * face_normal; 760 761 // -- Momentum 762 for (int j=0; j<3; j++) 763 v[j+1][i] -= wdetJb *(rho_inlet * face_normal * mean_velocity[j] + 764 norm[j] * P_inlet); 765 766 // -- Total Energy Density 767 v[4][i] -= wdetJb * face_normal * (E_inlet + P_inlet); 768 } 769 770 } // End Quadrature Point Loop 771 return 0; 772 } 773 774 // ***************************************************************************** 775 776 #endif // eulervortex_h 777