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