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 #ifndef __CUDACC__ 29 # include <math.h> 30 #endif 31 32 #ifndef M_PI 33 #define M_PI 3.14159265358979323846 34 #endif 35 36 #ifndef euler_context_struct 37 #define euler_context_struct 38 typedef struct EulerContext_ *EulerContext; 39 struct EulerContext_ { 40 CeedScalar center[3]; 41 CeedScalar curr_time; 42 CeedScalar vortex_strength; 43 CeedScalar mean_velocity[3]; 44 int euler_test; 45 bool implicit; 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 } 183 // Return 184 return 0; 185 } 186 187 // ***************************************************************************** 188 // This QFunction sets the initial conditions for Euler traveling vortex 189 // ***************************************************************************** 190 CEED_QFUNCTION(ICsEuler)(void *ctx, CeedInt Q, 191 const CeedScalar *const *in, CeedScalar *const *out) { 192 // Inputs 193 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 194 195 // Outputs 196 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 197 const EulerContext context = (EulerContext)ctx; 198 199 CeedPragmaSIMD 200 // Quadrature Point Loop 201 for (CeedInt i=0; i<Q; i++) { 202 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 203 CeedScalar q[5]; 204 205 Exact_Euler(3, context->curr_time, x, 5, q, ctx); 206 207 for (CeedInt j=0; j<5; j++) 208 q0[j][i] = q[j]; 209 } // End of Quadrature Point Loop 210 211 // Return 212 return 0; 213 } 214 215 // ***************************************************************************** 216 // This QFunction implements the following formulation of Euler equations 217 // with explicit time stepping method 218 // 219 // This is 3D Euler for compressible gas dynamics in conservation 220 // form with state variables of density, momentum density, and total 221 // energy density. 222 // 223 // State Variables: q = ( rho, U1, U2, U3, E ) 224 // rho - Mass Density 225 // Ui - Momentum Density, Ui = rho ui 226 // E - Total Energy Density, E = P / (gamma - 1) + rho (u u)/2 227 // 228 // Euler Equations: 229 // drho/dt + div( U ) = 0 230 // dU/dt + div( rho (u x u) + P I3 ) = 0 231 // dE/dt + div( (E + P) u ) = 0 232 // 233 // Equation of State: 234 // P = (gamma - 1) (E - rho (u u) / 2) 235 // 236 // Constants: 237 // cv , Specific heat, constant volume 238 // cp , Specific heat, constant pressure 239 // g , Gravity 240 // gamma = cp / cv, Specific heat ratio 241 // ***************************************************************************** 242 CEED_QFUNCTION(Euler)(void *ctx, CeedInt Q, 243 const CeedScalar *const *in, CeedScalar *const *out) { 244 // *INDENT-OFF* 245 // Inputs 246 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 247 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 248 // Outputs 249 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 250 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 251 252 const CeedScalar gamma = 1.4; 253 254 CeedPragmaSIMD 255 // Quadrature Point Loop 256 for (CeedInt i=0; i<Q; i++) { 257 // *INDENT-OFF* 258 // Setup 259 // -- Interp in 260 const CeedScalar rho = q[0][i]; 261 const CeedScalar u[3] = {q[1][i] / rho, 262 q[2][i] / rho, 263 q[3][i] / rho 264 }; 265 const CeedScalar E = q[4][i]; 266 267 // -- Interp-to-Interp q_data 268 const CeedScalar wdetJ = q_data[0][i]; 269 // -- Interp-to-Grad q_data 270 // ---- Inverse of change of coordinate matrix: X_i,j 271 // *INDENT-OFF* 272 const CeedScalar dXdx[3][3] = {{q_data[1][i], 273 q_data[2][i], 274 q_data[3][i]}, 275 {q_data[4][i], 276 q_data[5][i], 277 q_data[6][i]}, 278 {q_data[7][i], 279 q_data[8][i], 280 q_data[9][i]} 281 }; 282 // *INDENT-ON* 283 const CeedScalar 284 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 285 E_internal = E - E_kinetic, 286 P = E_internal * (gamma - 1); // P = pressure 287 288 // The Physics 289 // Zero v and dv so all future terms can safely sum into it 290 for (int j=0; j<5; j++) { 291 v[j][i] = 0; 292 for (int k=0; k<3; k++) 293 dv[k][j][i] = 0; 294 } 295 296 // -- Density 297 // ---- u rho 298 for (int j=0; j<3; j++) 299 dv[j][0][i] += wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 300 rho*u[2]*dXdx[j][2]); 301 // -- Momentum 302 // ---- rho (u x u) + P I3 303 for (int j=0; j<3; j++) 304 for (int k=0; k<3; k++) 305 dv[k][j+1][i] += wdetJ*((rho*u[j]*u[0] + (j==0?P:0))*dXdx[k][0] + 306 (rho*u[j]*u[1] + (j==1?P:0))*dXdx[k][1] + 307 (rho*u[j]*u[2] + (j==2?P:0))*dXdx[k][2]); 308 // -- Total Energy Density 309 // ---- (E + P) u 310 for (int j=0; j<3; j++) 311 dv[j][4][i] += wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 312 u[2]*dXdx[j][2]); 313 } // End Quadrature Point Loop 314 315 // Return 316 return 0; 317 } 318 319 // ***************************************************************************** 320 // This QFunction implements the Euler equations with (mentioned above) 321 // with implicit time stepping method 322 // 323 // ***************************************************************************** 324 CEED_QFUNCTION(IFunction_Euler)(void *ctx, CeedInt Q, 325 const CeedScalar *const *in, CeedScalar *const *out) { 326 // *INDENT-OFF* 327 // Inputs 328 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 329 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 330 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 331 // Outputs 332 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 333 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 334 335 const CeedScalar gamma = 1.4; 336 337 CeedPragmaSIMD 338 // Quadrature Point Loop 339 for (CeedInt i=0; i<Q; i++) { 340 // *INDENT-OFF* 341 // Setup 342 // -- Interp in 343 const CeedScalar rho = q[0][i]; 344 const CeedScalar u[3] = {q[1][i] / rho, 345 q[2][i] / rho, 346 q[3][i] / rho 347 }; 348 const CeedScalar E = q[4][i]; 349 350 // -- Interp-to-Interp q_data 351 const CeedScalar wdetJ = q_data[0][i]; 352 // -- Interp-to-Grad q_data 353 // ---- Inverse of change of coordinate matrix: X_i,j 354 // *INDENT-OFF* 355 const CeedScalar dXdx[3][3] = {{q_data[1][i], 356 q_data[2][i], 357 q_data[3][i]}, 358 {q_data[4][i], 359 q_data[5][i], 360 q_data[6][i]}, 361 {q_data[7][i], 362 q_data[8][i], 363 q_data[9][i]} 364 }; 365 // *INDENT-ON* 366 const CeedScalar 367 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 368 E_internal = E - E_kinetic, 369 P = E_internal * (gamma - 1); // P = pressure 370 371 // The Physics 372 // Zero v and dv so all future terms can safely sum into it 373 for (int j=0; j<5; j++) { 374 v[j][i] = 0; 375 for (int k=0; k<3; k++) 376 dv[k][j][i] = 0; 377 } 378 //-----mass matrix 379 for (int j=0; j<5; j++) 380 v[j][i] += wdetJ*q_dot[j][i]; 381 382 // -- Density 383 // ---- u rho 384 for (int j=0; j<3; j++) 385 dv[j][0][i] -= wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 386 rho*u[2]*dXdx[j][2]); 387 // -- Momentum 388 // ---- rho (u x u) + P I3 389 for (int j=0; j<3; j++) 390 for (int k=0; k<3; k++) 391 dv[k][j+1][i] -= wdetJ*((rho*u[j]*u[0] + (j==0?P:0))*dXdx[k][0] + 392 (rho*u[j]*u[1] + (j==1?P:0))*dXdx[k][1] + 393 (rho*u[j]*u[2] + (j==2?P:0))*dXdx[k][2]); 394 // -- Total Energy Density 395 // ---- (E + P) u 396 for (int j=0; j<3; j++) 397 dv[j][4][i] -= wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 398 u[2]*dXdx[j][2]); 399 } // End Quadrature Point Loop 400 401 // Return 402 return 0; 403 } 404 // ***************************************************************************** 405 // This QFunction sets the boundary conditions 406 // In this problem, only in/outflow BCs are implemented 407 // 408 // Inflow and outflow faces are determined based on 409 // sign(dot(mean_velocity, normal)): 410 // sign(dot(mean_velocity, normal)) > 0 : outflow BCs 411 // sign(dot(mean_velocity, normal)) < 0 : inflow BCs 412 // 413 // Outflow BCs: 414 // The validity of the weak form of the governing equations is 415 // extended to the outflow. 416 // 417 // Inflow BCs: 418 // Prescribed T_inlet and P_inlet are converted to conservative variables 419 // and applied weakly. 420 // 421 // ***************************************************************************** 422 CEED_QFUNCTION(Euler_Sur)(void *ctx, CeedInt Q, 423 const CeedScalar *const *in, 424 CeedScalar *const *out) { 425 // *INDENT-OFF* 426 // Inputs 427 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 428 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 429 // Outputs 430 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 431 // *INDENT-ON* 432 EulerContext context = (EulerContext)ctx; 433 const int euler_test = context->euler_test; 434 const bool implicit = context->implicit; 435 CeedScalar *mean_velocity = context->mean_velocity; 436 437 const CeedScalar gamma = 1.4; 438 const CeedScalar cv = 2.5; 439 const CeedScalar R = 1.; 440 CeedScalar T_inlet; 441 CeedScalar P_inlet; 442 443 // For test cases 1 and 3 the background velocity is zero 444 if (euler_test == 1 || euler_test == 3) 445 for (CeedInt i=0; i<3; i++) mean_velocity[i] = 0.; 446 447 // For test cases 1 and 2, T_inlet = T_inlet = 0.4 448 if (euler_test == 1 || euler_test == 2) T_inlet = P_inlet = .4; 449 else T_inlet = P_inlet = 1.; 450 451 CeedPragmaSIMD 452 // Quadrature Point Loop 453 for (CeedInt i=0; i<Q; i++) { 454 // Setup 455 // -- Interp in 456 const CeedScalar rho = q[0][i]; 457 const CeedScalar u[3] = {q[1][i] / rho, 458 q[2][i] / rho, 459 q[3][i] / rho 460 }; 461 const CeedScalar E = q[4][i]; 462 463 // -- Interp-to-Interp q_data 464 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 465 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 466 // We can effect this by swapping the sign on this weight 467 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 468 // ---- Normal vectors 469 const CeedScalar norm[3] = {q_data_sur[1][i], 470 q_data_sur[2][i], 471 q_data_sur[3][i] 472 }; 473 474 // face_normal = Normal vector of the face 475 const CeedScalar face_normal = norm[0]*mean_velocity[0] + 476 norm[1]*mean_velocity[1] + 477 norm[2]*mean_velocity[2]; 478 // The Physics 479 // Zero v so all future terms can safely sum into it 480 for (int j=0; j<5; j++) v[j][i] = 0; 481 482 // Implementing in/outflow BCs 483 if (face_normal > 0) { // outflow 484 const CeedScalar E_kinetic = (u[0]*u[0] + u[1]*u[1]) / 2.; 485 const CeedScalar P = (E - E_kinetic * rho) * (gamma - 1.); // pressure 486 const CeedScalar u_normal = norm[0]*u[0] + norm[1]*u[1] + 487 norm[2]*u[2]; // Normal velocity 488 // The Physics 489 // -- Density 490 v[0][i] -= wdetJb * rho * u_normal; 491 492 // -- Momentum 493 for (int j=0; j<3; j++) 494 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + norm[j] * P); 495 496 // -- Total Energy Density 497 v[4][i] -= wdetJb * u_normal * (E + P); 498 499 } else { // inflow 500 const CeedScalar rho_inlet = P_inlet/(R*T_inlet); 501 const CeedScalar E_kinetic_inlet = (mean_velocity[0]*mean_velocity[0] + 502 mean_velocity[1]*mean_velocity[1]) / 2.; 503 // incoming total energy 504 const CeedScalar E_inlet = rho_inlet * (cv * T_inlet + E_kinetic_inlet); 505 506 // The Physics 507 // -- Density 508 v[0][i] -= wdetJb * rho_inlet * face_normal; 509 510 // -- Momentum 511 for (int j=0; j<3; j++) 512 v[j+1][i] -= wdetJb *(rho_inlet * face_normal * mean_velocity[j] + 513 norm[j] * P_inlet); 514 515 // -- Total Energy Density 516 v[4][i] -= wdetJb * face_normal * (E_inlet + P_inlet); 517 } 518 519 } // End Quadrature Point Loop 520 return 0; 521 } 522 523 // ***************************************************************************** 524 525 #endif // eulervortex_h 526