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