Problem setup

A problem is defined by a call to Setup, which precomputes all grid quantities:

setup = Setup(;
    x = (range(0.0, 2.0, 129), range(0.0, 1.0, 65)),
    boundary_conditions = (;
        u = (
            (DirichletBC((1.0, 0.0)), PressureBC()), # x-left, x-right
            (DirichletBC(), DirichletBC()),          # y-bottom, y-top
        ),
    ),
)

Grid

The tuple x contains one vector of volume boundary coordinates per dimension (so a 2D grid with N × M volumes needs vectors of length N + 1 and M + 1). Uniform grids are created with range; non-uniform grids can be created with stretched_grid, cosine_grid, or tanh_grid, or any custom monotone vector. The grid can be visualized with plotgrid.

Note that some functionality requires a uniform grid in one or more directions: the spectral pressure solvers and the spectral quantities require uniformity in the periodic directions.

Boundary conditions

The available boundary condition types are PeriodicBC, DirichletBC, SymmetricBC, and PressureBC. They are given per field, dimension, and side: boundary_conditions.u[α][s] is the boundary condition for the velocity in direction α on the low (s = 1) or high (s = 2) side. Each boundary has exactly one type of boundary condition, and for periodic boundary conditions the opposite boundary must also be periodic. For simulations with a temperature equation, add a temp field with boundary conditions for the temperature.

Boundary conditions are enforced by filling the ghost volumes of the field arrays, see Spatial and temporal discretization.

Initial conditions

Fields are allocated with scalarfield and vectorfield, and initialized with velocityfield and temperaturefield from a given analytical expression. For periodic boxes, random_field creates a random velocity field with a prescribed energy spectrum that is exactly divergence free on the staggered grid.

API

Setup and grid

IncompressibleNavierStokes.Dimension Type

struct Dimension{N}

Represent an N-dimensional space. Returns N when called.

julia> d = Dimension(3)
Dimension{3}()

julia> d()
3

Fields

source

IncompressibleNavierStokes.Setup Method

Setup(; x, boundary_conditions, backend, workgroupsize)

