Help with outputs of horizontal derivatives of vertically integrated fields

[scott-conn]

Hi everyone,

I'm trying to output horizontal derivatives of vertically integrated fields. For example, something of the form:
$$\frac{\partial}{\partial x}\int_{-L_z}^0 u dz.$$

However, I'm having trouble when I add bathymetry to the simulation. The output is zero wherever the depth of the bathymetry is different to that of the underlying grid. I don't have this problem if I just output the vertically integrated field itself.

Here's an example where I tried adding a gaussian ridge to the baroclinic adjustment example from the docs.

using Oceananigans
using Oceananigans.Units
using Oceananigans.ImmersedBoundaries: GridFittedBottom

Lx = 1000kilometers
Ly = 1000kilometers
Lz = 1kilometers

underlying_grid = RectilinearGrid(size = (48, 48, 8),
                                x = (0, Lx),
                                y = (-Ly/2, Ly/2),
                                z = (-Lz, 0),
                                topology = (Periodic, Bounded, Bounded),
                                halo = (4,4,4))    

#add bathymetry: gaussian ridge
H = 300meters
δ = 200kilometers
h(x,y) = -Lz + H * exp(- (x - Lx/2)^2 / δ^2)
grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(h))

model = HydrostaticFreeSurfaceModel(; grid,
                                    coriolis = BetaPlane(latitude = -45),
                                    buoyancy = BuoyancyTracer(),
                                    tracers = :b,
                                    momentum_advection = WENO(),
                                    tracer_advection = WENO())


ramp(y, Δy) = min(max(0, y/Δy + 1/2), 1)

N² = 1e-5 # [s⁻²] buoyancy frequency / stratification
M² = 1e-7 # [s⁻²] horizontal buoyancy gradient

Δy = 100kilometers # width of the region of the front
Δb = Δy * M²       # buoyancy jump associated with the front
ϵb = 1e-2 * Δb     # noise amplitude

bᵢ(x, y, z) = N² * z + Δb * ramp(y, Δy) + ϵb * randn()

set!(model, b=bᵢ)

simulation = Simulation(model, Δt=20minutes, stop_time=20days)
conjure_time_step_wizard!(simulation, IterationInterval(20), cfl=0.2, max_Δt=20minutes)

#add output for ∂U/∂x
U = Field(Integral(model.velocities.u,dims=(3)))
∂x_U = ∂x(U)
filename = "baroclinic_adjustment"
save_fields_interval = 0.5day
simulation.output_writers[:U] = JLD2OutputWriter(model, (; ∂x_U);
                                                       filename = filename,
                                                       schedule = TimeInterval(save_fields_interval),
                                                       overwrite_existing = true)

@info "Running the simulation..."

run!(simulation)

@info "Simulation completed in " * prettytime(simulation.run_wall_time)

A snapshot of the output looks like:

You can see that the output is masked to zero in the centre region where the ridge is located. Any suggestions?

On a related note, I'd also like to output terms like
$$\frac{\partial}{\partial y}\left[f\int_{-L_z}^0 u dz\right].$$
I couldn't figure out how to reconstruct $f$ in the outputs. Is there a way to do this?

Thanks!

[glwagner]

Nice example! It looks kinda like a bug. I think @simone-silvestri might be able to help

[glwagner]

The code that tries to support reductions on immersed boundary grids is here:

Oceananigans.jl/src/ImmersedBoundaries/immersed_reductions.jl

