.. _vc_taylor_green_vortex:

====================
Taylor-Green Vortex
====================

The Taylor-Green vortex (TGV) is a canonical benchmark used to assess a solver's ability
to reproduce unsteady, spatially periodic flows with known reference data.
The initial conditions take the form of smooth trigonometric functions, making the TGV
an ideal test for equation-based initialization: the solver must evaluate :math:`\sin`
and :math:`\cos` expressions at every lattice node at startup and reproduce them to
floating-point precision before the time-marching begins.

Two complementary sub-cases are used here.
Case 00a is a 2-D TGV, which admits a closed-form analytical solution valid for all
time and therefore provides a stringent quantitative accuracy and convergence test.
Case 00b is a 3-D TGV, for which no analytical solution exists once the flow has
transitioned to turbulence; the dissipation rate is compared against the
pseudo-spectral DNS of :footcite:t:`vanRees2011`.


.. _vc_tgv_2d:

2-D Taylor-Green Vortex (Case 00a)
===================================

The 2-D TGV is a doubly periodic flow on the square domain
:math:`[0,\, 2\pi] \times [0,\, 2\pi]` initialised with a single Fourier mode.
Because the nonlinear advection terms are exactly zero for this initial condition,
the flow decays purely through viscous diffusion and the velocity field remains
a rescaled copy of the initial condition at all times.
This property yields an exact analytical solution that the LBM result can be
compared against directly.

.. _vc_tgv_2d_ic:

Initial conditions
------------------

The velocity and density fields at :math:`t = 0` are:

.. math::
    u_x(x, y, 0) &= U_0 \cos(k x) \sin(k y) \\
    u_y(x, y, 0) &= -U_0 \sin(k x) \cos(k y) \\
    \rho(x, y, 0) &= \rho_0
    :label: vc_tgv2d_ic

where :math:`k = 2\pi / L` is the wave number, :math:`L` is the domain length, and
:math:`U_0` is the peak velocity amplitude.
For a lattice of :math:`N` nodes per side with lattice spacing :math:`\Delta x = 1`,
the wave number in lattice units is :math:`k = 2\pi / N`.

These fields are prescribed via the ``models.initialization.equations`` block introduced
in PR #443.  Setting ``use_equation_init: true`` and providing symbolic expressions for
:math:`u_x`, :math:`u_y`, and :math:`\rho` causes the solver to evaluate them on the
GPU before the first time step, placing the initial condition to machine precision
without any approximation from interpolation or file I/O.

.. _vc_tgv_2d_analytical:

Analytical solution
-------------------

The exact solution for :math:`t > 0` is an exponential decay:

.. math::
    u_x(x, y, t) &= U_0 \exp\!\left(-2\nu k^2 t\right) \cos(k x) \sin(k y) \\
    u_y(x, y, t) &= -U_0 \exp\!\left(-2\nu k^2 t\right) \sin(k x) \cos(k y)
    :label: vc_tgv2d_analytical

The spatially integrated kinetic energy per unit area is:

.. math::
    E(t) = \frac{U_0^2}{4} \exp\!\left(-4\nu k^2 t\right) = E_0 \exp\!\left(-4\nu k^2 t\right)
    :label: vc_tgv2d_energy

The characteristic decay time is :math:`t_\nu = (4\nu k^2)^{-1}`.
For a well-resolved simulation the solver result should track
:math:numref:`vc_tgv2d_analytical` with a velocity error that decreases as
:math:`O(\Delta x^2)` when the grid is refined, consistent with the second-order
accuracy of the LBM streaming step.

.. _vc_tgv_2d_params:

Simulation parameters
---------------------

Four grid resolutions are run to assess convergence.
All cases use the D2Q9 velocity set with the RR-BGK collision operator
(the production 3rd-order Hermite operator).

.. list-table::
   :header-rows: 1
   :widths: 35 65

   * - Parameter
     - Value
   * - Velocity set
     - D2Q9
   * - Collision operator
     - RRBGK
   * - Relaxation time :math:`\tau`
     - 0.501 - 0.508 (varies with grid size)
   * - Peak velocity :math:`U_0`
     - 0.01 (lattice units)
   * - Reference density :math:`\rho_0`
     - 1.0
   * - Reynolds number
     - 240 (:math:`Re = U_0 N / \nu`, constant across grids)
   * - Domain
     - :math:`[0,\, 2\pi] \times [0,\, 2\pi]`, periodic in both directions
   * - Grid sizes :math:`N`
     - 8, 16, 32, 64 nodes per side
   * - Boundary conditions
     - Fully periodic (``periodic_dims: [true, true]``)
   * - Initialization
     - Equation-based (``use_equation_init: true``)

The Mach number :math:`\mathrm{Ma} = U_0 / c_s = U_0 \sqrt{3} \approx 0.017` is well
below the weakly-compressible limit of 0.1, so compressibility errors are negligible.

.. note::

   The convergence study runs four independent simulations (N = 8, 16, 32, 64).
   The expected slope on a log-log plot of the :math:`L_2` velocity error versus
   :math:`\Delta x` is 2, confirming second-order spatial accuracy.


