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 int euler_test; 43 bool implicit; 44 }; 45 #endif 46 47 // ***************************************************************************** 48 // This function sets the initial conditions 49 // 50 // Temperature: 51 // T = 1 - (gamma - 1) vortex_strength**2 exp(1 - r**2) / (8 gamma pi**2) 52 // Density: 53 // rho = (T/S_vortex)^(1 / (gamma - 1)) 54 // Pressure: 55 // P = rho * T 56 // Velocity: 57 // ui = 1 + vortex_strength exp((1 - r**2)/2.) [yc - y, x - xc] / (2 pi) 58 // r = sqrt( (x - xc)**2 + (y - yc)**2 ) 59 // Velocity/Momentum Density: 60 // Ui = rho ui 61 // Total Energy: 62 // E = P / (gamma - 1) + rho (u u)/2 63 // 64 // Constants: 65 // cv , Specific heat, constant volume 66 // cp , Specific heat, constant pressure 67 // vortex_strength , Strength of vortex 68 // center , Location of bubble center 69 // gamma = cp / cv, Specific heat ratio 70 // 71 // ***************************************************************************** 72 73 // ***************************************************************************** 74 // This helper function provides support for the exact, time-dependent solution 75 // (currently not implemented) and IC formulation for Euler traveling vortex 76 // ***************************************************************************** 77 CEED_QFUNCTION_HELPER int Exact_Euler(CeedInt dim, CeedScalar time, 78 const CeedScalar X[], CeedInt Nf, CeedScalar q[], 79 void *ctx) { 80 // Context 81 const EulerContext context = (EulerContext)ctx; 82 const CeedScalar vortex_strength = context->vortex_strength; 83 const CeedScalar *center = context->center; // Center of the domain 84 const CeedScalar *mean_velocity = context->mean_velocity; 85 86 // Setup 87 const CeedScalar gamma = 1.4; 88 const CeedScalar cv = 2.5; 89 const CeedScalar R = 1.; 90 const CeedScalar x = X[0], y = X[1]; // Coordinates 91 // Vortex center 92 const CeedScalar xc = center[0] + mean_velocity[0] * time; 93 const CeedScalar yc = center[1] + mean_velocity[1] * time; 94 95 const CeedScalar x0 = x - xc; 96 const CeedScalar y0 = y - yc; 97 const CeedScalar r = sqrt( x0*x0 + y0*y0 ); 98 const CeedScalar C = vortex_strength * exp((1. - r*r)/2.) / (2. * M_PI); 99 const CeedScalar delta_T = - (gamma - 1) * vortex_strength * vortex_strength * 100 exp(1 - r*r) / (8 * gamma * M_PI * M_PI); 101 const CeedScalar S_vortex = 1; // no perturbation in the entropy P / rho^gamma 102 const CeedScalar S_bubble = (gamma - 1.) * vortex_strength * vortex_strength / 103 (8.*gamma*M_PI*M_PI); 104 CeedScalar rho, P, T, E, u[3] = {0.}; 105 106 // Initial Conditions 107 switch (context->euler_test) { 108 case 0: // Traveling vortex 109 T = 1 + delta_T; 110 // P = rho * T 111 // P = S * rho^gamma 112 // Solve for rho, then substitute for P 113 rho = pow(T/S_vortex, 1 / (gamma - 1)); 114 P = rho * T; 115 u[0] = mean_velocity[0] - C*y0; 116 u[1] = mean_velocity[1] + C*x0; 117 118 // Assign exact solution 119 q[0] = rho; 120 q[1] = rho * u[0]; 121 q[2] = rho * u[1]; 122 q[3] = rho * u[2]; 123 q[4] = P / (gamma - 1.) + rho * (u[0]*u[0] + u[1]*u[1]) / 2.; 124 break; 125 case 1: // Constant zero velocity, density constant, total energy constant 126 rho = 1.; 127 E = 2.; 128 129 // Assign exact solution 130 q[0] = rho; 131 q[1] = rho * u[0]; 132 q[2] = rho * u[1]; 133 q[3] = rho * u[2]; 134 q[4] = E; 135 break; 136 case 2: // Constant nonzero velocity, density constant, total energy constant 137 rho = 1.; 138 E = 2.; 139 u[0] = mean_velocity[0]; 140 u[1] = mean_velocity[1]; 141 142 // Assign exact solution 143 q[0] = rho; 144 q[1] = rho * u[0]; 145 q[2] = rho * u[1]; 146 q[3] = rho * u[2]; 147 q[4] = E; 148 break; 149 case 3: // Velocity zero, pressure constant 150 // (so density and internal energy will be non-constant), 151 // but the velocity should stay zero and the bubble won't diffuse 152 // (for Euler, where there is no thermal conductivity) 153 P = 1.; 154 T = 1. - S_bubble * exp(1. - r*r); 155 rho = P / (R*T); 156 157 // Assign exact solution 158 q[0] = rho; 159 q[1] = rho * u[0]; 160 q[2] = rho * u[1]; 161 q[3] = rho * u[2]; 162 q[4] = rho * (cv * T + (u[0]*u[0] + u[1]*u[1])/2.); 163 break; 164 case 4: // Constant nonzero velocity, pressure constant 165 // (so density and internal energy will be non-constant), 166 // it should be transported across the domain, but velocity stays constant 167 P = 1.; 168 T = 1. - S_bubble * exp(1. - r*r); 169 rho = P / (R*T); 170 u[0] = mean_velocity[0]; 171 u[1] = mean_velocity[1]; 172 173 // Assign exact solution 174 q[0] = rho; 175 q[1] = rho * u[0]; 176 q[2] = rho * u[1]; 177 q[3] = rho * u[2]; 178 q[4] = rho * (cv * T + (u[0]*u[0] + u[1]*u[1])/2.); 179 break; 180 } 181 // Return 182 return 0; 183 } 184 185 // ***************************************************************************** 186 // This QFunction sets the initial conditions for Euler traveling vortex 187 // ***************************************************************************** 188 CEED_QFUNCTION(ICsEuler)(void *ctx, CeedInt Q, 189 const CeedScalar *const *in, CeedScalar *const *out) { 190 // Inputs 191 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 192 193 // Outputs 194 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 195 const EulerContext context = (EulerContext)ctx; 196 197 CeedPragmaSIMD 198 // Quadrature Point Loop 199 for (CeedInt i=0; i<Q; i++) { 200 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 201 CeedScalar q[5]; 202 203 Exact_Euler(3, context->curr_time, x, 5, q, ctx); 204 205 for (CeedInt j=0; j<5; j++) 206 q0[j][i] = q[j]; 207 } // End of Quadrature Point Loop 208 209 // Return 210 return 0; 211 } 212 213 // ***************************************************************************** 214 // This QFunction implements the following formulation of Euler equations 215 // with explicit time stepping method 216 // 217 // This is 3D Euler for compressible gas dynamics in conservation 218 // form with state variables of density, momentum density, and total 219 // energy density. 220 // 221 // State Variables: q = ( rho, U1, U2, U3, E ) 222 // rho - Mass Density 223 // Ui - Momentum Density, Ui = rho ui 224 // E - Total Energy Density, E = P / (gamma - 1) + rho (u u)/2 225 // 226 // Euler Equations: 227 // drho/dt + div( U ) = 0 228 // dU/dt + div( rho (u x u) + P I3 ) = 0 229 // dE/dt + div( (E + P) u ) = 0 230 // 231 // Equation of State: 232 // P = (gamma - 1) (E - rho (u u) / 2) 233 // 234 // Constants: 235 // cv , Specific heat, constant volume 236 // cp , Specific heat, constant pressure 237 // g , Gravity 238 // gamma = cp / cv, Specific heat ratio 239 // ***************************************************************************** 240 CEED_QFUNCTION(Euler)(void *ctx, CeedInt Q, 241 const CeedScalar *const *in, CeedScalar *const *out) { 242 // *INDENT-OFF* 243 // Inputs 244 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 245 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 246 // Outputs 247 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 248 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 249 250 const CeedScalar gamma = 1.4; 251 252 CeedPragmaSIMD 253 // Quadrature Point Loop 254 for (CeedInt i=0; i<Q; i++) { 255 // *INDENT-OFF* 256 // Setup 257 // -- Interp in 258 const CeedScalar rho = q[0][i]; 259 const CeedScalar u[3] = {q[1][i] / rho, 260 q[2][i] / rho, 261 q[3][i] / rho 262 }; 263 const CeedScalar E = q[4][i]; 264 265 // -- Interp-to-Interp q_data 266 const CeedScalar wdetJ = q_data[0][i]; 267 // -- Interp-to-Grad q_data 268 // ---- Inverse of change of coordinate matrix: X_i,j 269 // *INDENT-OFF* 270 const CeedScalar dXdx[3][3] = {{q_data[1][i], 271 q_data[2][i], 272 q_data[3][i]}, 273 {q_data[4][i], 274 q_data[5][i], 275 q_data[6][i]}, 276 {q_data[7][i], 277 q_data[8][i], 278 q_data[9][i]} 279 }; 280 // *INDENT-ON* 281 const CeedScalar 282 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 283 E_internal = E - E_kinetic, 284 P = E_internal * (gamma - 1); // P = pressure 285 286 // The Physics 287 // Zero v and dv so all future terms can safely sum into it 288 for (int j=0; j<5; j++) { 289 v[j][i] = 0; 290 for (int k=0; k<3; k++) 291 dv[k][j][i] = 0; 292 } 293 294 // -- Density 295 // ---- u rho 296 for (int j=0; j<3; j++) 297 dv[j][0][i] += wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 298 rho*u[2]*dXdx[j][2]); 299 // -- Momentum 300 // ---- rho (u x u) + P I3 301 for (int j=0; j<3; j++) 302 for (int k=0; k<3; k++) 303 dv[k][j+1][i] += wdetJ*((rho*u[j]*u[0] + (j==0?P:0))*dXdx[k][0] + 304 (rho*u[j]*u[1] + (j==1?P:0))*dXdx[k][1] + 305 (rho*u[j]*u[2] + (j==2?P:0))*dXdx[k][2]); 306 // -- Total Energy Density 307 // ---- (E + P) u 308 for (int j=0; j<3; j++) 309 dv[j][4][i] += wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 310 u[2]*dXdx[j][2]); 311 } // End Quadrature Point Loop 312 313 // Return 314 return 0; 315 } 316 317 // ***************************************************************************** 318 // This QFunction implements the Euler equations with (mentioned above) 319 // with implicit time stepping method 320 // 321 // ***************************************************************************** 322 CEED_QFUNCTION(IFunction_Euler)(void *ctx, CeedInt Q, 323 const CeedScalar *const *in, CeedScalar *const *out) { 324 // *INDENT-OFF* 325 // Inputs 326 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 327 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 328 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 329 // Outputs 330 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 331 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 332 333 const CeedScalar gamma = 1.4; 334 335 CeedPragmaSIMD 336 // Quadrature Point Loop 337 for (CeedInt i=0; i<Q; i++) { 338 // *INDENT-OFF* 339 // Setup 340 // -- Interp in 341 const CeedScalar rho = q[0][i]; 342 const CeedScalar u[3] = {q[1][i] / rho, 343 q[2][i] / rho, 344 q[3][i] / rho 345 }; 346 const CeedScalar E = q[4][i]; 347 348 // -- Interp-to-Interp q_data 349 const CeedScalar wdetJ = q_data[0][i]; 350 // -- Interp-to-Grad q_data 351 // ---- Inverse of change of coordinate matrix: X_i,j 352 // *INDENT-OFF* 353 const CeedScalar dXdx[3][3] = {{q_data[1][i], 354 q_data[2][i], 355 q_data[3][i]}, 356 {q_data[4][i], 357 q_data[5][i], 358 q_data[6][i]}, 359 {q_data[7][i], 360 q_data[8][i], 361 q_data[9][i]} 362 }; 363 // *INDENT-ON* 364 const CeedScalar 365 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 366 E_internal = E - E_kinetic, 367 P = E_internal * (gamma - 1); // P = pressure 368 369 // The Physics 370 // Zero v and dv so all future terms can safely sum into it 371 for (int j=0; j<5; j++) { 372 v[j][i] = 0; 373 for (int k=0; k<3; k++) 374 dv[k][j][i] = 0; 375 } 376 //-----mass matrix 377 for (int j=0; j<5; j++) 378 v[j][i] += wdetJ*q_dot[j][i]; 379 380 // -- Density 381 // ---- u rho 382 for (int j=0; j<3; j++) 383 dv[j][0][i] -= wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 384 rho*u[2]*dXdx[j][2]); 385 // -- Momentum 386 // ---- rho (u x u) + P I3 387 for (int j=0; j<3; j++) 388 for (int k=0; k<3; k++) 389 dv[k][j+1][i] -= wdetJ*((rho*u[j]*u[0] + (j==0?P:0))*dXdx[k][0] + 390 (rho*u[j]*u[1] + (j==1?P:0))*dXdx[k][1] + 391 (rho*u[j]*u[2] + (j==2?P:0))*dXdx[k][2]); 392 // -- Total Energy Density 393 // ---- (E + P) u 394 for (int j=0; j<3; j++) 395 dv[j][4][i] -= wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 396 u[2]*dXdx[j][2]); 397 } // End Quadrature Point Loop 398 399 // Return 400 return 0; 401 } 402 // ***************************************************************************** 403 // This QFunction sets the boundary conditions 404 // In this problem, only in/outflow BCs are implemented 405 // 406 // Inflow and outflow faces are determined based on 407 // sign(dot(mean_velocity, normal)): 408 // sign(dot(mean_velocity, normal)) > 0 : outflow BCs 409 // sign(dot(mean_velocity, normal)) < 0 : inflow BCs 410 // 411 // Outflow BCs: 412 // The validity of the weak form of the governing equations is 413 // extended to the outflow. 414 // 415 // Inflow BCs: 416 // Prescribed T_inlet and P_inlet are converted to conservative variables 417 // and applied weakly. 418 // 419 // ***************************************************************************** 420 CEED_QFUNCTION(Euler_Sur)(void *ctx, CeedInt Q, 421 const CeedScalar *const *in, 422 CeedScalar *const *out) { 423 // *INDENT-OFF* 424 // Inputs 425 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 426 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 427 // Outputs 428 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 429 // *INDENT-ON* 430 EulerContext context = (EulerContext)ctx; 431 const int euler_test = context->euler_test; 432 const bool implicit = context->implicit; 433 CeedScalar *mean_velocity = context->mean_velocity; 434 435 const CeedScalar gamma = 1.4; 436 const CeedScalar cv = 2.5; 437 const CeedScalar R = 1.; 438 CeedScalar T_inlet; 439 CeedScalar P_inlet; 440 441 // For test cases 1 and 3 the background velocity is zero 442 if (euler_test == 1 || euler_test == 3) 443 for (CeedInt i=0; i<3; i++) mean_velocity[i] = 0.; 444 445 // For test cases 1 and 2, T_inlet = T_inlet = 0.4 446 if (euler_test == 1 || euler_test == 2) T_inlet = P_inlet = .4; 447 else T_inlet = P_inlet = 1.; 448 449 CeedPragmaSIMD 450 // Quadrature Point Loop 451 for (CeedInt i=0; i<Q; i++) { 452 // Setup 453 // -- Interp in 454 const CeedScalar rho = q[0][i]; 455 const CeedScalar u[3] = {q[1][i] / rho, 456 q[2][i] / rho, 457 q[3][i] / rho 458 }; 459 const CeedScalar E = q[4][i]; 460 461 // -- Interp-to-Interp q_data 462 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 463 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 464 // We can effect this by swapping the sign on this weight 465 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 466 // ---- Normal vectors 467 const CeedScalar norm[3] = {q_data_sur[1][i], 468 q_data_sur[2][i], 469 q_data_sur[3][i] 470 }; 471 472 // face_normal = Normal vector of the face 473 const CeedScalar face_normal = norm[0]*mean_velocity[0] + 474 norm[1]*mean_velocity[1] + 475 norm[2]*mean_velocity[2]; 476 // The Physics 477 // Zero v so all future terms can safely sum into it 478 for (int j=0; j<5; j++) v[j][i] = 0; 479 480 // Implementing in/outflow BCs 481 if (face_normal > 0) { // outflow 482 const CeedScalar E_kinetic = (u[0]*u[0] + u[1]*u[1]) / 2.; 483 const CeedScalar P = (E - E_kinetic * rho) * (gamma - 1.); // pressure 484 const CeedScalar u_normal = norm[0]*u[0] + norm[1]*u[1] + 485 norm[2]*u[2]; // Normal velocity 486 // The Physics 487 // -- Density 488 v[0][i] -= wdetJb * rho * u_normal; 489 490 // -- Momentum 491 for (int j=0; j<3; j++) 492 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + norm[j] * P); 493 494 // -- Total Energy Density 495 v[4][i] -= wdetJb * u_normal * (E + P); 496 497 } else { // inflow 498 const CeedScalar rho_inlet = P_inlet/(R*T_inlet); 499 const CeedScalar E_kinetic_inlet = (mean_velocity[0]*mean_velocity[0] + 500 mean_velocity[1]*mean_velocity[1]) / 2.; 501 // incoming total energy 502 const CeedScalar E_inlet = rho_inlet * (cv * T_inlet + E_kinetic_inlet); 503 504 // The Physics 505 // -- Density 506 v[0][i] -= wdetJb * rho_inlet * face_normal; 507 508 // -- Momentum 509 for (int j=0; j<3; j++) 510 v[j+1][i] -= wdetJb *(rho_inlet * face_normal * mean_velocity[j] + 511 norm[j] * P_inlet); 512 513 // -- Total Energy Density 514 v[4][i] -= wdetJb * face_normal * (E_inlet + P_inlet); 515 } 516 517 } // End Quadrature Point Loop 518 return 0; 519 } 520 521 // ***************************************************************************** 522 523 #endif // eulervortex_h 524