Background tracer distributions and advecting velocity fields

[glwagner]

We may able to provide an interface for specifying background tracer distributions and velocity fields (which perform advection only?) using FunctionFields and the existing advection operators.

I think we would want this functionality to assume that the linear balances between background terms are somehow separately satisfied. With that assumption we can include just the two nonlinear terms associated with 1) advection of the resolved field by the background field and 2) advection of the background field by the resolved fields.

For this to be easily implemented we need to change the function signature of the advection operators for momentum from

div_ũu(i, j, k, grid, advection, U)

to

div_ũu(i, j, k, grid, advection, U, u)

which then enables us to include background terms via

div_ũu(i, j, k, grid, advection, U_background, u) + div_ũu(i, j, k, grid, advection, U, u_background)

for example.

Tracers are fine, since the current function is

div_uc(i, j, k, grid, advection, U, c)

which will become

div_uc(i, j, k, grid, advection, U_background, c) + div_uc(i, j, k, grid, advection, U, c_background)

[glwagner]

This might require adding a new property to IncompressibleModel called background_fields (or something). Then users would write background_fields = (u=background_u, T=background_T)... We might also want to extend FunctionField to take parameters. @ali-ramadhan thoughts?

[glwagner]

This would hopefully make implementing the eady_turbulence.jl example a lot easier. Right now it requires quite a bit of expert knowledge:

Oceananigans.jl/examples/eady_turbulence.jl

Lines 92 to 132 in e808a82

	# The $x$-momentum forcing
	#
	# $ F_u = - α w - α (z + H) ∂ₓu $
	#
	# is applied at location `(f, c, c)`.
	function Fu_eady_func(i, j, k, grid, clock, state, p)
	return @inbounds (- p.α * ℑxzᶠᵃᶜ(i, j, k, grid, state.velocities.w)
	- p.α * (grid.zC[k] + p.H) * ∂xᶠᵃᵃ(i, j, k, grid, ℑxᶜᵃᵃ, state.velocities.u))
	end
	# The $y$-momentum forcing
	#
	# $ F_v = - α (z + H) ∂ₓv
	#
	# is applied at location `(c, f, c)`.
	function Fv_eady_func(i, j, k, grid, clock, state, p)
	return @inbounds -p.α * (grid.zC[k] + p.H) * ∂xᶜᵃᵃ(i, j, k, grid, ℑxᶠᵃᵃ, state.velocities.v)
	end
	# The $z$-momentum forcing
	#
	# $ F_w = - α (z + H) ∂ₓw
	#
	# is applied at location `(c, c, f)`.
	function Fw_eady_func(i, j, k, grid, clock, state, p)
	return @inbounds -p.α * (grid.zF[k] + p.H) * ∂xᶜᵃᵃ(i, j, k, grid, ℑxᶠᵃᵃ, state.velocities.w)
	end
	# The buoyancy forcing
	#
	# $ F_b = - α (z + H) ∂ₓb + α f v $
	#
	# is applied at location `(c, c, c)`.
	function Fb_eady_func(i, j, k, grid, clock, state, p)
	return @inbounds (- p.α * (grid.zC[k] + p.H) * ∂xᶜᵃᵃ(i, j, k, grid, ℑxᶠᵃᵃ, state.tracers.b)
	+ p.f * p.α * ℑyᵃᶜᵃ(i, j, k, grid, state.velocities.v))
	end

[ali-ramadhan]

This abstraction would also simplify the ECCO/SOSE LESbrary.jl setup although it's also messy due to other terms and me jamming in Interpolations.jl objects: https://github.com/CliMA/LESbrary.jl/blob/master/scripts/pilot_simulation.jl#L117-L145

It seems like some complex setups could benefit a lot from a background fields abstraction. It shouldn't affect performance since the two interaction terms should results in no-ops.

I think we would want this functionality to assume that the linear balances between background terms are somehow separately satisfied.

Agree with this. Users are responsible for making sure that they impose sensible background fields (just like they are responsible for setting sensible fields).