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