Use total tracer in diffusive flux calculation

[glwagner]

cc @hdrake @liuchihl

[glwagner]

Probably before merging this we also want to use the total velocities in the stress computation

[glwagner]

But I wonder also if we should add a feature to BackgroundFields that allows this to be toggled on and off

[glwagner]

Ok @hdrake @liuchihl I have added the option to include or not the background field when computing closure fluxes.

They are noit included by default. If you want to include them you need to build the BackgroundFields explicitly by writing something like

background_fields = BackgroundFields(; background_closure_fluxes=true, b=B)

where B is the background buoyancy field as before. Then pass this to the model constructor instead of passing a NamedTuple.

Let me know if this seems like a good interface and also if it works.

[hdrake]

Looks like a good interface to me.

But is it on purpose that there is only support for background fields in the NonhydrostaticModel and not for the HydrostaticFreeSurfaceModel?

@liuchihl will test it in our configurations.

[glwagner]

Looks like a good interface to me.

But is it on purpose that there is only support for background fields in the NonhydrostaticModel and not for the HydrostaticFreeSurfaceModel?

@liuchihl will test it in our configurations.

Well yes, it's substantial effort to support background fields. So we implemented it in the nonhydrostatic model first. Nobody has requested having background fields for the hydrostatic model. It's not impossible but might require some thinking if it's going to work with the more complicated turbulence closures (like CATKE or k-epsilon) that sometimes get used for hydrostatic applications.

Since the nonhydrostatic model is fast (at least on one GPU) the hydrostatic model is mostly important for simulations on the sphere (although this statement needs to be evaluated more carefully for complex domains when we have a proper nonhydrostatic solver).

[hdrake]

@glwagner, can you help us understand how the diffusive flux calculations here interact with gradient boundary conditions?

We've run some simulations and seem to be getting the expected behavior in the interior, where buoyancy tendencies are due to the convergence of both the perturbation and background diffusive fluxes.

However, we are not seeing the expected behavior from the flux convergence due to boundary conditions (on the domain, not even considering immersed boundary conditions yet). My expectation was that the GradientBoundaryCondition would be applied to whichever tracer fields are passed along to the tendency functions:

             - ∇_dot_qᶜ(i, j, k, grid, closure, diffusivities, val_tracer_index, closure_c, clock, model_fields, buoyancy)
             - immersed_∇_dot_qᶜ(i, j, k, grid, closure_c, c_immersed_bc, closure, diffusivities, val_tracer_index, clock, model_fields)

which in our case should be the sum of the perturbation and background tracer fields. Instead, it seems that our solutions are behaving as though the GradientBoundaryCondition is only being applied to the perturbation fluxes at the boundaries. Indeed, when we remove the background field from the gradient values passed to GradientBoundaryCondition, we get the behavior we are looking for.

[glwagner]

@glwagner, can you help us understand how the diffusive flux calculations here interact with gradient boundary conditions?

We've run some simulations and seem to be getting the expected behavior in the interior, where buoyancy tendencies are due to the convergence of both the perturbation and background diffusive fluxes.

However, we are not seeing the expected behavior from the flux convergence due to boundary conditions (on the domain, not even considering immersed boundary conditions yet). My expectation was that the GradientBoundaryCondition would be applied to whichever tracer fields are passed along to the tendency functions:

             - ∇_dot_qᶜ(i, j, k, grid, closure, diffusivities, val_tracer_index, closure_c, clock, model_fields, buoyancy)
             - immersed_∇_dot_qᶜ(i, j, k, grid, closure_c, c_immersed_bc, closure, diffusivities, val_tracer_index, clock, model_fields)

which in our case should be the sum of the perturbation and background tracer fields. Instead, it seems that our solutions are behaving as though the GradientBoundaryCondition is only being applied to the perturbation fluxes at the boundaries. Indeed, when we remove the background field from the gradient values passed to GradientBoundaryCondition, we get the behavior we are looking for.

GradientBoundaryCondition isn't applied to the fields at all. The tracer gradient across the boundary is the quantity we need in order to compute fluxes across boundaries. Therefore when the immersed boundary flux is computed, we call

Oceananigans.jl/src/ImmersedBoundaries/immersed_boundary_condition.jl

Lines 100 to 101 in 738d172

	@inline right_gradient(i, j, k, ibg, κ, Δ, bc::GBC, c, clock, fields) = getbc(bc, i, j, k, ibg, clock, fields)
	@inline left_gradient(i, j, k, ibg, κ, Δ, bc::GBC, c, clock, fields) = getbc(bc, i, j, k, ibg, clock, fields)

Contrast this with the routine for ValueBoundaryCondition which is a bit more involved:

Oceananigans.jl/src/ImmersedBoundaries/immersed_boundary_condition.jl

Lines 103 to 107 in 738d172

	@inline function right_gradient(i, j, k, ibg, κ, Δ, bc::VBC, c, clock, fields)
	cᵇ = getbc(bc, i, j, k, ibg, clock, fields)
	cⁱʲᵏ = @inbounds c[i, j, k]
	return 2 * (cᵇ - cⁱʲᵏ) / Δ
	end

