Lines Matching +full:- +full:t
1 (example-petsc-navier-stokes)=
3 # Compressible Navier-Stokes mini-app
6 …-dependent Navier-Stokes equations of compressible gas dynamics in a static Eulerian three-dimensi…
7 Moreover, the Navier-Stokes example has been developed using PETSc, so that the pointwise physics (…
9 ## Running the mini-app
12 :start-after: <!-- fluids-inclusion -->
14 ## The Navier-Stokes equations
17 The compressible Navier-Stokes equations in conservative form are
21 \frac{\partial \rho}{\partial t} + \nabla \cdot \bm{U} &= 0 \\
22 …rac{\partial \bm{U}}{\partial t} + \nabla \cdot \left( \frac{\bm{U} \otimes \bm{U}}{\rho} + P \bm{…
23 …frac{\partial E}{\partial t} + \nabla \cdot \left( \frac{(E + P)\bm{U}}{\rho} -\bm{u} \cdot \bm{\s…
25 $$ (eq-ns)
27 …})^T + \lambda (\nabla \cdot \bm{u})\bm{I}_3)$ is the Cauchy (symmetric) stress tensor, with $\mu$…
28 …-ns`, $\rho$ represents the volume mass density, $U$ the momentum density (defined as $\bm{U}=\rho…
31 P = \left( {c_p}/{c_v} -1\right) \left( E - {\bm{U}\cdot\bm{U}}/{(2 \rho)} \right) \, ,
32 $$ (eq-state)
36 The system {eq}`eq-ns` can be rewritten in vector form
39 \frac{\partial \bm{q}}{\partial t} + \nabla \cdot \bm{F}(\bm{q}) -S(\bm{q}) = 0 \, ,
40 $$ (eq-vector-ns)
42 for the state variables 5-dimensional vector
60 - \bm{\sigma} \\
61 - \bm{u} \cdot \bm{\sigma} - k \nabla T
70 $$ (eq-ns-flux)
77 \bm{q}_N (\bm{x},t)^{(e)} = \sum_{k=1}^{P}\psi_k (\bm{x})\bm{q}_k^{(e)}
81 We use tensor-product bases $\psi_{kji} = h_i(X_0)h_j(X_1)h_k(X_2)$.
83 To obtain a finite element discretization, we first multiply the strong form {eq}`eq-vector-ns` by …
86 …mega} \bm v \cdot \left(\frac{\partial \bm{q}_N}{\partial t} + \nabla \cdot \bm{F}(\bm{q}_N) - \bm…
95 \int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t} - \bm{S}(\bm{q}_N) \right) \…
96 - \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\
100 $$ (eq-weak-vector-ns)
111 #### Explicit time-stepping method
113 …formulation is solved with the adaptive Runge-Kutta-Fehlberg (RKF4-5) method by default (any expli…
116 \bm{q}_N^{n+1} = \bm{q}_N^n + \Delta t \sum_{i=1}^{s} b_i k_i \, ,
123 k_1 &= f(t^n, \bm{q}_N^n)\\
124 k_2 &= f(t^n + c_2 \Delta t, \bm{q}_N^n + \Delta t (a_{21} k_1))\\
125 k_3 &= f(t^n + c_3 \Delta t, \bm{q}_N^n + \Delta t (a_{31} k_1 + a_{32} k_2))\\
127 k_i &= f\left(t^n + c_i \Delta t, \bm{q}_N^n + \Delta t \sum_{j=1}^s a_{ij} k_j \right)\\
134 f(t^n, \bm{q}_N^n) = - [\nabla \cdot \bm{F}(\bm{q}_N)]^n + [S(\bm{q}_N)]^n \, .
137 #### Implicit time-stepping method
139 …d using the option `-implicit` is solved with Backward Differentiation Formula (BDF) method by def…
143 \bm f(\bm q_N) \equiv \bm g(t^{n+1}, \bm{q}_N, \bm{\dot{q}}_N) = 0 \, ,
144 $$ (eq-ts-implicit-ns)
152 …e integration scheme (backward difference formulas, generalized alpha, implicit Runge-Kutta, etc.).
153 Each nonlinear system {eq}`eq-ts-implicit-ns` will correspond to a weak form, as explained below.
154 …ing how difficult a given problem is to solve, we consider the Jacobian of {eq}`eq-ts-implicit-ns`,
160 …$\Delta t$, so small time steps result in the Jacobian being dominated by the second term, which i…
162 …Both terms are significant for time-accurate simulation and the setup costs of strong precondition…
167 We solve {eq}`eq-weak-vector-ns` using a Galerkin discretization (default) or a stabilized method, …
169 …-dominated problems (any time the cell Péclet number is larger than 1), and those tend to blow up …
170 …hesetal2010`, which offers a comprehensive review of stabilization and shock-capturing methods for…
172 - **SUPG** (streamline-upwind/Petrov-Galerkin)
174 …ighted residual of the strong form {eq}`eq-vector-ns` is added to the Galerkin formulation {eq}`eq…
179 …\int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t} - \bm{S}(\bm{q}_N) \right) …
180 - \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\
182 … F_{\text{adv}}}{\partial \bm q}\right) \bm\tau \left( \frac{\partial \bm{q}_N}{\partial t} \, + \,
183 \nabla \cdot \bm{F} \, (\bm{q}_N) - \bm{S}(\bm{q}_N) \right) \,dV &= 0
186 $$ (eq-weak-vector-ns-supg)
188 This stabilization technique can be selected using the option `-stab supg`.
190 - **SU** (streamline-upwind)
192 …This method is a simplified version of *SUPG* {eq}`eq-weak-vector-ns-supg` which is developed for …
196 …\int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t} - \bm{S}(\bm{q}_N) \right) …
197 - \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\
202 $$ (eq-weak-vector-ns-su)
204 This stabilization technique can be selected using the option `-stab su`.
206 In both {eq}`eq-weak-vector-ns-su` and {eq}`eq-weak-vector-ns-supg`, $\bm\tau \in \mathbb R^{5\time…
207 … be explained via an ansatz for subgrid state fluctuations $\tilde{\bm q} = -\bm\tau \bm r$ where …
208 …riational form can be readily expressed by differentiating $\bm F_{\text{adv}}$ of {eq}`eq-ns-flux`
215 (\diff\bm U \otimes \bm U + \bm U \otimes \diff\bm U)/\rho - (\bm U \otimes \bm U)/\rho^2 \diff\rho…
216 (E + P)\diff\bm U/\rho + (\diff E + \diff P)\bm U/\rho - (E + P) \bm U/\rho^2 \diff\rho
221 where $\diff P$ is defined by differentiating {eq}`eq-state`.
224 …m X} = \nabla_{\bm x}\bm X \cdot \bm u$, with units of reference length (non-dimensional) per seco…
228 …it vector $\hat{\bm n}$ is given by $\lVert \bigl(\nabla_{\bm X} \bm x\bigr)^T \hat{\bm n} \rVert$.
237 $$ (eq-peclet)
239 For scalar advection-diffusion, the stabilization is a scalar
243 $$ (eq-tau-advdiff)
245 where $\xi(\mathrm{Pe}) = \coth \mathrm{Pe} - 1/\mathrm{Pe}$ approaches 1 at large local Péclet num…
246 Note that $\tau$ has units of time and, in the transport-dominated limit, is proportional to elemen…
247 For advection-diffusion, $\bm F(q) = \bm u q$, and thus the SU stabilization term is
251 $$ (eq-su-stabilize-advdiff)
253 where the term in parentheses is a rank-1 diffusivity tensor that has been pulled back to the refer…
254 See {cite}`hughesetal2010` equations 15-17 and 34-36 for further discussion of this formulation.
256 For the Navier-Stokes and Euler equations, {cite}`whiting2003hierarchical` defines a $5\times 5$ di…
261 The Navier-Stokes code in this example uses the following formulation for $\tau_c$, $\tau_m$, $\tau…
273 \mathcal{F} = \sqrt{ \rho^2 \left [ \left(\frac{2C_t}{\Delta t}\right)^2
278 where $\bm g = \nabla_{\bm x} \bm{X}^T \cdot \nabla_{\bm x} \bm{X}$ is the metric tensor and $\Vert…
281 For Advection-Diffusion, we use a modified version of the formulation for Navier-Stokes:
284 \tau = \left [ \left(\frac{2 C_t}{\Delta t}\right)^2
286 + \frac{\kappa^2 \Vert \bm g \Vert_F ^2}{C_d} \right]^{-1/2}
289 Otherwise, $C_a$ is set via `-Ctau_a` and $C_t$ via `-Ctau_t`.
295 $$ (eq-tau-conservative)
302 \Lambda_i = [u_i - a, u_i, u_i, u_i, u_i+a],
303 $$ (eq-eigval-advdiff)
311 $$ (eq-wavespeed)
317 …three types of problems/physical models that can be selected at run time via the option `-problem`.
318 …-advection`, the problem of the transport of energy in a uniform vector velocity field, {ref}`prob…
321 For scale-resolving simulations (such as LES and DNS), statistics for a simulation are more often u…
328 …gle \phi \rangle(x,y) = \frac{1}{L_z + (T_f - T_0)}\int_0^{L_z} \int_{T_0}^{T_f} \phi(x, y, z, t) …
337 The function $\langle \phi \rangle (x,y)$ is represented on a 2-D finite element grid, taken from t…
345 To represent these higher-order functions on the parent FEM space, we perform an $L^2$ projection.
349 \langle \phi \rangle_z(x,y,t) = \frac{1}{L_z} \int_0^{L_z} \phi(x, y, z, t) \mathrm{d}z
362 …_N = \int_0^{L_x} \int_0^{L_y} \left [\frac{1}{L_z} \int_0^{L_z} \phi(x,y,z,t) \mathrm{d}z \right …
368 \bm M [\langle \phi \rangle_z]_N = \frac{1}{L_z} \int_\Omega \phi(x,y,z,t) \psi^\mathrm{parent}_N(x…
377 To calculate the temporal integral, we do a running average using left-rectangle rule.
379 Periodically, the integral is updated using left-rectangle rule:
381 $$\overline{\phi}_\mathrm{new} = \overline{\phi}_{\mathrm{old}} + \phi(t_\mathrm{new}) \Delta T$$
382 where $\phi(t_\mathrm{new})$ is the statistic at the current time and $\Delta T$ is the time since …
383 When stats are written out to file, this running sum is then divided by $T_f - T_0$ to get the time…
388 … \rangle]_N = \frac{1}{L_z + (T_f - T_0)} \int_\Omega \int_{T_0}^{T_f} \phi(x,y,z,t) \psi^\mathrm{…
390 where the integral $\int_{T_0}^{T_f} \phi(x,y,z,t) \mathrm{d}t$ is calculated on a running basis.
395 This running average is only updated at the interval specified by `-ts_monitor_turbulence_spanstats…
396 …is only solved when statistics are written to file, which is controlled by `-ts_monitor_turbulence…
406 | ----------------- | -------- |
411 | $\langle \rho T \rangle$ | MeanDensityTemperature |
412 | $\langle \rho T u_i \rangle$ | MeanDensityTemperatureFlux[$i$] |
420 To get second-order statistics from these terms, simply use the identity:
423 \langle \phi' \theta' \rangle = \langle \phi \theta \rangle - \langle \phi \rangle \langle \theta \…
426 (differential-filtering)=
434 \overline{\phi} - \nabla \cdot (\beta (\bm{D}\bm{\Delta})^2 \nabla \overline{\phi} ) = \phi
437 …lution field, $\bm{\Delta} \in \mathbb{R}^{3 \times 3}$ a symmetric positive-definite rank 2 tenso…
442 - \cancel{\int_{\partial \Omega} \beta v \nabla \overline \phi \cdot (\bm{D}\bm{\Delta})^2 \bm{\hat…
446 The boundary integral resulting from integration-by-parts is crossed out, as we assume that $(\bm{D…
466 This is set via `-diff_filter_grid_based_width`.
468 … filter width tensor is most conveniently defined by $\bm{\Delta} = \bm{g}^{-1/2}$ where $\bm g = …
472 The coefficients for that anisotropic scaling are given by `-diff_filter_width_scaling`, denoted he…
490 \zeta = 1 - \exp\left(-\frac{y^+}{A^+}\right)
493 where $y^+$ is the wall-friction scaled wall-distance ($y^+ = y u_\tau / \nu = y/\delta_\nu$), $A^+…
494 … we assume that $\delta_\nu$ is constant across the wall and is defined by `-diff_filter_friction_…
495 $A^+$ is defined by `-diff_filter_damping_constant`.
497 To apply this scalar damping coefficient to the filter width tensor, we construct the wall-damping …
499 The wall-normal filter width is allowed to be damped to a zero filter width.
500 It is currently assumed that the second component of the filter width tensor is in the wall-normal …
517 $\beta$ can be set via `-diff_filter_kernel_scaling`.
519 (problem-advection)=
520 ## Advection-Diffusion
522 A simplified version of system {eq}`eq-ns`, only accounting for the transport of total energy, is g…
525 \frac{\partial E}{\partial t} + \nabla \cdot (\bm{u} E ) - \kappa \nabla E = 0 \, ,
526 $$ (eq-advection)
531 - **Rotation**
534 …We have solved {eq}`eq-advection` applying zero energy density $E$, and no-flux for $\bm{u}$ on th…
536 - **Translation**
540 …inflow boundaries such that the weak form boundary integral in {eq}`eq-weak-vector-ns` is defined …
547 …The weak form boundary integral in {eq}`eq-weak-vector-ns` for outflow boundary conditions is defi…
553 (problem-euler-vortex)=
557 Three-dimensional Euler equations, which are simplified and nondimensionalized version of system {e…
561 \frac{\partial \rho}{\partial t} + \nabla \cdot \bm{U} &= 0 \\
562 \frac{\partial \bm{U}}{\partial t} + \nabla \cdot \left( \frac{\bm{U} \otimes \bm{U}}{\rho} + P \bm…
563 \frac{\partial E}{\partial t} + \nabla \cdot \left( \frac{(E + P)\bm{U}}{\rho} \right) &= 0 \, , \\
565 $$ (eq-euler)
567 …=1$, $P=1$, $T=P/\rho= 1$ (Specific Gas Constant, $R$, is 1), and $\bm{u}=(u_1,u_2,0)$ while the p…
570 …2 \pi} \, e^{0.5(1-r^2)} \, (-\bar{y}, \, \bar{x}) \, , \\ \delta T &= - \frac{(\gamma-1) \, \epsi…
573 where $(\bar{x}, \, \bar{y}) = (x-x_c, \, y-y_c)$, $(x_c, \, y_c)$ represents the center of the dom…
576 (problem-shock-tube)=
580 …se for discontinuity capturing in one dimension. For this problem, the three-dimensional Euler equ…
603 h_{SHOCK} = 2 \left( C_{YZB} \,|\, \bm p \,|\, \right)^{-1}
614 (problem-density-current)=
617 …d to test non-reflecting/Riemann boundary conditions. It's primarily intended for Euler equations,…
623 \rho &= \rho_\infty\left(1+A\exp\left(\frac{-(\bar{x}^2 + \bar{y}^2)}{2\sigma^2}\right)\right) \\
625 E &= \frac{p_\infty}{\gamma -1}\left(1+A\exp\left(\frac{-(\bar{x}^2 + \bar{y}^2)}{2\sigma^2}\right)…
629 …nd width of the perturbation, respectively, and $(\bar{x}, \bar{y}) = (x-x_e, y-y_e)$ is the dista…
632 …ing an HLL (Harten, Lax, van Leer) Riemann solver {cite}`toro2009` (option `-freestream_riemann hl…
633 … such as HLLC {cite}`toro2009` (option `-freestream_riemann hllc`, which is default), which is a l…
635 ## Vortex Shedding - Flow past Cylinder
637 …ameter $D=1$ is centered at $(0,0)$ in a computational domain $-4.5 \leq x \leq 15.5$, $-4.5 \leq …
641 At time $t=0$, this domain is subjected to freestream boundary conditions at the inflow (left) and …
642 … at the top and bottom boundaries $(y = \pm 4.5)$ (zero normal velocity component, zero heat-flux).
643 The cylinder wall is an adiabatic (no heat flux) no-slip boundary condition.
663 …em (from {cite}`straka1993numerical`), we solve the full Navier-Stokes equations {eq}`eq-ns`, for …
664 …defined in terms of the Exner pressure, $\pi(\bm{x},t)$, and potential temperature, $\theta(\bm{x}…
667 …o &= \frac{P_0}{( c_p - c_v)\theta(\bm{x},t)} \pi(\bm{x},t)^{\frac{c_v}{ c_p - c_v}} \, , \\ e &= …
671 For this problem, we have used no-slip and non-penetration boundary conditions for $\bm{u}$, and no…
678 $$ u_1 = u_{\max} \left [ 1 - \left ( \frac{x_2}{H}\right)^2 \right] \quad \quad u_2 = u_3 = 0$$
679 $$T = T_w \left [ 1 + \frac{Pr \hat{E}c}{3} \left \{1 - \left(\frac{x_2}{H}\right)^4 \right \} \ri…
680 $$p = p_0 - \frac{2\rho_0 u_{\max}^2 x_1}{Re_H H}$$
682 where $H$ is the channel half-height, $u_{\max}$ is the center velocity, $T_w$ is the temperature a…
684 Boundary conditions are periodic in the streamwise direction, and no-slip and non-penetration bound…
689 ### Laminar Boundary Layer - Blasius
696 allowed to float and temperature is set constant. At the outlet, a user-set
698 flux terms use interior solution values). The wall is a no-slip,
699 no-penetration, no-heat flux condition. The top of the domain is treated as an
715 {cite}`shurSTG2014`. Below follows a re-description of the formulation to match
721 \bm{u}(\bm{x}, t) = \bm{\overline{u}}(\bm{x}) + \bm{C}(\bm{x}) \cdot \bm{v}'
726 …qrt{q^n(\bm{x})} \bm{\sigma}^n \cos(\kappa^n \bm{d}^n \cdot \bm{\hat{x}}^n(\bm{x}, t) + \phi^n ) \\
727 \bm{\hat{x}}^n &= \left[(x - U_0 t)\max(2\kappa_{\min}/\kappa^n, 0.1) , y, z \right]^T
733 tensor $\bm{C}$ (such that $\bm{R} = \bm{CC}^T$ ), bulk velocity $U_0$,
748 \kappa^n = \kappa_{\min} (1 + \alpha)^{n-1} \ , \quad \forall n=1, 2, ... , N
754 …appa^n}{\sum^N_{n=1} E(\kappa^n)\Delta \kappa^n} \ ,\quad \Delta \kappa^n = \kappa^n - \kappa^{n-1}
760 f_\eta = \exp \left[-(12\kappa /\kappa_\eta)^2 \right], \quad
761 f_\mathrm{cut} = \exp \left( - \left [ \frac{4\max(\kappa-0.9\kappa_\mathrm{cut}, 0)}{\kappa_\mathr…
765 (\nu^3/\varepsilon)^{-1/4}$ with $\nu$ the kinematic viscosity and
776 or temperature using the the `-weakT` flag.
804 lt --Calc-->ke --Calc-->kn
805 y --Calc-->ke
811 u0 --Copy--> u0C[U0]
817 ubar --Copy--> ubarC;
818 y --Copy--> yC;
819 lt --Copy--> ltC;
820 eps --Copy--> epsC;
822 rand --Copy--> randC;
823 rand --> N --Calc--> kn;
824 Rij --Calc--> Cij[C_ij]
827 This is done once at runtime. The spatially-varying terms are then evaluated at
828 each quadrature point on-the-fly, either by interpolation (for $l_t$,
853 | ----------------- | -------- | -------------- | --------- |
873 S(\bm{q}) = -\sigma(\bm{x})\left.\frac{\partial \bm{q}}{\partial \bm{Y}}\right\rvert_{\bm{q}} \bm{Y…
876 …- P_\mathrm{ref}, \bm{0}, 0]^T$, and $\sigma(\bm{x})$ is a linear ramp starting at `-idl_start` wi…
877 The damping is defined in terms of a pressure-primitive anomaly $\bm Y'$ converted to conservative …
878 $P_\mathrm{ref}$ has a default value equal to `-reference_pressure` flag, with an optional flag `-i…
883 …y the nodal layout of the default, equispaced box mesh and are enabled via `-mesh_transform platem…
885 The angle of this tilt is controlled by `-platemesh_top_angle`.
892 are placed such that `-platemesh_Ndelta` elements are within
893 `-platemesh_refine_height` of the wall. They are placed such that the element
894 height matches a geometric growth ratio defined by `-platemesh_growth`. The
895 remaining elements are then distributed from `-platemesh_refine_height` to the
903 specified via `-platemesh_y_node_locs_path`. If this flag is given an empty
906 ## Taylor-Green Vortex
908 This problem is really just an initial condition, the [Taylor-Green Vortex](https://en.wikipedia.or…
913 v &= -V_0 \cos(\hat x) \sin(\hat y) \sin(\hat z) \\
923 This initial condition is traditionally given for the incompressible Navier-Stokes equations.
924 …e reference state is selected using the `-reference_{velocity,pressure,temperature}` flags (Euclid…