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