These gradients are then used to compute the flux, for example

Oceananigans.jl/src/ImmersedBoundaries/immersed_boundary_condition.jl

Lines 117 to 122 in 738d172

	@inline function _west_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields)
	Δ = Δx(index_left(i, LX), j, k, ibg, LX, LY, LZ)
	κ = h_diffusivity(i, j, k, ibg, flip(LX), LY, LZ, closure, K, id, clock)
	∇c = left_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields)
	return - κ * ∇c
	end

Instead, it seems that our solutions are behaving as though the GradientBoundaryCondition is only being applied to the perturbation fluxes at the boundaries.

I don't know exactly what it means for a gradient to be applied to the field. Can you please clarify?

[glwagner]

Additionally can you help me understand why you are using GradientBoundaryCondition at all? I had thought that it's usually preferred to use FluxBoundaryCondition since this is often the preferred parameterization for large scale motions and basically we only use ValueBoundaryCondition for DNS (and GradientBoundaryCondition almost never). Curious to understand the class of problems that benefit from GradientBoundaryCondition.

[liuchihl]

@glwagner Thanks for implementing the total tracer diffusive flux at a high level. After running several tests, I found it to work exceptionally well! I conducted a series of tests: 1) comparing 1D vs 3D, 2) with and without the Coriolis force, and 3) with and without the immersed boundary. Everything looks great! Here are some simple examples on a rotated coordinate:

• 1D test with a small f:

nonconstantdiffusivity250days-theta.0.002_Nx4_Ny4_smallf_zlargerf.mp4

• 3D simulation with immersed grids:

nonconstantdiffusivity8days-theta.0.2_Nx4_Ny4_immersed_3Dfields_withcrossflux.mp4

The only caveat mentioned by @hdrake is that GradientBoundaryCondition is only being applied to the perturbation fluxes at the boundaries, i.e., GradientBoundaryCondition(-N^2*cos(θ)) is needed to make the total buoyancy gradient to be 0.

[hdrake]

Additionally can you help me understand why you are using GradientBoundaryCondition at all? [- @glwagner]

Flux boundary conditions are both conceptually and practically easier to implement when the boundary fluxes are zero or constant. They can be trickier when they depend on interior flow variables. In our case, for example, the boundary condition on the perturbation variable $b'$ (the buoyancy tracer in Oceananigans) is that the total diffusive flux should vanish, which means the perturbation diffusive flux needs to be minus the background diffusive flux.

Imposing a flux boundary condition requires knowing the diffusivity $\kappa$ right at the boundary. It is obvious how to implement this if the diffusivity is a constant, because the background diffusive flux is also a known constant, but less obvious how to do it when using a subgrid turbulence closure that yields a diffusivity that varies in space and time. @tomchor pointed out to use that we can sidestep this complexity if we just divide both sides of the boundary condition by $\kappa$, because then the boundary condition simply becomes that the buoyancy gradient across the boundary should just be equal to minus the background buoyancy gradient—a known constant in our problem.

[hdrake]

I don't know exactly what it means for a gradient to be applied to the field. Can you please clarify?

I just meant where in the code the gradient boundary conditions get imposed, which you've shown us is in the calculation of the gradients that feed into the downgradient diffusive fluxes that are used in the diffusive flux divergence contribution to the tracer tendencies. Thanks!

[glwagner]

I don't know exactly what it means for a gradient to be applied to the field. Can you please clarify?

I just meant where in the code the gradient boundary conditions get imposed, which you've shown us is in the calculation of the gradients that feed into the downgradient diffusive fluxes that are used in the diffusive flux divergence contribution to the tracer tendencies. Thanks!

Okay great. I would only add, I think it's clearer to think of the gradient as being used to diagnose the cross-boundary flux (rather than imposed).

I guess the point here is that there is actually an apparent flux of tracer into the perturbation field because of the presence of the background. So we are imagining that the background is being maintained by some large scale circulation which is ultimately the source of this apparent flux. FluxBoundaryCondition can be used if you know the diffusivity a priori but otherwise it does look like it will be simpler to use GradientBoundaryCondition to add this apparent flux contribution to the evolution of the perturbation.

[glwagner]

@glwagner Thanks for implementing the total tracer diffusive flux at a high level. After running several tests, I found it to work exceptionally well! I conducted a series of tests: 1) comparing 1D vs 3D, 2) with and without the Coriolis force, and 3) with and without the immersed boundary. Everything looks great! Here are some simple examples on a rotated coordinate:

• 1D test with a small f:

nonconstantdiffusivity250days-theta.0.002_Nx4_Ny4_smallf_zlargerf.mp4

• 3D simulation with immersed grids:

nonconstantdiffusivity8days-theta.0.2_Nx4_Ny4_immersed_3Dfields_withcrossflux.mp4
The only caveat mentioned by @hdrake is that GradientBoundaryCondition is only being applied to the perturbation fluxes at the boundaries, i.e., GradientBoundaryCondition(-N^2*cos(θ)) is needed to make the total buoyancy gradient to be 0.

Thanks @liuchihl !