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 /// Advection initial condition and operator for Navier-Stokes example using PETSc 19 20 #ifndef advection_h 21 #define advection_h 22 23 #include <math.h> 24 25 #ifndef setup_context_struct 26 #define setup_context_struct 27 typedef struct SetupContext_ *SetupContext; 28 struct SetupContext_ { 29 CeedScalar theta0; 30 CeedScalar thetaC; 31 CeedScalar P0; 32 CeedScalar N; 33 CeedScalar cv; 34 CeedScalar cp; 35 CeedScalar Rd; 36 CeedScalar g; 37 CeedScalar rc; 38 CeedScalar lx; 39 CeedScalar ly; 40 CeedScalar lz; 41 CeedScalar center[3]; 42 CeedScalar dc_axis[3]; 43 CeedScalar wind[3]; 44 CeedScalar time; 45 int wind_type; // See WindType: 0=ROTATION, 1=TRANSLATION 46 int bubble_type; // See BubbleType: 0=SPHERE, 1=CYLINDER 47 int bubble_continuity_type; // See BubbleContinuityType: 0=SMOOTH, 1=BACK_SHARP 2=THICK 48 }; 49 #endif 50 51 #ifndef advection_context_struct 52 #define advection_context_struct 53 typedef struct AdvectionContext_ *AdvectionContext; 54 struct AdvectionContext_ { 55 CeedScalar CtauS; 56 CeedScalar strong_form; 57 CeedScalar E_wind; 58 bool implicit; 59 int stabilization; // See StabilizationType: 0=none, 1=SU, 2=SUPG 60 }; 61 #endif 62 63 // ***************************************************************************** 64 // This QFunction sets the initial conditions and the boundary conditions 65 // for two test cases: ROTATION and TRANSLATION 66 // 67 // -- ROTATION (default) 68 // Initial Conditions: 69 // Mass Density: 70 // Constant mass density of 1.0 71 // Momentum Density: 72 // Rotational field in x,y 73 // Energy Density: 74 // Maximum of 1. x0 decreasing linearly to 0. as radial distance 75 // increases to (1.-r/rc), then 0. everywhere else 76 // 77 // Boundary Conditions: 78 // Mass Density: 79 // 0.0 flux 80 // Momentum Density: 81 // 0.0 82 // Energy Density: 83 // 0.0 flux 84 // 85 // -- TRANSLATION 86 // Initial Conditions: 87 // Mass Density: 88 // Constant mass density of 1.0 89 // Momentum Density: 90 // Constant rectilinear field in x,y 91 // Energy Density: 92 // Maximum of 1. x0 decreasing linearly to 0. as radial distance 93 // increases to (1.-r/rc), then 0. everywhere else 94 // 95 // Boundary Conditions: 96 // Mass Density: 97 // 0.0 flux 98 // Momentum Density: 99 // 0.0 100 // Energy Density: 101 // Inflow BCs: 102 // E = E_wind 103 // Outflow BCs: 104 // E = E(boundary) 105 // Both In/Outflow BCs for E are applied weakly in the 106 // QFunction "Advection_Sur" 107 // 108 // ***************************************************************************** 109 110 // ***************************************************************************** 111 // This helper function provides support for the exact, time-dependent solution 112 // (currently not implemented) and IC formulation for 3D advection 113 // ***************************************************************************** 114 CEED_QFUNCTION_HELPER int Exact_Advection(CeedInt dim, CeedScalar time, 115 const CeedScalar X[], CeedInt Nf, CeedScalar q[], void *ctx) { 116 const SetupContext context = (SetupContext)ctx; 117 const CeedScalar rc = context->rc; 118 const CeedScalar lx = context->lx; 119 const CeedScalar ly = context->ly; 120 const CeedScalar lz = context->lz; 121 const CeedScalar *wind = context->wind; 122 123 // Setup 124 const CeedScalar x0[3] = {0.25*lx, 0.5*ly, 0.5*lz}; 125 const CeedScalar center[3] = {0.5*lx, 0.5*ly, 0.5*lz}; 126 127 // -- Coordinates 128 const CeedScalar x = X[0]; 129 const CeedScalar y = X[1]; 130 const CeedScalar z = X[2]; 131 132 // -- Energy 133 CeedScalar r = 0.; 134 switch (context->bubble_type) { 135 // original sphere 136 case 0: { // (dim=3) 137 r = sqrt(pow((x - x0[0]), 2) + 138 pow((y - x0[1]), 2) + 139 pow((z - x0[2]), 2)); 140 } break; 141 // cylinder (needs periodicity to work properly) 142 case 1: { // (dim=2) 143 r = sqrt(pow((x - x0[0]), 2) + 144 pow((y - x0[1]), 2) ); 145 } break; 146 } 147 148 // Initial Conditions 149 switch (context->wind_type) { 150 case 0: // Rotation 151 q[0] = 1.; 152 q[1] = -(y - center[1]); 153 q[2] = (x - center[0]); 154 q[3] = 0; 155 break; 156 case 1: // Translation 157 q[0] = 1.; 158 q[1] = wind[0]; 159 q[2] = wind[1]; 160 q[3] = wind[2]; 161 break; 162 } 163 164 switch (context->bubble_continuity_type) { 165 // original continuous, smooth shape 166 case 0: { 167 q[4] = r <= rc ? (1.-r/rc) : 0.; 168 } break; 169 // discontinuous, sharp back half shape 170 case 1: { 171 q[4] = ((r <= rc) && (y<center[1])) ? (1.-r/rc) : 0.; 172 } break; 173 // attempt to define a finite thickness that will get resolved under grid refinement 174 case 2: { 175 q[4] = ((r <= rc) 176 && (y<center[1])) ? (1.-r/rc)*fmin(1.0,(center[1]-y)/1.25) : 0.; 177 } break; 178 } 179 return 0; 180 } 181 182 // ***************************************************************************** 183 // This QFunction sets the initial conditions for 3D advection 184 // ***************************************************************************** 185 CEED_QFUNCTION(ICsAdvection)(void *ctx, CeedInt Q, 186 const CeedScalar *const *in, 187 CeedScalar *const *out) { 188 // Inputs 189 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 190 // Outputs 191 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 192 193 CeedPragmaSIMD 194 // Quadrature Point Loop 195 for (CeedInt i=0; i<Q; i++) { 196 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 197 CeedScalar q[5] = {}; 198 199 Exact_Advection(3, 0., x, 5, q, ctx); 200 for (CeedInt j=0; j<5; j++) q0[j][i] = q[j]; 201 } // End of Quadrature Point Loop 202 203 // Return 204 return 0; 205 } 206 207 // ***************************************************************************** 208 // This QFunction implements the following formulation of the advection equation 209 // 210 // This is 3D advection given in two formulations based upon the weak form. 211 // 212 // State Variables: q = ( rho, U1, U2, U3, E ) 213 // rho - Mass Density 214 // Ui - Momentum Density , Ui = rho ui 215 // E - Total Energy Density 216 // 217 // Advection Equation: 218 // dE/dt + div( E u ) = 0 219 // 220 // ***************************************************************************** 221 CEED_QFUNCTION(Advection)(void *ctx, CeedInt Q, 222 const CeedScalar *const *in, CeedScalar *const *out) { 223 // Inputs 224 // *INDENT-OFF* 225 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 226 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 227 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 228 229 // Outputs 230 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 231 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 232 // *INDENT-ON* 233 234 // Context 235 AdvectionContext context = (AdvectionContext)ctx; 236 const CeedScalar CtauS = context->CtauS; 237 const CeedScalar strong_form = context->strong_form; 238 239 CeedPragmaSIMD 240 // Quadrature Point Loop 241 for (CeedInt i=0; i<Q; i++) { 242 // Setup 243 // -- Interp in 244 const CeedScalar rho = q[0][i]; 245 const CeedScalar u[3] = {q[1][i] / rho, 246 q[2][i] / rho, 247 q[3][i] / rho 248 }; 249 const CeedScalar E = q[4][i]; 250 // -- Grad in 251 const CeedScalar drho[3] = {dq[0][0][i], 252 dq[1][0][i], 253 dq[2][0][i] 254 }; 255 // *INDENT-OFF* 256 const CeedScalar du[3][3] = {{(dq[0][1][i] - drho[0]*u[0]) / rho, 257 (dq[1][1][i] - drho[1]*u[0]) / rho, 258 (dq[2][1][i] - drho[2]*u[0]) / rho}, 259 {(dq[0][2][i] - drho[0]*u[1]) / rho, 260 (dq[1][2][i] - drho[1]*u[1]) / rho, 261 (dq[2][2][i] - drho[2]*u[1]) / rho}, 262 {(dq[0][3][i] - drho[0]*u[2]) / rho, 263 (dq[1][3][i] - drho[1]*u[2]) / rho, 264 (dq[2][3][i] - drho[2]*u[2]) / rho} 265 }; 266 // *INDENT-ON* 267 const CeedScalar dE[3] = {dq[0][4][i], 268 dq[1][4][i], 269 dq[2][4][i] 270 }; 271 // -- Interp-to-Interp q_data 272 const CeedScalar wdetJ = q_data[0][i]; 273 // -- Interp-to-Grad q_data 274 // ---- Inverse of change of coordinate matrix: X_i,j 275 // *INDENT-OFF* 276 const CeedScalar dXdx[3][3] = {{q_data[1][i], 277 q_data[2][i], 278 q_data[3][i]}, 279 {q_data[4][i], 280 q_data[5][i], 281 q_data[6][i]}, 282 {q_data[7][i], 283 q_data[8][i], 284 q_data[9][i]} 285 }; 286 // *INDENT-ON* 287 // The Physics 288 // Note with the order that du was filled and the order that dXdx was filled 289 // du[j][k]= du_j / dX_K (note cap K to be clear this is u_{j,xi_k}) 290 // dXdx[k][j] = dX_K / dx_j 291 // X_K=Kth reference element coordinate (note cap X and K instead of xi_k} 292 // x_j and u_j are jth physical position and velocity components 293 294 // No Change in density or momentum 295 for (CeedInt f=0; f<4; f++) { 296 for (CeedInt j=0; j<3; j++) 297 dv[j][f][i] = 0; 298 v[f][i] = 0; 299 } 300 301 // -- Total Energy 302 // Evaluate the strong form using div(E u) = u . grad(E) + E div(u) 303 // or in index notation: (u_j E)_{,j} = u_j E_j + E u_{j,j} 304 CeedScalar div_u = 0, u_dot_grad_E = 0; 305 for (CeedInt j=0; j<3; j++) { 306 CeedScalar dEdx_j = 0; 307 for (CeedInt k=0; k<3; k++) { 308 div_u += du[j][k] * dXdx[k][j]; // u_{j,j} = u_{j,K} X_{K,j} 309 dEdx_j += dE[k] * dXdx[k][j]; 310 } 311 u_dot_grad_E += u[j] * dEdx_j; 312 } 313 CeedScalar strong_conv = E*div_u + u_dot_grad_E; 314 315 // Weak Galerkin convection term: dv \cdot (E u) 316 for (CeedInt j=0; j<3; j++) 317 dv[j][4][i] = (1 - strong_form) * wdetJ * E * (u[0]*dXdx[j][0] + 318 u[1]*dXdx[j][1] + 319 u[2]*dXdx[j][2]); 320 v[4][i] = 0; 321 322 // Strong Galerkin convection term: - v div(E u) 323 v[4][i] = -strong_form * wdetJ * strong_conv; 324 325 // Stabilization requires a measure of element transit time in the velocity 326 // field u. 327 CeedScalar uX[3]; 328 for (CeedInt j=0; j<3; 329 j++) uX[j] = dXdx[j][0]*u[0] + dXdx[j][1]*u[1] + dXdx[j][2]*u[2]; 330 const CeedScalar TauS = CtauS / sqrt(uX[0]*uX[0] + uX[1]*uX[1] + uX[2]*uX[2]); 331 for (CeedInt j=0; j<3; j++) 332 dv[j][4][i] -= wdetJ * TauS * strong_conv * uX[j]; 333 } // End Quadrature Point Loop 334 335 return 0; 336 } 337 338 // ***************************************************************************** 339 // This QFunction implements 3D (mentioned above) with 340 // implicit time stepping method 341 // 342 // ***************************************************************************** 343 CEED_QFUNCTION(IFunction_Advection)(void *ctx, CeedInt Q, 344 const CeedScalar *const *in, 345 CeedScalar *const *out) { 346 // *INDENT-OFF* 347 // Inputs 348 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 349 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 350 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 351 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 352 // Outputs 353 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 354 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 355 // *INDENT-ON* 356 AdvectionContext context = (AdvectionContext)ctx; 357 const CeedScalar CtauS = context->CtauS; 358 const CeedScalar strong_form = context->strong_form; 359 360 CeedPragmaSIMD 361 // Quadrature Point Loop 362 for (CeedInt i=0; i<Q; i++) { 363 // Setup 364 // -- Interp in 365 const CeedScalar rho = q[0][i]; 366 const CeedScalar u[3] = {q[1][i] / rho, 367 q[2][i] / rho, 368 q[3][i] / rho 369 }; 370 const CeedScalar E = q[4][i]; 371 // -- Grad in 372 const CeedScalar drho[3] = {dq[0][0][i], 373 dq[1][0][i], 374 dq[2][0][i] 375 }; 376 // *INDENT-OFF* 377 const CeedScalar du[3][3] = {{(dq[0][1][i] - drho[0]*u[0]) / rho, 378 (dq[1][1][i] - drho[1]*u[0]) / rho, 379 (dq[2][1][i] - drho[2]*u[0]) / rho}, 380 {(dq[0][2][i] - drho[0]*u[1]) / rho, 381 (dq[1][2][i] - drho[1]*u[1]) / rho, 382 (dq[2][2][i] - drho[2]*u[1]) / rho}, 383 {(dq[0][3][i] - drho[0]*u[2]) / rho, 384 (dq[1][3][i] - drho[1]*u[2]) / rho, 385 (dq[2][3][i] - drho[2]*u[2]) / rho} 386 }; 387 // *INDENT-ON* 388 const CeedScalar dE[3] = {dq[0][4][i], 389 dq[1][4][i], 390 dq[2][4][i] 391 }; 392 // -- Interp-to-Interp q_data 393 const CeedScalar wdetJ = q_data[0][i]; 394 // -- Interp-to-Grad q_data 395 // ---- Inverse of change of coordinate matrix: X_i,j 396 // *INDENT-OFF* 397 const CeedScalar dXdx[3][3] = {{q_data[1][i], 398 q_data[2][i], 399 q_data[3][i]}, 400 {q_data[4][i], 401 q_data[5][i], 402 q_data[6][i]}, 403 {q_data[7][i], 404 q_data[8][i], 405 q_data[9][i]} 406 }; 407 // *INDENT-ON* 408 // The Physics 409 // Note with the order that du was filled and the order that dXdx was filled 410 // du[j][k]= du_j / dX_K (note cap K to be clear this is u_{j,xi_k} ) 411 // dXdx[k][j] = dX_K / dx_j 412 // X_K=Kth reference element coordinate (note cap X and K instead of xi_k} 413 // x_j and u_j are jth physical position and velocity components 414 415 // No Change in density or momentum 416 for (CeedInt f=0; f<4; f++) { 417 for (CeedInt j=0; j<3; j++) 418 dv[j][f][i] = 0; 419 v[f][i] = wdetJ * q_dot[f][i]; //K Mass/transient term 420 } 421 422 // -- Total Energy 423 // Evaluate the strong form using div(E u) = u . grad(E) + E div(u) 424 // or in index notation: (u_j E)_{,j} = u_j E_j + E u_{j,j} 425 CeedScalar div_u = 0, u_dot_grad_E = 0; 426 for (CeedInt j=0; j<3; j++) { 427 CeedScalar dEdx_j = 0; 428 for (CeedInt k=0; k<3; k++) { 429 div_u += du[j][k] * dXdx[k][j]; // u_{j,j} = u_{j,K} X_{K,j} 430 dEdx_j += dE[k] * dXdx[k][j]; 431 } 432 u_dot_grad_E += u[j] * dEdx_j; 433 } 434 CeedScalar strong_conv = E*div_u + u_dot_grad_E; 435 CeedScalar strong_res = q_dot[4][i] + strong_conv; 436 437 v[4][i] = wdetJ * q_dot[4][i]; // transient part (ALWAYS) 438 439 // Weak Galerkin convection term: -dv \cdot (E u) 440 for (CeedInt j=0; j<3; j++) 441 dv[j][4][i] = -wdetJ * (1 - strong_form) * E * (u[0]*dXdx[j][0] + 442 u[1]*dXdx[j][1] + 443 u[2]*dXdx[j][2]); 444 445 // Strong Galerkin convection term: v div(E u) 446 v[4][i] += wdetJ * strong_form * strong_conv; 447 448 // Stabilization requires a measure of element transit time in the velocity 449 // field u. 450 CeedScalar uX[3]; 451 for (CeedInt j=0; j<3; 452 j++) uX[j] = dXdx[j][0]*u[0] + dXdx[j][1]*u[1] + dXdx[j][2]*u[2]; 453 const CeedScalar TauS = CtauS / sqrt(uX[0]*uX[0] + uX[1]*uX[1] + uX[2]*uX[2]); 454 455 for (CeedInt j=0; j<3; j++) 456 switch (context->stabilization) { 457 case 0: 458 break; 459 case 1: dv[j][4][i] += wdetJ * TauS * strong_conv * uX[j]; //SU 460 break; 461 case 2: dv[j][4][i] += wdetJ * TauS * strong_res * uX[j]; //SUPG 462 break; 463 } 464 } // End Quadrature Point Loop 465 466 return 0; 467 } 468 469 // ***************************************************************************** 470 // This QFunction implements consistent outflow and inflow BCs 471 // for 3D advection 472 // 473 // Inflow and outflow faces are determined based on sign(dot(wind, normal)): 474 // sign(dot(wind, normal)) > 0 : outflow BCs 475 // sign(dot(wind, normal)) < 0 : inflow BCs 476 // 477 // Outflow BCs: 478 // The validity of the weak form of the governing equations is extended 479 // to the outflow and the current values of E are applied. 480 // 481 // Inflow BCs: 482 // A prescribed Total Energy (E_wind) is applied weakly. 483 // 484 // ***************************************************************************** 485 CEED_QFUNCTION(Advection_Sur)(void *ctx, CeedInt Q, 486 const CeedScalar *const *in, 487 CeedScalar *const *out) { 488 // *INDENT-OFF* 489 // Inputs 490 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 491 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 492 // Outputs 493 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 494 // *INDENT-ON* 495 AdvectionContext context = (AdvectionContext)ctx; 496 const CeedScalar E_wind = context->E_wind; 497 const CeedScalar strong_form = context->strong_form; 498 const bool implicit = context->implicit; 499 500 CeedPragmaSIMD 501 // Quadrature Point Loop 502 for (CeedInt i=0; i<Q; i++) { 503 // Setup 504 // -- Interp in 505 const CeedScalar rho = q[0][i]; 506 const CeedScalar u[3] = {q[1][i] / rho, 507 q[2][i] / rho, 508 q[3][i] / rho 509 }; 510 const CeedScalar E = q[4][i]; 511 512 // -- Interp-to-Interp q_data 513 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 514 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 515 // We can effect this by swapping the sign on this weight 516 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 517 518 // ---- Normal vectors 519 const CeedScalar norm[3] = {q_data_sur[1][i], 520 q_data_sur[2][i], 521 q_data_sur[3][i] 522 }; 523 // Normal velocity 524 const CeedScalar u_normal = norm[0]*u[0] + norm[1]*u[1] + norm[2]*u[2]; 525 526 // No Change in density or momentum 527 for (CeedInt j=0; j<4; j++) { 528 v[j][i] = 0; 529 } 530 // Implementing in/outflow BCs 531 if (u_normal > 0) { // outflow 532 v[4][i] = -(1 - strong_form) * wdetJb * E * u_normal; 533 } else { // inflow 534 v[4][i] = -(1 - strong_form) * wdetJb * E_wind * u_normal; 535 } 536 } // End Quadrature Point Loop 537 return 0; 538 } 539 // ***************************************************************************** 540 541 #endif // advection_h 542