.. _vc_tgv_3d:

3-D Taylor-Green Vortex (Case 00b)
===================================

The 3-D TGV starts from a smooth, low-amplitude initial condition that is linearly
unstable.  At moderate Reynolds numbers the flow develops a cascade of smaller vortical
structures and eventually reaches a regime of decaying turbulence.
Because an exact solution does not exist for :math:`t > 0`, validation is performed by
comparing the kinetic energy dissipation rate :math:`\varepsilon(t) = -\mathrm{d}E/\mathrm{d}t`
against the spectral DNS results of :footcite:t:`vanRees2011`.

.. _vc_tgv_3d_ic:

Initial conditions
------------------

The velocity and pressure fields at :math:`t = 0` are the standard TGV initial condition
on the periodic cube :math:`[0,\, 2\pi]^3`:

.. math::
    u_x(x, y, z, 0) &= V_0 \sin(k x) \cos(k y) \cos(k z) \\
    u_y(x, y, z, 0) &= -V_0 \cos(k x) \sin(k y) \cos(k z) \\
    u_z(x, y, z, 0) &= 0 \\
    \rho(x, y, z, 0) &= \rho_0 + \frac{\rho_0 V_0^2}{16 c_s^2}
                        \left[\cos(2kz) + 2\right]
                        \left[\cos(2kx) + \cos(2ky)\right]
    :label: vc_tgv3d_ic

where :math:`k = 2\pi / L` and :math:`V_0` is the peak velocity.
The density expression in :math:numref:`vc_tgv3d_ic` is the weakly-compressible
approximation of the incompressible pressure distribution.
In the current simulation :math:`\rho_0 = 1` is used uniformly (i.e. the density
correction is omitted) because at :math:`\mathrm{Ma} \approx 0.07` the maximum
perturbation is of order :math:`10^{-3}` and does not affect the velocity dynamics.

As with Case 00a, these fields are set via ``models.initialization.equations``,
evaluating the trigonometric expressions on the GPU at startup.

.. _vc_tgv_3d_validation:

Validation approach
-------------------

The volume-averaged kinetic energy is:

.. math::
    E(t) = \frac{1}{N^3} \sum_{i,j,k} \frac{1}{2} \rho(i,j,k,t)
           \left(u_x^2 + u_y^2 + u_z^2\right)
    :label: vc_tgv3d_energy

The dissipation rate is estimated by finite difference:

.. math::
    \varepsilon(t) = -\frac{\mathrm{d}E}{\mathrm{d}t} \approx
                    -\frac{E(t + \Delta T) - E(t - \Delta T)}{2\,\Delta T}
    :label: vc_tgv3d_dissipation

where :math:`\Delta T` is the interval between saved snapshots.

The DNS reference of :footcite:t:`vanRees2011` was computed with a pseudospectral code
at :math:`Re_\Gamma = 1/\nu = 1600` on grids up to :math:`768^3`.
The reference dissipation curve is digitised from Figure 8a and 8b of the paper
(768\ :sup:`3` PS results, the most converged resolution).
The Reynolds number is defined with characteristic length :math:`L^* = 1/k_0 = 1`
(the inverse wavenumber), so in lattice units :math:`Re = V_0 N / (2\pi\nu)`.
The peak dissipation time and magnitude are well established and serve as the primary
quantitative targets.

.. _vc_tgv_3d_params:

Simulation parameters
---------------------

.. list-table::
   :header-rows: 1
   :widths: 35 65

   * - Parameter
     - Value
   * - Velocity set
     - D3Q27
   * - Collision operator
     - RRBGK
   * - Relaxation time :math:`\tau`
     - 0.503056
   * - Peak velocity :math:`V_0`
     - 0.04 (lattice units)
   * - Reference density :math:`\rho_0`
     - 1.0
   * - Domain
     - :math:`[0,\, 2\pi]^3`, periodic in all three directions
   * - Grid size :math:`N`
     - 256 nodes per side
   * - Boundary conditions
     - Fully periodic (``periodic_dims: [true, true, true]``)
   * - Initialization
     - Equation-based (``use_equation_init: true``)
   * - Reference DNS
     - :footcite:t:`vanRees2011`, :math:`Re_\Gamma = 1600`

.. note::

   The Reynolds number follows the van Rees convention :math:`Re_\Gamma = 1/\nu`
   with characteristic length :math:`L^* = 1` (inverse wavenumber) and
   :math:`V_0 = 1`. In lattice units this becomes
   :math:`Re = V_0 N / (2\pi\nu)`, giving :math:`\nu \approx 0.001019`
   and :math:`\tau = 3\nu + 0.5 \approx 0.503056`.


.. toctree::
   :maxdepth: 1
   :hidden:

   2-D Taylor-Green Vortex<00a_taylor_green_vortex_2d.ipynb>
   3-D Taylor-Green Vortex<00b_taylor_green_vortex_3d.ipynb>

.. footbibliography::