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