PAR interpolation incorrect above top grid point

[jagoosw]

@syou83syou83 has noticed that interpolating the PAR field at a point above the top center point gives a value lower than that at the top point. This is probably because we've not been filling the halo region so its interpolating between [top value] and 0, the default for the halo region.

This should be easily fixed by filling the top halo points with the surface value, or surface ^ 2 / top point so that interpoalting to the surface will give the correct value. I will do this tomorrow.

[jagoosw]

This is slightly more complicated than I thought because we want to set the top point to $PAR_0(2 - A_0)$, and then have the normal fill_halo_regions! behaviour so I think I need to make the PAR fields its own type and then over load fill_halo_regions!

[jagoosw]

Probably easiest to give the PAR field a special boundary condition and then overload _fill_top_halo!?

[johnryantaylor]

Hi Both, Thanks for spotting this Si.  Jago: does it work to just set the PAR in the halo cells to the value of PAR at the first grid point?  I don’t think that It should matter very much whether we use the exponentially decaying form of the function because PAR is very unlikely to be a limiting factor at the surface (i.e. the growth function is very likely to be saturated in terms of PAR). John

…

On Feb 7, 2023 at 10:20 PM +0000, Jago Strong-Wright ***@***.***>, wrote: Probably easiest to give the PAR field a special boundary condition and then overload _fill_top_halo!? — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[glwagner]

Can you use ValueBoundaryCondition?

[jagoosw]

Can you use ValueBoundaryCondition?

Yeah I realised this this morning, thanks!

[jagoosw]

@glwagner have you got any idea what I'm doing wrong here? If I set the field to be:

field = CenterField(grid; boundary_conditions = regularize_field_boundary_conditions(FieldBoundaryConditions(top = ValueBoundaryCondition(surface_PAR)), grid, :PAR))

where surface_PAR is a function of x, y, t. When I call fill_halo_regions! I get this error:

nested task error: MethodError: objects of type Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Center, Center, Nothing, Oceananigans.BoundaryConditions.RightBoundary, typeof(surface_PAR), Nothing, Tuple{}, Tuple{}, Tuple{}} are not callable

