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