Create useful quantities for Cartesian box mesh ``x[1] \times \dots \times x[d]with boundary conditionsboundary_conditions. Return a named tuple ([α]` denotes a tuple index) with the following fields:

• N[α]: Number of finite volumes in direction α, including ghost volumes

• Nu[α][β]: Number of u[α] velocity DOFs in direction β

• Np[α]: Number of pressure DOFs in direction α

• Iu[α]: Cartesian index range of u[α] velocity DOFs

• Ip: Cartesian index range of pressure DOFs

• xlims[α]: Tuple containing the limits of the physical domain (not grid) in the direction α

• x[α]: α-coordinates of all volume boundaries, including the left point of the first ghost volume

• xu[α][β]: β-coordinates of u[α] velocity points

• xp[α]: α-coordinates of pressure points

• Δ[α]: All volume widths in direction α

• Δu[α]: Distance between pressure points in direction α

• A_stag[α]: Interpolation weights from volume centers to volume faces in the direction α

Note that the memory footprint of the redundant 1D-arrays above is negligible compared to the memory footprint of the 2D/3D-fields used in the code.

source

IncompressibleNavierStokes.cosine_grid Method

cosine_grid(a, b, N) -> Any

Create a nonuniform grid of N + 1 points from a to b using a cosine profile, i.e.

$$x_i=a+\frac{1}{2}(1-\cos \left(\pi \frac{i}{n}\right))(b-a),i=0,\dots ,N$$

See also stretched_grid.

source

IncompressibleNavierStokes.stretched_grid Function

stretched_grid(a, b, N) -> Any
stretched_grid(a, b, N, s) -> Any

Create a nonuniform grid of N + 1 points from a to b with a stretch factor of s. If s = 1, return a uniform spacing from a to b. Otherwise, return a vector $x\in {\mathbb{R}}^{N+1}$ such that $x_n=a+\sum _{i=1}^n{s}^{i-1}h$ for $n=0,\dots ,N$. Setting $x_N=b$ then gives $h=(b-a)\frac{1-s}{1-s^N}$, resulting in

$$x_n=a+(b-a)\frac{1-s^n}{1-s^N},n=0,\dots ,N.$$

Note that stretched_grid(a, b, N, s)[n] corresponds to $x_{n-1}$.

See also cosine_grid.

source

IncompressibleNavierStokes.tanh_grid Function

tanh_grid(a, b, N) -> Any
tanh_grid(a, b, N, γ) -> Any

Create a nonuniform grid of N + 1 points from a to b, as proposed by Trias et al. [16].

source

Boundary conditions

IncompressibleNavierStokes.AbstractBC Type

abstract type AbstractBC

Boundary condition for one side of the domain.

Fields

source

IncompressibleNavierStokes.DirichletBC Type

struct DirichletBC{U} <: IncompressibleNavierStokes.AbstractBC

Dirichlet boundary conditions for the velocity. The value u is one of:

• nothing (default): no-slip boundary conditions, where all velocity components are zero;

• a tuple of d constants (u1_BC, ..., ud_BC), one per velocity component, where d is the dimension;

• a function (dim, x..., t) -> u_BC returning the boundary value of velocity component dim at the point x... and time t.

Fields

• u: Boundary condition

source

IncompressibleNavierStokes.PeriodicBC Type

struct PeriodicBC <: IncompressibleNavierStokes.AbstractBC

Periodic boundary conditions. Must be periodic on both sides.

Fields

source

IncompressibleNavierStokes.PressureBC Type

struct PressureBC <: IncompressibleNavierStokes.AbstractBC

Pressure boundary conditions. The pressure is prescribed on the boundary (usually an "outlet"). The velocity has zero Neumann conditions.

Note: Currently, the pressure is prescribed with the constant value of zero on the entire boundary.

Fields

source

IncompressibleNavierStokes.SymmetricBC Type

struct SymmetricBC <: IncompressibleNavierStokes.AbstractBC

Symmetric boundary conditions. The parallel velocity and pressure is the same at each side of the boundary. The normal velocity is zero.

Fields

source

IncompressibleNavierStokes.apply_bc_p! Method

apply_bc_p!(p, t, setup; kwargs...) -> Any

Apply pressure boundary conditions (in-place version).

source

IncompressibleNavierStokes.apply_bc_p Method

apply_bc_p(p, t, setup; kwargs...) -> Any

Apply pressure boundary conditions (differentiable version).

source

IncompressibleNavierStokes.apply_bc_temp! Method

apply_bc_temp!(temp, t, setup; kwargs...) -> Any

Apply temperature boundary conditions (in-place version).

source

IncompressibleNavierStokes.apply_bc_temp Method

apply_bc_temp(temp, t, setup; kwargs...) -> Any

Apply temperature boundary conditions (differentiable version).

source

IncompressibleNavierStokes.apply_bc_u! Method

apply_bc_u!(u, t, setup; kwargs...) -> Any

Apply velocity boundary conditions (in-place version).

source

IncompressibleNavierStokes.apply_bc_u Method

apply_bc_u(u, t, setup; kwargs...) -> Any

Apply velocity boundary conditions (differentiable version).

source

IncompressibleNavierStokes.boundary Method

boundary(β, N, I, isright) -> Any

Get boundary indices of boundary layer normal to β. The CartesianIndices given by I should contain those of the inner DOFs, typically Ip or Iu[α]. The boundary layer is then just outside those.

source

IncompressibleNavierStokes.offset_p Function

offset_p(bc, isright)

Number of non-DOF pressure components at boundary. If isright, it is at the end/right/rear/top boundary, otherwise beginning.

source

IncompressibleNavierStokes.offset_u Function

offset_u(bc, isright, isnormal)

Number of non-DOF velocity components at boundary. If isnormal, then the velocity is normal to the boundary, else parallel. If isright, it is at the end/right/rear/top boundary, otherwise beginning.

source

Field initializers

IncompressibleNavierStokes.orlandi_profile Method

orlandi_profile(k; kpeak) -> Any

Initial energy spectrum profile $E(k)\propto k^4{e}^{-2(k/k_p{)}^2}$ with peak wavenumber kpeak, as in Orlandi [17].

source

IncompressibleNavierStokes.random_field Function

random_field(setup; ...) -> Any
random_field(
    setup,
    t;
    profile,
    totalenergy,
    rng,
    kwargs...
) -> Any

Create a random divergence-free velocity field with a prescribed energy spectrum profile. The energy in the wavenumber shell κ ≤ |k| < κ + 1 is totalenergy * profile(κ; kwargs...) / p, where p normalizes the profile such that the mean kinetic energy density ⟨u_i u_i⟩ / 2 is exactly totalenergy. By default, the Orlandi form orlandi_profile is used (pass e.g. kpeak = 5 to move the peak).

The field is constructed in spectral space from white noise, projected onto the divergence-free space of the staggered grid (using the modified wavenumbers of the staggered divergence operator), and rescaled shell-wise. The resulting field is thus exactly divergence-free on the staggered grid and has exactly the prescribed spectrum.

source

IncompressibleNavierStokes.scalarfield Method

scalarfield(setup) -> Any

Create empty scalar field.

source

IncompressibleNavierStokes.symmetric_tensorfield Method

symmetric_tensorfield(
    setup
) -> Union{NamedTuple{(:xx, :xy, :xz, :yy, :yz, :zz), <:NTuple{6, Any}}, NamedTuple{(:xx, :xy, :yy), <:Tuple{Any, Any, Any}}}

Symmetric tensor field, stored as a named tuple σ.ij.

source

IncompressibleNavierStokes.temperaturefield Function

temperaturefield(setup, tempfunc) -> Any
temperaturefield(setup, tempfunc, t) -> Any

Create temperature field from function with boundary conditions at time t.

source

IncompressibleNavierStokes.tensorfield Method

tensorfield(
    setup
) -> Union{NamedTuple{(:xx, :yx, :zx, :xy, :yy, :zy, :xz, :yz, :zz), <:NTuple{9, Any}}, NamedTuple{(:xx, :yx, :xy, :yy), <:NTuple{4, Any}}}

Non-symmetric tensor field, stored as a named tuple σ.ij.

source

IncompressibleNavierStokes.vectorfield Method

vectorfield(setup) -> Any

Create empty vector field.

source

IncompressibleNavierStokes.velocityfield Function

velocityfield(setup, ufunc; ...) -> Any
velocityfield(setup, ufunc, t; psolver, doproject) -> Any

Create divergence free velocity field u with boundary conditions at time t. The initial conditions of u[α] are specified by the function ufunc(α, x...).

source