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