ERROR: TaskFailedException
Stacktrace:
 [1] wait
   @ ./task.jl:345 [inlined]
 [2] wait
   @ ~/.julia/packages/KernelAbstractions/3ZHln/src/cpu.jl:65 [inlined]
 [3] wait(::KernelAbstractions.CPU, ev::KernelAbstractions.CPUEvent)
   @ KernelAbstractions ~/.julia/packages/KernelAbstractions/3ZHln/src/cpu.jl:64
 [4] fill_halo_event!(::Int64, ::Tuple{Vector{Function}, Vector{BoundaryCondition{C, Nothing} where C<:Oceananigans.BoundaryConditions.AbstractBoundaryConditionClassification}, Vector{BoundaryCondition}}, ::OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, ::Tuple{Colon, Colon, Colon}, ::Tuple{Center, Center, Center}, ::CPU, ::KernelAbstractions.NoneEvent, ::RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}; kwargs::Base.Pairs{Symbol, Tuple{}, Tuple{Symbol}, NamedTuple{(:reduced_dimensions,), Tuple{Tuple{}}}})
   @ Oceananigans.BoundaryConditions ~/.julia/packages/Oceananigans/dXuua/src/BoundaryConditions/fill_halo_regions.jl:59
 [5] #fill_halo_regions!#18
   @ ~/.julia/packages/Oceananigans/dXuua/src/BoundaryConditions/fill_halo_regions.jl:43 [inlined]
 [6] #fill_halo_regions!#74
   @ ~/.julia/packages/Oceananigans/dXuua/src/Fields/field.jl:705 [inlined]
 [7] fill_halo_regions!(::Field{Center, Center, Center, Nothing, RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, Float64, FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Value, Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Center, Center, Nothing, Oceananigans.BoundaryConditions.RightBoundary, typeof(surface_PAR), Nothing, Tuple{}, Tuple{}, Tuple{}}}, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing}})
   @ Oceananigans.Fields ~/.julia/packages/Oceananigans/dXuua/src/Fields/field.jl:693
 [8] top-level scope
   @ REPL[141]:1
 [9] top-level scope
   @ ~/.julia/packages/CUDA/Ey3w2/src/initialization.jl:52

    nested task error: MethodError: objects of type Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Center, Center, Nothing, Oceananigans.BoundaryConditions.RightBoundary, typeof(surface_PAR), Nothing, Tuple{}, Tuple{}, Tuple{}} are not callable
    Stacktrace:
      [1] call
        @ ~/.julia/packages/Cassette/34vIw/src/context.jl:456 [inlined]
      [2] fallback
        @ ~/.julia/packages/Cassette/34vIw/src/context.jl:454 [inlined]
      [3] _overdub_fallback(::Any, ::Vararg{Any})
        @ ~/.julia/packages/Cassette/34vIw/src/overdub.jl:586 [inlined]
      [4] overdub
        @ ~/.julia/packages/Cassette/34vIw/src/overdub.jl:586 [inlined]
      [5] getbc(::BoundaryCondition{Oceananigans.BoundaryConditions.Value, Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Center, Center, Nothing, Oceananigans.BoundaryConditions.RightBoundary, typeof(surface_PAR), Nothing, Tuple{}, Tuple{}, Tuple{}}}, ::Int64, ::Int64, ::RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU})
        @ ~/.julia/packages/Oceananigans/dXuua/src/BoundaryConditions/boundary_condition.jl:106 [inlined]
      [6] overdub
        @ ~/.julia/packages/Oceananigans/dXuua/src/BoundaryConditions/boundary_condition.jl:106 [inlined]
      [7] right_gradient(::BoundaryCondition{Oceananigans.BoundaryConditions.Value, Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Center, Center, Nothing, Oceananigans.BoundaryConditions.RightBoundary, typeof(surface_PAR), Nothing, Tuple{}, Tuple{}, Tuple{}}}, ::Float64, ::Float64, ::Int64, ::Int64, ::RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU})
        @ ~/.julia/packages/Oceananigans/dXuua/src/BoundaryConditions/fill_halo_regions_value_gradient.jl:13 [inlined]
      [8] overdub
        @ ~/.julia/packages/Oceananigans/dXuua/src/BoundaryConditions/fill_halo_regions_value_gradient.jl:13 [inlined]
      [9] overdub
        @ ~/.julia/packages/Oceananigans/dXuua/src/BoundaryConditions/fill_halo_regions_value_gradient.jl:97 [inlined]
     [10] macro expansion
        @ ~/.julia/packages/Oceananigans/dXuua/src/BoundaryConditions/fill_halo_regions.jl:135 [inlined]
     [11] overdub
        @ ~/.julia/packages/KernelAbstractions/3ZHln/src/macros.jl:266 [inlined]
     [12] __thread_run(tid::Int64, len::Int64, rem::Int64, obj::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_bottom_and_top_halo!)}, ndrange::Nothing, iterspace::KernelAbstractions.NDIteration.NDRange{2, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Nothing, Nothing}, args::Tuple{OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Value, Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Center, Center, Nothing, Oceananigans.BoundaryConditions.RightBoundary, typeof(surface_PAR), Nothing, Tuple{}, Tuple{}, Tuple{}}}, Tuple{Int64, Int64}, Tuple{Center, Center, Center}, RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}}, dynamic::KernelAbstractions.NDIteration.NoDynamicCheck)
        @ KernelAbstractions ~/.julia/packages/KernelAbstractions/3ZHln/src/cpu.jl:157
     [13] __run(obj::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_bottom_and_top_halo!)}, ndrange::Nothing, iterspace::KernelAbstractions.NDIteration.NDRange{2, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Nothing, Nothing}, args::Tuple{OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Value, Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Center, Center, Nothing, Oceananigans.BoundaryConditions.RightBoundary, typeof(surface_PAR), Nothing, Tuple{}, Tuple{}, Tuple{}}}, Tuple{Int64, Int64}, Tuple{Center, Center, Center}, RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}}, dynamic::KernelAbstractions.NDIteration.NoDynamicCheck)
        @ KernelAbstractions ~/.julia/packages/KernelAbstractions/3ZHln/src/cpu.jl:130
     [14] (::KernelAbstractions.var"#37#38"{Tuple{KernelAbstractions.NoneEvent}, Nothing, typeof(KernelAbstractions.__run), Tuple{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_bottom_and_top_halo!)}, Nothing, KernelAbstractions.NDIteration.NDRange{2, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Nothing, Nothing}, Tuple{OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Value, Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Center, Center, Nothing, Oceananigans.BoundaryConditions.RightBoundary, typeof(surface_PAR), Nothing, Tuple{}, Tuple{}, Tuple{}}}, Tuple{Int64, Int64}, Tuple{Center, Center, Center}, RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}}, KernelAbstractions.NDIteration.NoDynamicCheck}})()
        @ KernelAbstractions ~/.julia/packages/KernelAbstractions/3ZHln/src/cpu.jl:22

I've tried a few different things but can't get this to work with a function.

[jagoosw]

Never mind, worked this out now!

[iuryt]

Never mind, worked this out now!

What did you have to change to make it work?
I thought it was because you were using this function regularize_field_boundary_conditions

[jagoosw]

Never mind, worked this out now!

What did you have to change to make it work? I thought it was because you were using this function regularize_field_boundary_conditions

If I don't use that it fails to build the field because it tries to set a default boundary condition in the periodic directions. The issue was that I wasn't using fill_halo_regions! correctly (was doing fill_halo_regions!(field) not fill_halo_regions!(field, model.clock, fields(model)).

[glwagner]

It might be more robust to use a discrete form boundary condition here so you don't have to regularize

[glwagner]

The other issue is that regularize_field_boundary_conditions is not supposed to be user facing so it could change...

More broadly I think this use case motivates rethinking how we initialize boundary conditions in Oceananigans to make user code more robust / easier to write