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