Lines 52 to 86 in 7a4b3f0

	# ImmersedReducedFields
	const XIRF = AbstractField{Nothing, <:Any, <:Any, <:ImmersedBoundaryGrid}
	const YIRF = AbstractField{<:Any, Nothing, <:Any, <:ImmersedBoundaryGrid}
	const ZIRF = AbstractField{<:Any, <:Any, Nothing, <:ImmersedBoundaryGrid}
	const YZIRF = AbstractField{<:Any, Nothing, Nothing, <:ImmersedBoundaryGrid}
	const XZIRF = AbstractField{Nothing, <:Any, Nothing, <:ImmersedBoundaryGrid}
	const XYIRF = AbstractField{Nothing, Nothing, <:Any, <:ImmersedBoundaryGrid}
	const XYZIRF = AbstractField{Nothing, Nothing, Nothing, <:ImmersedBoundaryGrid}
	const IRF = Union{XIRF, YIRF, ZIRF, YZIRF, XZIRF, XYIRF, XYZIRF}
	@inline condition_operand(func::Function, op::IRF, cond, mask) = ConditionalOperation(op; func, condition=NotImmersedColumn(immersed_column(op), cond ), mask)
	@inline condition_operand(func::Function, op::IRF, ::Nothing, mask) = ConditionalOperation(op; func, condition=NotImmersedColumn(immersed_column(op), truefunc), mask)
	@inline condition_operand(func::typeof(identity), op::IRF, ::Nothing, mask) = ConditionalOperation(op; func, condition=NotImmersedColumn(immersed_column(op), truefunc), mask)
	@inline function immersed_column(field::IRF)
	grid = field.grid
	reduced_dims = reduced_dimensions(field)
	LX, LY, LZ = map(center_to_nothing, location(field))
	one_field = ConditionalOperation{LX, LY, LZ}(OneField(Int), identity, grid, NotImmersed(truefunc), zero(grid))
	return sum(one_field, dims=reduced_dims)
	end
	@inline center_to_nothing(::Type{Face}) = Face
	@inline center_to_nothing(::Type{Center}) = Center
	@inline center_to_nothing(::Type{Nothing}) = Center
	@inline function evaluate_condition(condition::NotImmersedColumn, i, j, k, ibg, co::ConditionalOperation, args...)
	LX, LY, LZ = location(co)
	return evaluate_condition(condition.func, i, j, k, ibg, args...) & !(is_immersed_column(i, j, k, condition.immersed_column))
	end
	@inline is_immersed_column(i, j, k, column) = @inbounds column[i, j, k] == 0

I think specifically this conditional:

Oceananigans.jl/src/ImmersedBoundaries/immersed_reductions.jl

Line 86 in 7a4b3f0

@inline is_immersed_column(i, j, k, column) = @inbounds column[i, j, k] == 0

may be evaluated wrong.

Possibly could test by manually setting that to true to see if it fixes your bug @scott-conn eg

is_immersed_column(i, j, k, column) = true

[glwagner]

I couldn't figure out how to reconstruct in the outputs. Is there a way to do this?

There isn't right now --- you have to manually compute f = f0 + beta * y. You could raise an issue with a feature request for a Field that corresponds to the Coriolis parameter. I think that would be pretty easy to design.

[simone-silvestri]

Interesting, it looks like the U field is computed correctly, but the derivative excludes immersed regions on 1. I ll do some digging

[simone-silvestri]

Oh, I see why this happens. The derivative internally calls ∂xᶜᶜᶜ which evaluates to zero in an immersed boundary. Good bug catching we can fix it

[glwagner]

Seems like we need to use k=nothing or something rather than k=1 so it can be ignored by conditionals.

[glwagner]

Or actually, we may need a special method for ∂xᶜᶜᶜ that dispatches on reduced fields.

This will work because (unlike non-reduced fields), ReducedField retain their location information on the GPU.

[simone-silvestri]

That would require a slight change because ∂xᶜᶜᶜ is a local method while operations on reduced fields will probably require information of the whole column. We can make ∂xᶜᶜᶜ work by looping over the reduced dimensions though.

[glwagner]

I'm not sure I grasp what you're saying. Are you saying that ReducedField needs compute the reduced immersed grid that effectively applies to it? What is the problem with doing that in the constructor for a reduced field?

[glwagner]

@scott-conn as a stop gap, I think you can try putting this somewhere in your script:

import Oceananigans.Operators: ∂xᶜᶜᶜ
using Oceananigans.ImmersedBoundaries: XYIRF, ImmersedBoundaryGrid
@inline ∂xᶜᶜᶜ(i, j, k, ibg::ImmersedBoundaryGrid, U::XYIRF) = ∂xᶜᶜᶜ(i, j, k, ibg.underlying_grid, U)

This will also compute the derivative within the inactive regions. So you'll have to ignore those in output. That said we don't have a great system for masking inactive regions yet anyways, heh. Curious if this works for you.

[simone-silvestri]

can you post the error message you get?

[scott-conn]

Ah sorry if I wasn't clear, there isn't an error message. The code runs fine but the output still defaults to zero in the immersed region.

[glwagner]

I think I know what's wrong. You can try this

import Oceananigans.Operators: ∂xᶜᶜᶜ
using Oceananigans.Fields: ZReducedField
@inline ∂xᶜᶜᶜ(i, j, k, ibg::ImmersedBoundaryGrid, U::ZReducedField) = ∂xᶜᶜᶜ(i, j, k, ibg.underlying_grid, U)

it's no less correct either. (On the GPU, we toss the grid from the Field so it's just a "regular" reduced field, not an immersed reduced field.)

[scott-conn]

That did the trick, thanks!

[glwagner]

Excellent. That should also work on the CPU.