Remove Tullio dependency; take #2 and implement universally accurate CFL time-scale

[simone-silvestri]

Take two on PR #2252
i.e. use AbstractOperations to calculate the CFL instead of Tullio.jl

[glwagner]

@navidcy and I implemented this CFL operation a bit differently but I can't remember where ... maybe @navidcy can point you to it

[navidcy]

Oceananigans.jl/src/Advection/cell_advection_timescale.jl

Lines 1 to 18 in bcc34f0

	using Oceananigans.AbstractOperations: KernelFunctionOperation
	using Oceananigans.Grids: topology, min_Δx, min_Δy, min_Δz
	function cell_advection_timescale(grid, velocities)
	u, v, w = velocities
	τ = KernelFunctionOperation{Center, Center, Center}(cell_advection_timescaleᶜᶜᶜ, grid, u, v, w)
	return minimum(τ)
	end
	@inline function cell_advection_timescaleᶜᶜᶜ(i, j, k, grid, u, v, w)
	Δx = Δxᶠᶜᶜ(i, j, k, grid)
	Δy = Δyᶜᶠᶜ(i, j, k, grid)
	Δz = Δzᶜᶜᶠ(i, j, k, grid)
	return @inbounds min(Δx / abs(u[i, j, k]),
	Δy / abs(v[i, j, k]),
	Δz / abs(w[i, j, k]))
	end

[simone-silvestri]

Ok, it's doing the same thing, I'll just remove Tullio and use this for the CFL

[simone-silvestri]

I am going to change a bit the definition to be in line with what we had before, so instead of

@inline function cell_advection_timescaleᶜᶜᶜ(i, j, k, grid, u, v, w)
    Δx = Δxᶠᶜᶜ(i, j, k, grid)
    Δy = Δyᶜᶠᶜ(i, j, k, grid)
    Δz = Δzᶜᶜᶠ(i, j, k, grid)

    return @inbounds min(Δx / abs(u[i, j, k]),
                         Δy / abs(v[i, j, k]),
                         Δz / abs(w[i, j, k])) 
end

take all 3 dimensions

@inline function cell_advection_timescaleᶜᶜᶜ(i, j, k, grid, u, v, w)
    Δx = Δxᶠᶜᶜ(i, j, k, grid)
    Δy = Δyᶜᶠᶜ(i, j, k, grid)
    Δz = Δzᶜᶜᶠ(i, j, k, grid)

    return @inbounds 1 / (abs(u[i, j, k]) / Δx +
                          abs(v[i, j, k]) / Δy +
                          abs(w[i, j, k]) / Δz) 
end

it is a small difference but it is a more conservative definition (and usually used in n-dimensional cases)

[navidcy]

Looks good.

Do we want this cfl.jl to live in Diagnostics or in Advection?

[navidcy]

I guess there is diffusive cfl...?

[navidcy]

is this the proper way to compute CFL? I'll trust you on that... That's like the geometric mean of the timescales in each direction?

[simone-silvestri]

Imagine a 2D situation, we start from
$$\frac{\partial c}{\partial t} = u\frac{\partial c}{\partial x} + v\frac{\partial c}{\partial y}$$
We have advection in x and y, by extension from the 1D CFL condition we can approximate the CFL condition as
$$|u|\frac{\Delta t}{\Delta x} + |v|\frac{\Delta t}{\Delta y} < CFL (\text{ maximum 1 })$$
This means that
$$\Delta t / \tau < CFL$$
where
$$\tau = \left(\frac{|u|}{\Delta x} + \frac{|v|}{\Delta y} \right)^{-1}$$
The extension to 3D includes also the advection in z

This is just a choice, as to account for diagonal motion probably you can get away with a square root of 2 somewhere. But it is the usual definition of the CFL condition
https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition

[navidcy]

Sounds good! Why don't we include all this in the docstring of the cell_advection?

[simone-silvestri]

Ok, I was a bit unsure which CFL definition is the best, so I checked the 2D stability and I got a slightly better proof.
Let's start with the discretized version of the 2D advection equation (forward euler and 2nd order centered finite volume on a staggered regular cartesian grid)
$$\frac{c^{n+1} - c^n}{\Delta t} = u \left( \frac{c^n_i - c^n_{i-1}}{\Delta x}\right) + v \left( \frac{c^n_j - c^n_{j-1}}{\Delta y}\right)$$
Let's assume that c is a two-dimensional wave depending on an x wavenumber $\kappa$, a y wavenumber $\mathcal{l}$, and a time-dependent amplification factor $\xi(t)$, then
$c^n_{ij} = \xi^n \exp{(- \Im \kappa \cdot i \Delta x - \Im \mathcal{l} \cdot j \Delta y)}$. (because of overlap with the x-index $i$, I defined $\Im$ as the imaginary number $\Im = \sqrt{-1}$)
Substituting this definition of $c$ and dividing through by $c^n_{ij}$ we get
$$\frac{\xi^{n+1} / \xi^n - 1}{\Delta t} = u \left(\frac{1 - \exp{(- \Im \kappa \Delta x)}}{\Delta x}\right) + v \left(\frac{1 - \exp{(- \Im \mathcal{l} \Delta y)}}{\Delta y}\right)$$
we can make use of $\exp{\Im \theta} = \cos{\theta} + \Im \sin{\theta}$ and rewrite a bit:
$$\frac{\xi^{n+1}}{\xi^n} = 1 + \Delta t \cdot \left[ \frac{u}{\Delta x} \left( 1 - \cos{\kappa \Delta x} + \Im \sin{\kappa\Delta x}\right) + \frac{v}{\Delta y} \left( 1 - \cos{\mathcal{l} \Delta y} + \Im \sin{\mathcal{l} \Delta y} \right) \right]$$
Now, to ensure stability, the real part of $\xi^{n+1} / \xi^n$ should be bounded, so we have to ensure that
$$\left| \Re \left( \frac{\xi^{n+1}}{\xi^n} \right) \right| < 1$$
This yields
$$-2 < \Delta t \cdot \left[ \frac{u}{\Delta x} \left( 1 - \cos{\kappa \Delta x} \right) + \frac{v}{\Delta y} \left( 1 - \cos{\mathcal{l} \Delta y} \right) \right] < 0$$
The right inequality does not limit $\Delta t$, but the left does:
$$\Delta t \cdot \left[ \frac{u}{\Delta x} \left( 1 - \cos{\kappa \Delta x} \right) + \frac{v}{\Delta y} \left( 1 - \cos{\mathcal{l} \Delta y} \right) \right] < 2$$
The worst-case scenario occurs when both cosines evaluate to -1, to hit this condition it is enough to have grid-scale noise, which has the maximum expressible wavenumber of $\kappa = \pi / \Delta x$. In this case we have
$$\Delta t \cdot \left( \frac{u}{\Delta x} + \frac{v}{\Delta y} \right)< 1$$
Since the direction is arbitrary you can substitute $u$ and $v$ with their absolute values and you get
$$\Delta t < \left( \frac{|u|}{\Delta x} + \frac{|v|}{\Delta y} \right)^{-1}$$
Same analysis in 3D yields the above CFL definition

[glwagner]

Why we would keep the function accurate_cell_advection_timescale as a duplicate of cell_advection_timescale? The word "accurate" is meaningless now.

[navidcy]

Agree!

[navidcy]

merge?