1# Examples 2 3## Compressible Navier-Stokes Equations 4 5### Gaussian Wave 6This test case is taken/inspired by that presented in {cite}`mengaldoCompressibleBC2014`. It is intended to test non-reflecting/Riemann boundary conditions. It's primarily intended for Euler equations, but has been implemented for the Navier-Stokes equations here for flexibility. 7 8The problem has a perturbed initial condition and lets it evolve in time. The initial condition contains a Gaussian perturbation in the pressure field: 9 10$$ 11\begin{aligned} 12\rho &= \rho_\infty\left(1+A\exp\left(\frac{-(\bar{x}^2 + \bar{y}^2)}{2\sigma^2}\right)\right) \\ 13\bm{U} &= \bm U_\infty \\ 14E &= \frac{p_\infty}{\gamma -1}\left(1+A\exp\left(\frac{-(\bar{x}^2 + \bar{y}^2)}{2\sigma^2}\right)\right) + \frac{\bm U_\infty \cdot \bm U_\infty}{2\rho_\infty}, 15\end{aligned} 16$$ 17 18where $A$ and $\sigma$ are the amplitude and width of the perturbation, respectively, and $(\bar{x}, \bar{y}) = (x-x_e, y-y_e)$ is the distance to the epicenter of the perturbation, $(x_e, y_e)$. 19The simulation produces a strong acoustic wave and leaves behind a cold thermal bubble that advects at the fluid velocity. 20 21The boundary conditions are freestream in the x and y directions. When using an HLL (Harten, Lax, van Leer) Riemann solver {cite}`toro2009` (option `-freestream_riemann hll`), the acoustic waves exit the domain cleanly, but when the thermal bubble reaches the boundary, it produces strong thermal oscillations that become acoustic waves reflecting into the domain. 22This problem can be fixed using a more sophisticated Riemann solver such as HLLC {cite}`toro2009` (option `-freestream_riemann hllc`, which is default), which is a linear constant-pressure wave that transports temperature and transverse momentum at the fluid velocity. 23 24#### Running 25:::{list-table} Gaussian Wave Runtime Options 26:header-rows: 1 27 28* - Option 29 - Description 30 - Default value 31 - Unit 32 33* - `-epicenter` 34 - Coordinates of center of perturbation 35 - `0,0,0` 36 - `m` 37 38* - `-amplitude` 39 - Amplitude of the perturbation 40 - `0.1` 41 - 42 43* - `-width` 44 - Width parameter of the perturbation 45 - `0.002` 46 - `m` 47 48::: 49 50This problem can be run with the `examples/gaussianwave.yaml` file via: 51 52``` 53./build/navierstokes -options_file examples/gaussianwave.yaml 54``` 55 56`examples/gaussianwave.yaml`: 57```{literalinclude} ../examples/gaussianwave.yaml 58:language: yaml 59``` 60 61### Vortex Shedding - Flow past Cylinder 62This test case, based on {cite}`shakib1991femcfd`, is an example of using an externally provided mesh from Gmsh. 63A cylinder with diameter $D=1$ is centered at $(0,0)$ in a computational domain $-4.5 \leq x \leq 15.5$, $-4.5 \leq y \leq 4.5$. 64We solve this as a 3D problem with (default) one element in the $z$ direction. 65The domain is filled with an ideal gas at rest (zero velocity) with temperature 24.92 and pressure 7143. 66The viscosity is 0.01 and thermal conductivity is 14.34 to maintain a Prandtl number of 0.71, which is typical for air. 67At time $t=0$, this domain is subjected to freestream boundary conditions at the inflow (left) and Riemann-type outflow on the right, with exterior reference state at velocity $(1, 0, 0)$ giving Reynolds number $100$ and Mach number $0.01$. 68A symmetry (adiabatic free slip) condition is imposed at the top and bottom boundaries $(y = \pm 4.5)$ (zero normal velocity component, zero heat-flux). 69The cylinder wall is an adiabatic (no heat flux) no-slip boundary condition. 70As we evolve in time, eddies appear past the cylinder leading to a vortex shedding known as the vortex street, with shedding period of about 6. 71 72The Gmsh input file, `examples/meshes/cylinder.geo` is parametrized to facilitate experimenting with similar configurations. 73The Strouhal number (nondimensional shedding frequency) is sensitive to the size of the computational domain and boundary conditions. 74 75Forces on the cylinder walls are computed using the "reaction force" method, which is variationally consistent with the volume operator. 76Given the force components $\bm F = (F_x, F_y, F_z)$ and surface area $S = \pi D L_z$ where $L_z$ is the spanwise extent of the domain, we define the coefficients of lift and drag as 77 78$$ 79\begin{aligned} 80C_L &= \frac{2 F_y}{\rho_\infty u_\infty^2 S} \\ 81C_D &= \frac{2 F_x}{\rho_\infty u_\infty^2 S} \\ 82\end{aligned} 83$$ 84 85where $\rho_\infty, u_\infty$ are the freestream (inflow) density and velocity respectively. 86 87#### Running 88 89The initial condition is taken from `-reference_temperature` and `-reference_pressure`. 90To run this problem, first generate a mesh: 91 92```console 93$ make -C examples/meshes 94``` 95 96Then run by building the executable and running: 97 98```console 99$ make -j 100$ mpiexec -n 6 build/navierstokes -options_file examples/vortexshedding.yaml 101``` 102 103The vortex shedding period is roughly 5.6 and this problem runs until time 100 (2000 time steps). 104The above run writes a file named `force.csv` (see `ts_monitor_wall_force` in `examples/vortexshedding.yaml`), which can be postprocessed by running to create a figure showing lift and drag coefficients over time. 105 106```console 107$ python postprocess/vortexshedding.py 108``` 109 110`examples/vortexshedding.yaml`: 111```{literalinclude} ../examples/vortexshedding.yaml 112:language: yaml 113``` 114 115(example-density-current)= 116### Density Current 117 118For this test problem (from {cite}`straka1993numerical`), we solve the full Navier-Stokes equations {eq}`eq-ns`, for which a cold air bubble (of radius $r_c$) drops by convection in a neutrally stratified atmosphere. 119Its initial condition is defined in terms of the Exner pressure, $\pi(\bm{x},t)$, and potential temperature, $\theta(\bm{x},t)$, that relate to the state variables via 120 121$$ 122\begin{aligned} \rho &= \frac{P_0}{( c_p - c_v)\theta(\bm{x},t)} \pi(\bm{x},t)^{\frac{c_v}{ c_p - c_v}} \, , \\ e &= c_v \theta(\bm{x},t) \pi(\bm{x},t) + \bm{u}\cdot \bm{u} /2 + g z \, , \end{aligned} 123$$ 124 125where $P_0$ is the atmospheric pressure. 126For this problem, we have used no-slip and non-penetration boundary conditions for $\bm{u}$, and no-flux for mass and energy densities. 127 128#### Running 129 130:::{list-table} Density Current Runtime Options 131:header-rows: 1 132 133* - Option 134 - Description 135 - Default value 136 - Unit 137 138* - `-center` 139 - Location of bubble center 140 - `(lx,ly,lz)/2` 141 - `(m,m,m)` 142 143* - `-dc_axis` 144 - Axis of density current cylindrical anomaly, or `(0,0,0)` for spherically symmetric 145 - `(0,0,0)` 146 - 147 148* - `-rc` 149 - Characteristic radius of thermal bubble 150 - `1000` 151 - `m` 152 153* - `-theta0` 154 - Reference potential temperature 155 - `300` 156 - `K` 157 158* - `-thetaC` 159 - Perturbation of potential temperature 160 - `-15` 161 - `K` 162 163* - `-P0` 164 - Atmospheric pressure 165 - `1E5` 166 - `Pa` 167 168* - `-N` 169 - Brunt-Vaisala frequency 170 - `0.01` 171 - `1/s` 172::: 173 174This example can be run with: 175 176``` 177./navierstokes -problem density_current -dm_plex_box_faces 16,1,8 -degree 1 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 2000,125,1000 -dm_plex_dim 3 -rc 400. -bc_wall 1,2,5,6 -wall_comps 1,2,3 -bc_symmetry_y 3,4 -mu 75 178``` 179 180### Channel 181 182A compressible channel flow. Analytical solution given in 183{cite}`whitingStabilizedFEM1999`: 184 185$$ u_1 = u_{\max} \left [ 1 - \left ( \frac{x_2}{H}\right)^2 \right] \quad \quad u_2 = u_3 = 0$$ 186$$T = T_w \left [ 1 + \frac{Pr \hat{E}c}{3} \left \{1 - \left(\frac{x_2}{H}\right)^4 \right \} \right]$$ 187$$p = p_0 - \frac{2\rho_0 u_{\max}^2 x_1}{Re_H H}$$ 188 189where $H$ is the channel half-height, $u_{\max}$ is the center velocity, $T_w$ is the temperature at the wall, $Pr=\frac{\mu}{c_p \kappa}$ is the Prandlt number, $\hat E_c = \frac{u_{\max}^2}{c_p T_w}$ is the modified Eckert number, and $Re_h = \frac{u_{\max}H}{\nu}$ is the Reynolds number. 190 191Boundary conditions are periodic in the streamwise direction, and no-slip and non-penetration boundary conditions at the walls. 192The flow is driven by a body force determined analytically from the fluid properties and setup parameters $H$ and $u_{\max}$. 193 194#### Running 195 196:::{list-table} Channel Runtime Options 197:header-rows: 1 198 199* - Option 200 - Description 201 - Default value 202 - Unit 203 204* - `-umax` 205 - Maximum/centerline velocity of the flow 206 - `10` 207 - `m/s` 208 209* - `-theta0` 210 - Reference potential temperature 211 - `300` 212 - `K` 213 214* - `-P0` 215 - Atmospheric pressure 216 - `1E5` 217 - `Pa` 218 219* - `-body_force_scale` 220 - Multiplier for body force (`-1` for flow reversal) 221 - 1 222 - 223::: 224 225This problem can be run with the `examples/channel.yaml` file via: 226 227``` 228./build/navierstokes -options_file examples/channel.yaml 229``` 230`examples/channel.yaml:` 231```{literalinclude} ../examples/channel.yaml 232:language: yaml 233``` 234 235### Flat Plate Boundary Layer 236 237#### Meshing 238 239The flat plate boundary layer example has custom meshing features to better resolve the flow when using a generated box mesh. 240These meshing features modify the nodal layout of the default, equispaced box mesh and are enabled via `-mesh_transform platemesh`. 241One of those is tilting the top of the domain, allowing for it to be a outflow boundary condition. 242The angle of this tilt is controlled by `-platemesh_top_angle`. 243 244The primary meshing feature is the ability to grade the mesh, providing better resolution near the wall. 245There are two methods to do this; algorithmically, or specifying the node locations via a file. 246Algorithmically, a base node distribution is defined at the inlet (assumed to be $\min(x)$) and then linearly stretched/squeezed to match the slanted top boundary condition. 247Nodes are placed such that `-platemesh_Ndelta` elements are within `-platemesh_refine_height` of the wall. 248They are placed such that the element height matches a geometric growth ratio defined by `-platemesh_growth`. 249The remaining elements are then distributed from `-platemesh_refine_height` to the top of the domain linearly in logarithmic space. 250 251Alternatively, a file may be specified containing the locations of each node. 252The file should be newline delimited, with the first line specifying the number of points and the rest being the locations of the nodes. 253The node locations used exactly at the inlet (assumed to be $\min(x)$) and linearly stretched/squeezed to match the slanted top boundary condition. 254The file is specified via `-platemesh_y_node_locs_path`. 255If this flag is given an empty string, then the algorithmic approach will be performed. 256 257:::{list-table} Boundary Layer Meshing Runtime Options 258:header-rows: 1 259 260* - Option 261 - Description 262 - Default value 263 - Unit 264 265* - `-platemesh_modify_mesh` 266 - Whether to modify the mesh using the given options below. 267 - `false` 268 - 269 270* - `-platemesh_refine_height` 271 - Height at which `-platemesh_Ndelta` number of elements should refined into 272 - `5.9E-4` 273 - `m` 274 275* - `-platemesh_Ndelta` 276 - Number of elements to keep below `-platemesh_refine_height` 277 - `45` 278 - 279 280* - `-platemesh_growth` 281 - Growth rate of the elements in the refinement region 282 - `1.08` 283 - 284 285* - `-platemesh_top_angle` 286 - Downward angle of the top face of the domain. This face serves as an outlet. 287 - `5` 288 - `degrees` 289 290* - `-platemesh_y_node_locs_path` 291 - Path to file with y node locations. If empty, will use mesh warping instead. 292 - `""` 293 - 294::: 295 296(example-blasius)= 297#### Laminar Boundary Layer - Blasius 298 299Simulation of a laminar boundary layer flow, with the inflow being prescribed by a [Blasius similarity solution](https://en.wikipedia.org/wiki/Blasius_boundary_layer). 300At the inflow, the velocity is prescribed by the Blasius soution profile, density is set constant, and temperature is allowed to float. 301Using `weakT: true`, density is allowed to float and temperature is set constant. 302At the outlet, a user-set pressure is used for pressure in the inviscid flux terms (all other inviscid flux terms use interior solution values). 303The wall is a no-slip, no-penetration, no-heat flux condition. 304The top of the domain is treated as an outflow and is tilted at a downward angle to ensure that flow is always exiting it. 305 306#### Running 307:::{list-table} Blasius Runtime Options 308:header-rows: 1 309 310* - Option 311 - Description 312 - Default value 313 - Unit 314 315* - `-velocity_infinity` 316 - Freestream velocity 317 - `40` 318 - `m/s` 319 320* - `-temperature_infinity` 321 - Freestream temperature 322 - `288` 323 - `K` 324 325* - `-pressure_infinity` 326 - Atmospheric pressure, also sets IDL reference pressure 327 - `1.01E5` 328 - `Pa` 329 330* - `-temperature_wall` 331 - Wall temperature 332 - `288` 333 - `K` 334 335* - `-delta0` 336 - Boundary layer height at the inflow 337 - `4.2e-3` 338 - `m` 339 340* - `-n_chebyshev` 341 - Number of Chebyshev terms 342 - `20` 343 - 344 345* - `-chebyshev_` 346 - Prefix for Chebyshev snes solve 347 - 348 - 349::: 350 351This problem can be run with the `examples/blasius.yaml` file via: 352 353``` 354./build/navierstokes -options_file examples/blasius.yaml 355``` 356 357`examples/blasius.yaml`: 358```{literalinclude} ../examples/blasius.yaml 359:language: yaml 360``` 361 362#### Turbulent Boundary Layer 363 364Simulating a turbulent boundary layer without modeling the turbulence requires resolving the turbulent flow structures. 365These structures may be introduced into the simulations either by allowing a laminar boundary layer naturally transition to turbulence, or imposing turbulent structures at the inflow. 366The latter approach has been taken here, specifically using a *synthetic turbulence generation* (STG) method. 367See {ref}`bc-stg` for details on STG. 368 369#### Running 370This problem can be run with the `examples/blasius.yaml` file via: 371 372``` 373./build/navierstokes -options_file examples/blasius.yaml -stg_use true 374``` 375 376Note the added `-stg_use true` flag 377This overrides the `stg: use: false` setting in the `examples/blasius.yaml` file, enabling the use of the STG inflow. 378 379### Taylor-Green Vortex 380 381This problem is really just an initial condition, the [Taylor-Green Vortex](https://en.wikipedia.org/wiki/Taylor%E2%80%93Green_vortex): 382 383$$ 384\begin{aligned} 385u &= V'_x + V \sin(\hat x) \cos(\hat y) \sin(\hat z) \\ 386v &= V'_y - V \cos(\hat x) \sin(\hat y) \sin(\hat z) \\ 387w &= V'_z\\ 388p &= p_0 + \frac{\rho_0 V_0^2}{16} \left ( \cos(2 \hat x) + \cos(2 \hat y)\right) \left( \cos(2 \hat z) + 2 \right) \\ 389\rho &= \frac{p}{R T_0} \\ 390\end{aligned} 391$$ 392 393where $\hat x = 2 \pi x / L$ for $L$ the length of the domain in that specific direction, $V$ is the magnitude of the vortex velocity, and $V'_i$ is the background velocity that (optionally) advects the initial condition. 394The coordinate modification is done to transform a given grid onto a domain of $x,y,z \in [0, 2\pi)$. 395 396This initial condition is traditionally given for the incompressible Navier-Stokes equations. 397The reference state is selected using the `-reference_{velocity,pressure,temperature}` flags (Euclidean norm of `-reference_velocity` is used for $V$). 398 399:::{list-table} Taylor-Green Vortex Runtime Options 400:header-rows: 1 401 402* - Option 403 - Description 404 - Default value 405 - Unit 406 407* - `-taylorgreen_background_velocity` 408 - Background velocity added to Taylor-Green initial condition 409 - `0,0,0` 410 - `m/s` 411::: 412 413## Compressible Euler Equations 414 415(problem-euler-vortex)= 416 417### Isentropic Vortex 418 419Three-dimensional Euler equations, which are simplified and nondimensionalized version of system {eq}`eq-ns` and account only for the convective fluxes, are given by 420 421$$ 422\begin{aligned} 423\frac{\partial \rho}{\partial t} + \nabla \cdot \bm{U} &= 0 \\ 424\frac{\partial \bm{U}}{\partial t} + \nabla \cdot \left( \frac{\bm{U} \otimes \bm{U}}{\rho} + P \bm{I}_3 \right) &= 0 \\ 425\frac{\partial E}{\partial t} + \nabla \cdot \left( \frac{(E + P)\bm{U}}{\rho} \right) &= 0 \, , \\ 426\end{aligned} 427$$ (eq-euler) 428 429Following the setup given in {cite}`zhang2011verification`, the mean flow for this problem is $\rho=1$, $P=1$, $T=P/\rho= 1$ (Specific Gas Constant, $R$, is 1), and $\bm{u}=(u_1,u_2,0)$ while the perturbation $\delta \bm{u}$, and $\delta T$ are defined as 430 431$$ 432\begin{aligned} (\delta u_1, \, \delta u_2) &= \frac{\epsilon}{2 \pi} \, e^{0.5(1-r^2)} \, (-\bar{y}, \, \bar{x}) \, , \\ \delta T &= - \frac{(\gamma-1) \, \epsilon^2}{8 \, \gamma \, \pi^2} \, e^{1-r^2} \, , \\ \end{aligned} 433$$ 434 435where $(\bar{x}, \, \bar{y}) = (x-x_c, \, y-y_c)$, $(x_c, \, y_c)$ represents the center of the domain, $r^2=\bar{x}^2 + \bar{y}^2$, and $\epsilon$ is the vortex strength ($\epsilon$ < 10). 436There is no perturbation in the entropy $S=P/\rho^\gamma$ ($\delta S=0)$. 437 438#### Running 439:::{list-table} Isentropic Vortex Runtime Options 440:header-rows: 1 441 442* - Option 443 - Description 444 - Default value 445 - Unit 446 447* - `-center` 448 - Location of vortex center 449 - `(lx,ly,lz)/2` 450 - `(m,m,m)` 451 452* - `-mean_velocity` 453 - Background velocity vector 454 - `(1,1,0)` 455 - 456 457* - `-vortex_strength` 458 - Strength of vortex < 10 459 - `5` 460 - 461 462* - `-c_tau` 463 - Stabilization constant 464 - `0.5` 465 - 466::: 467 468This problem can be run with: 469 470``` 471./navierstokes -problem euler_vortex -dm_plex_box_faces 20,20,1 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1000,1000,50 -dm_plex_dim 3 -bc_inflow 4,6 -bc_outflow 3,5 -bc_symmetry_z 1,2 -mean_velocity .5,-.8,0. 472``` 473 474(problem-shock-tube)= 475### Shock Tube 476 477This test problem is based on Sod's Shock Tube (from{cite}`sodshocktubewiki`), a canonical test case for discontinuity capturing in one dimension. For this problem, the three-dimensional Euler equations are formulated exactly as in the Isentropic Vortex problem. The default initial conditions are $P=1$, $\rho=1$ for the driver section and $P=0.1$, $\rho=0.125$ for the driven section. The initial velocity is zero in both sections. Symmetry boundary conditions are applied to the side walls and wall boundary conditions are applied at the end walls. 478 479SU upwinding and discontinuity capturing have been implemented into the explicit timestepping operator for this problem. Discontinuity capturing is accomplished using a modified version of the $YZ\beta$ operator described in {cite}`tezduyar2007yzb`. This discontinuity capturing scheme involves the introduction of a dissipation term of the form 480 481$$ 482\int_{\Omega} \nu_{SHOCK} \nabla \bm v \!:\! \nabla \bm q dV 483$$ 484 485The shock capturing viscosity is implemented following the first formulation described in {cite}`tezduyar2007yzb`. The characteristic velocity $u_{cha}$ is taken to be the acoustic speed while the reference density $\rho_{ref}$ is just the local density. Shock capturing viscosity is defined by the following 486 487$$ 488\nu_{SHOCK} = \tau_{SHOCK} u_{cha}^2 489$$ 490 491where, 492 493$$ 494\tau_{SHOCK} = \frac{h_{SHOCK}}{2u_{cha}} \left( \frac{ \,|\, \nabla \rho \,|\, h_{SHOCK}}{\rho_{ref}} \right)^{\beta} 495$$ 496 497$\beta$ is a tuning parameter set between 1 (smoother shocks) and 2 (sharper shocks). The parameter $h_{SHOCK}$ is a length scale that is proportional to the element length in the direction of the density gradient unit vector. This density gradient unit vector is defined as $\hat{\bm j} = \frac{\nabla \rho}{|\nabla \rho|}$. The original formulation of Tezduyar and Senga relies on the shape function gradient to define the element length scale, but this gradient is not available to qFunctions in libCEED. To avoid this problem, $h_{SHOCK}$ is defined in the current implementation as 498 499$$ 500h_{SHOCK} = 2 \left( C_{YZB} \,|\, \bm p \,|\, \right)^{-1} 501$$ 502 503where 504 505$$ 506p_k = \hat{j}_i \frac{\partial \xi_i}{x_k} 507$$ 508 509The constant $C_{YZB}$ is set to 0.1 for piecewise linear elements in the current implementation. Larger values approaching unity are expected with more robust stabilization and implicit timestepping. 510 511#### Running 512:::{list-table} Shock Tube Runtime Options 513:header-rows: 1 514 515* - Option 516 - Description 517 - Default value 518 - Unit 519 520* - `-yzb` 521 - Use YZB discontinuity capturing 522 - `none` 523 - 524 525* - `-stab` 526 - Stabilization method (`none`, `su`, or `supg`) 527 - `none` 528 - 529 530* - `-c_tau` 531 - Stabilization constant 532 - `0.5` 533 - 534::: 535 536This problem can be run with: 537 538``` 539./navierstokes -problem shocktube -yzb -stab su -bc_symmetry_z 3,4 -bc_symmetry_y 1,2 -bc_wall 5,6 -dm_plex_dim 3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1000,100,100 -dm_plex_box_faces 200,1,1 -units_second 0.1 540``` 541 542(problem-advection)= 543## Advection-Diffusion 544 545There is a reduced mode in HONEE for pure advection-diffusion, which holds density $\rho$ and momentum density $\rho \bm u$ constant while advecting "total energy density" $E$. 546This reduced mode is given by 547 548$$ 549\frac{\partial E}{\partial t} + \nabla \cdot (\bm{u} E ) - \kappa \nabla E = 0 \, , 550$$ (eq-advection) 551 552with $\bm{u}$ the vector velocity field and $\kappa$ the diffusion coefficient. 553 554### Advection Field Options 555There are three different definitions for $\bm{u}$: 556 557- **Rotation** 558 559 A uniform circular velocity field transports the blob of total energy. 560 We have solved {eq}`eq-advection` applying zero energy density $E$, and no-flux for $\bm{u}$ on the boundaries. 561 562- **Translation** 563 564 In this case, a background wind with a constant rectilinear velocity field, enters the domain and transports the blob of total energy out of the domain. 565 566- **Boundary Layer** 567 568 This case has a linear velocity profile with only the y component set: $u_y = y / L_y$. 569 It starts at 0 for $y=0$ and then increases to 1 at the top of the domain. 570 571### Initial Condition options 572There are also several different definitions for initial conditions. 573Some require specific advection profiles, other's can be used with multiple. 574- **Bubble (Sphere and Cylinder)** 575 576 These are simple initial conditions with a controllable radius and center point. 577 They use an magnitude of one within the bubble and zero outside and have different smoothing options for the boundary of the bubble with the rest of the domain. 578 The difference between sphere and cylinder is whether the radius is applied in all 3 dimensions (sphere) or just in the x and y directions. 579 580- **Cosine Hill** 581 582 This is similar to the bubble ICs, but uses a cosine wave to define the bubble and it's radius is set to half the width of the domain (so the bubble fills the entire domain). 583 584- **Skew** 585 586 This IC is meant for for the translation advection profile only. 587 This IC features a line discontinuity intersecting the midpoint of the lower edge of the box and in the same direction as the advection velocity 588 The solution is either 0 or 1 on either side of the discontinuity. 589 590- **Wave** 591 592 This IC is meant for for the translation advection profile only. 593 This either a sine or square wave that oscillates in the direction of advection velocity. 594 The frequency and phase of the wave is controllable. 595 596- **Boundary Layer** 597 598 This IC is meant to be paired with the boundary layer advection profile. 599 This initial condition features a linear profile in the y direction up to a height set by the user. 600 601For the inflow boundary conditions, a prescribed $E_{wind}$ is applied weakly on the inflow boundaries such that the weak form boundary integral in {eq}`eq-weak-vector-ns` is defined as 602 603$$ 604\int_{\partial \Omega_{inflow}} \bm v \cdot \bm{F}(\bm q_N) \cdot \widehat{\bm{n}} \,dS = \int_{\partial \Omega_{inflow}} \bm v \, E_{wind} \, \bm u \cdot \widehat{\bm{n}} \,dS \, , 605$$ 606 607For the outflow boundary conditions, we have used the current values of $E$, following {cite}`papanastasiou1992outflow` which extends the validity of the weak form of the governing equations to the outflow instead of replacing them with unknown essential or natural boundary conditions. 608The weak form boundary integral in {eq}`eq-weak-vector-ns` for outflow boundary conditions is defined as 609 610$$ 611\int_{\partial \Omega_{outflow}} \bm v \cdot \bm{F}(\bm q_N) \cdot \widehat{\bm{n}} \,dS = \int_{\partial \Omega_{outflow}} \bm v \, E \, \bm u \cdot \widehat{\bm{n}} \,dS \, , 612$$ 613 614The advection problems can be run in both 2D and 3D, based on the DM defined for the problem. 615The following additional command-line options are available: 616 617### Running 618 619:::{list-table} Advection Runtime Options 620:header-rows: 1 621 622* - Option 623 - Description 624 - Default value 625 - Unit 626 627* - `-strong_form` 628 - Strong (1) or weak/integrated by parts (0) advection term of the residual 629 - `0` 630 - 631 632* - `-stab` 633 - Stabilization method (`none`, `su`, or `supg`) 634 - `none` 635 - 636 637* - `-stab_tau` 638 - Formulation for $\tau$ in stabilization (`ctau`, `advdiff_shakib`) 639 - `ctau` 640 - 641 642* - `-Ctau_t` 643 - Scaling factor on the temporal portion of the $\tau$ formulation 644 - 0. 645 - 646 647* - `-Ctau_a` 648 - Scaling factor on the advection portion of the $\tau$ formulation 649 - $P^2$ 650 - 651 652* - `-Ctau_d` 653 - Scaling factor on the diffusion portion of the $\tau$ formulation 654 - $P^4$ 655 - 656 657* - `-CtauS` 658 - Scale coefficient for stabilization tau (nondimensional) 659 - `0` 660 - 661 662* - `-diffusion_coeff` 663 - Diffusion coefficient 664 - `0` 665 - 666 667* - `-wind_type` 668 - Wind type in Advection (`rotation`, `translation`, `boundary_layer`) 669 - `rotation` 670 - 671 672* - `-wind_translation` 673 - Constant wind vector when `-wind_type translation` 674 - `1,0,0` 675 - 676 677* - `-E_wind` 678 - Total energy of inflow wind when `-wind_type translation` 679 - `1E6` 680 - `J` 681 682* - `-advection_ic_type` 683 - Initial condition type, (`sphere`, `cylinder`, `cosine_hill`, `skew`, `wave`, `boundary_layer`) 684 - `sphere` 685 - 686 687* - `-advection_ic_bubble_rc` 688 - For `sphere` or `cylinder` IC, characteristic radius of thermal bubble 689 - `1000` 690 - `m` 691 692* - `-advection_ic_bubble_continuity` 693 - For `sphere` or `cylinder` IC, different shapes of bubble, (`smooth`, `back_sharp`, `thick`, `cosine`) 694 - `smooth` 695 - 696 697* - `-advection_ic_wave_type` 698 - For `wave` IC, the wave form used for `-advection_ic_type wave` (`sine`, `square`) 699 - `sine` 700 - 701 702* - `-advection_ic_wave_frequency` 703 - For `wave` IC, frequency of the wave 704 - $2\pi$ 705 - `1/s` 706 707* - `-advection_ic_wave_phase` 708 - For `wave` IC, phase angle of the wave 709 - 0 710 - 711 712* - `-advection_ic_bl_height_factor` 713 - For `boundary_layer` IC, sets the height of the linear boundary layer initial condition in proportion to the domain height 714 - 1 715 - 716::: 717 718For 3D advection, an example of the `rotation` mode can be run with: 719 720``` 721./navierstokes -problem advection -dm_plex_box_faces 10,10,10 -dm_plex_dim 3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 8000,8000,8000 -bc_wall 1,2,3,4,5,6 -wall_comps 4 -wind_type rotation -implicit -stab su 722``` 723 724and the `translation` mode with: 725 726``` 727./navierstokes -problem advection -dm_plex_box_faces 10,10,10 -dm_plex_dim 3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 8000,8000,8000 -wind_type translation -wind_translation .5,-1,0 -bc_inflow 1,2,3,4,5,6 728``` 729 730For 2D advection, an example of the `rotation` mode can be run with: 731 732``` 733./navierstokes -problem advection -dm_plex_box_faces 20,20 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1000,1000 -bc_wall 1,2,3,4 -wall_comps 4 -wind_type rotation -implicit -stab supg 734``` 735 736and the `translation` mode with: 737 738``` 739./navierstokes -problem advection -dm_plex_box_faces 20,20 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1000,1000 -units_meter 1e-4 -wind_type translation -wind_translation 1,-.5 -bc_inflow 1,2,3,4 740``` 741Note the lengths in `-dm_plex_box_upper` are given in meters, and will be nondimensionalized according to `-units_meter`. 742 743 744The boundary layer problem can be run with: 745 746``` 747./build/navierstokes -options_file examples/advection_bl.yaml 748``` 749 750`examples/advection_bl.yaml`: 751```{literalinclude} ../examples/advection_bl.yaml 752:language: yaml 753``` 754 755The wave advection problem can be run with: 756 757``` 758./build/navierstokes -options_file examples/advection_wave.yaml 759``` 760 761`examples/advection_wave.yaml`: 762```{literalinclude} ../examples/advection_wave.yaml 763:language: yaml 764``` 765Note that the wave frequency, velocity direction, and domain size are set specifically to allow bi-periodic boundary conditions. 766