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:

```julia 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 ```

```julia 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