Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix _immered_cell and adapt for PartialCellBottom #3682

Merged
merged 25 commits into from
Sep 27, 2024
Merged

Conversation

simone-silvestri
Copy link
Collaborator

closes #3681

We can probably use this PR to add also a validation for the PartialCellBottom and the bug-fixes in the implementation.
What do people think?

cc @jm-c

@glwagner
Copy link
Member

glwagner commented Aug 5, 2024

@ali-ramadhan

@glwagner
Copy link
Member

glwagner commented Aug 5, 2024

closes #3681

We can probably use this PR to add also a validation for the PartialCellBottom and the bug-fixes in the implementation. What do people think?

cc @jm-c

presumably we should always use PartialCellBottom in examples with bathymetry. Is there a reason to ever use GridFittedBottom?

@navidcy navidcy added bug 🐞 Even a perfect program still has bugs GPU 👾 Where Oceananigans gets its powers from labels Aug 9, 2024
@glwagner
Copy link
Member

How does this PR change adapt?

Comment on lines 98 to 104
# Face node below current cell
z = znode(i, j, k, underlying_grid, c, c, f)
h = @inbounds ib.bottom_height[i, j, 1]
ϵ = ib.minimum_fractional_cell_height
# z + Δz is equal to the face above the current cell
Δz = Δzᶜᶜᶜ(i, j, k, ibg.underlying_grid)
return z + Δz * ϵ ≤ h
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we expect this change to fix the weird behaviors we were seeing in #3529 ?

cc @francispoulin

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not sure, @jm-c is the author of this fix, and he said that with this fix the immersed boundary example did not NaN immediately but it still NaNed after a little bit so we believe there is more bugs to be squashed. You can try that example on this and see if it makes any difference. It will most likely change the solution.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am happy to try the example from #3529 it on this branch, when you think it's ready for testing.

@glwagner
Copy link
Member

What are we waiting for to merge this PR?

@simone-silvestri
Copy link
Collaborator Author

I guess to have the partial cells working properly, @jm-c was saying there is still some fixes to be done.

@glwagner
Copy link
Member

The first thing maybe is to add a test that catches the bug.

@simone-silvestri
Copy link
Collaborator Author

@ali-ramadhan does this PR fix your issue? I will try the internal tide example on this PR, add some tests for PartialCellBottom, and then we can probably merge it. I am not sure whether this solves all the issues with PartialCellBottom, probably more validation is required, but this is a first step.

@simone-silvestri
Copy link
Collaborator Author

This is the status of the internal_tide example on this PR, not producing NaNs anymore, but still not really working.

internal_tide.mp4

@ali-ramadhan
Copy link
Member

ali-ramadhan commented Aug 29, 2024

@simone-silvestri I think another error has now come up with the original MWE in issue #3681. I tried the main branch and this branch on two different systems and I'm getting a different CUDA illegal memory access error. Since it's probably not related to this PR I'll open a new issue about it. I'm guessing something is out of bounds in an @inbounds. I'll try catching it with --check-bounds=yes.

MWE:

using Oceananigans
using Oceananigans.ImmersedBoundaries: PartialCellBottom

underlying_grid = RectilinearGrid(GPU(), topology=(Periodic, Flat, Bounded), size=(4, 4), x=(0, 1), z=(0, 1))
slope(x) = x
grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(slope))

model = HydrostaticFreeSurfaceModel(; grid)

Error:

ERROR: CUDA error: an illegal memory access was encountered (code 700, ERROR_ILLEGAL_ADDRESS)
Stacktrace:
  [1] throw_api_error(res::CUDA.cudaError_enum)
    @ CUDA ~/.julia/packages/CUDA/Tl08O/lib/cudadrv/libcuda.jl:30
  [2] check
    @ ~/.julia/packages/CUDA/Tl08O/lib/cudadrv/libcuda.jl:37 [inlined]
  [3] cuStreamGetCaptureInfo
    @ ~/.julia/packages/CUDA/Tl08O/lib/utils/call.jl:34 [inlined]
  [4] capture_status(stream::CUDA.CuStream)
    @ CUDA ~/.julia/packages/CUDA/Tl08O/lib/cudadrv/graph.jl:174
  [5] is_capturing
    @ ~/.julia/packages/CUDA/Tl08O/lib/cudadrv/graph.jl:179 [inlined]
  [6] convert(::Type{CUDA.CuPtr{Float64}}, managed::CUDA.Managed{CUDA.DeviceMemory})
    @ CUDA ~/.julia/packages/CUDA/Tl08O/src/memory.jl:539
  [7] unsafe_convert
    @ ~/.julia/packages/CUDA/Tl08O/src/array.jl:434 [inlined]
  [8] #pointer#1123
    @ ~/.julia/packages/CUDA/Tl08O/src/array.jl:392 [inlined]
  [9] pointer (repeats 2 times)
    @ ~/.julia/packages/CUDA/Tl08O/src/array.jl:384 [inlined]
 [10] unsafe_convert(::Type{CUDA.CuDeviceArray{Float64, 3, 1}}, a::CUDA.CuArray{Float64, 3, CUDA.DeviceMemory})
    @ CUDA ~/.julia/packages/CUDA/Tl08O/src/array.jl:454
 [11] adapt_storage
    @ ~/.julia/packages/CUDA/Tl08O/src/compiler/execution.jl:162 [inlined]
 [12] adapt_structure
    @ ~/.julia/packages/Adapt/7T9au/src/Adapt.jl:57 [inlined]
 [13] adapt
    @ ~/.julia/packages/Adapt/7T9au/src/Adapt.jl:40 [inlined]
 [14] #1
    @ ~/.julia/packages/OffsetArrays/hwmnB/ext/OffsetArraysAdaptExt.jl:9 [inlined]
 [15] parent_call
    @ ~/.julia/packages/OffsetArrays/hwmnB/src/OffsetArrays.jl:315 [inlined]
 [16] adapt_structure
    @ ~/.julia/packages/OffsetArrays/hwmnB/ext/OffsetArraysAdaptExt.jl:9 [inlined]
 [17] adapt
    @ ~/.julia/packages/Adapt/7T9au/src/Adapt.jl:40 [inlined]
 [18] cudaconvert
    @ ~/.julia/packages/CUDA/Tl08O/src/compiler/execution.jl:198 [inlined]
 [19] map
    @ ./tuple.jl:294 [inlined]
 [20] map(f::typeof(CUDA.cudaconvert), t::Tuple{KernelAbstractions.CompilerMetadata{…}, OffsetArrays.OffsetArray{…}, Nothing, Nothing, Tuple{…}, ImmersedBoundaryGrid{…}, Tuple{}})
    @ Base ./tuple.jl:294
 [21] macro expansion
    @ ~/.julia/packages/CUDA/Tl08O/src/compiler/execution.jl:110 [inlined]
 [22] (::KernelAbstractions.Kernel{…})(::OffsetArrays.OffsetArray{…}, ::Vararg{…}; ndrange::Nothing, workgroupsize::Nothing)
    @ CUDA.CUDAKernels ~/.julia/packages/CUDA/Tl08O/src/CUDAKernels.jl:103
 [23] (::KernelAbstractions.Kernel{CUDA.CUDAKernels.CUDABackend, KernelAbstractions.NDIteration.StaticSize{…}, KernelAbstractions.NDIteration.StaticSize{…}, typeof(Oceananigans.BoundaryConditions.gpu__no_fill!)})(::OffsetArrays.OffsetArray{Float64, 3, CUDA.CuArray{…}}, ::Vararg{Any})
    @ CUDA.CUDAKernels ~/.julia/packages/CUDA/Tl08O/src/CUDAKernels.jl:89
 [24] launch!(::GPU, ::ImmersedBoundaryGrid{…}, ::Tuple{…}, ::typeof(Oceananigans.BoundaryConditions._no_fill!), ::OffsetArrays.OffsetArray{…}, ::Vararg{…}; include_right_boundaries::Bool, reduced_dimensions::Tuple{}, location::Nothing, active_cells_map::Nothing, kwargs::@Kwargs{})
    @ Oceananigans.Utils /global/u1/a/alir/Oceananigans.jl/src/Utils/kernel_launching.jl:168
 [25] launch!
    @ /global/u1/a/alir/Oceananigans.jl/src/Utils/kernel_launching.jl:154 [inlined]
 [26] #fill_open_boundary_regions!#48
    @ /global/u1/a/alir/Oceananigans.jl/src/BoundaryConditions/fill_halo_regions_open.jl:20 [inlined]
 [27] fill_open_boundary_regions!
    @ /global/u1/a/alir/Oceananigans.jl/src/BoundaryConditions/fill_halo_regions_open.jl:9 [inlined]
 [28] fill_halo_regions!(::OffsetArrays.OffsetArray{…}, ::FieldBoundaryConditions{…}, ::Tuple{…}, ::Tuple{…}, ::ImmersedBoundaryGrid{…}; fill_boundary_normal_velocities::Bool, kwargs::@Kwargs{})
    @ Oceananigans.BoundaryConditions /global/u1/a/alir/Oceananigans.jl/src/BoundaryConditions/fill_halo_regions.jl:53
 [29] fill_halo_regions!
    @ /global/u1/a/alir/Oceananigans.jl/src/BoundaryConditions/fill_halo_regions.jl:48 [inlined]
 [30] #fill_halo_regions!#61
    @ /global/u1/a/alir/Oceananigans.jl/src/Fields/field.jl:773 [inlined]
 [31] fill_halo_regions!(::Field{Face, Center, Nothing, Nothing, ImmersedBoundaryGrid{…}, Tuple{…}, OffsetArrays.OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{…}})
    @ Oceananigans.Fields /global/u1/a/alir/Oceananigans.jl/src/Fields/field.jl:770
 [32] fill_halo_regions!(::Tuple{Field{…}, Field{…}}; kwargs::@Kwargs{})
    @ Oceananigans.Fields /global/u1/a/alir/Oceananigans.jl/src/Fields/field_tuples.jl:64
 [33] fill_halo_regions!(::Tuple{Field{…}, Field{…}})
    @ Oceananigans.Fields /global/u1/a/alir/Oceananigans.jl/src/Fields/field_tuples.jl:56
 [34] Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitAuxiliaryFields(grid::ImmersedBoundaryGrid{Float64, Periodic, Flat, Bounded, RectilinearGrid{…}, PartialCellBottom{…}, Nothing, Nothing, GPU})
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels /global/u1/a/alir/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl:239
 [35] materialize_free_surface(free_surface::SplitExplicitFreeSurface{…}, velocities::@NamedTuple{}, grid::ImmersedBoundaryGrid{…})
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels /global/u1/a/alir/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl:121
 [36] HydrostaticFreeSurfaceModel(; grid::ImmersedBoundaryGrid{…}, clock::Clock{…}, momentum_advection::Centered{…}, tracer_advection::Centered{…}, buoyancy::Nothing, coriolis::Nothing, free_surface::SplitExplicitFreeSurface{…}, tracers::Nothing, forcing::@NamedTuple{}, closure::Nothing, boundary_conditions::@NamedTuple{}, particles::Nothing, biogeochemistry::Nothing, velocities::Nothing, pressure::Nothing, diffusivity_fields::Nothing, auxiliary_fields::@NamedTuple{})
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels /global/u1/a/alir/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl:177
 [37] top-level scope
    @ REPL[6]:1
Some type information was truncated. Use `show(err)` to see complete types.

@simone-silvestri
Copy link
Collaborator Author

Hmmm, I cannot reproduce it on the CPU with --check-bounds=yes. Perhaps it was fixed in the latest three commits?

@inline Δzᶠᶜᶜ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶜᶜ(i, j, k, ibg)
@inline Δzᶠᶜᶠ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶜᶠ(i, j, k, ibg)
@inline Δzᶜᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶜᶜᶜ(i, j, k, ibg)
@inline Δzᶜᶠᶠ(i, j, k, ibg::YFlatPCBIBG) = Δzᶜᶜᶠ(i, j, k, ibg)
@inline Δzᶠᶠᶜ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶠᶜ(i, j, k, ibg)
@inline Δzᶠᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶠᶜᶜ(i, j, k, ibg)

@ali-ramadhan
Copy link
Member

ali-ramadhan commented Aug 29, 2024

I was able to reproduce on the CPU with --check-bounds=yes but on the main branch. Lemme try on this branch.


Stacktrace from main for reference:

julia> using Oceananigans

julia> using Oceananigans.ImmersedBoundaries: PartialCellBottom

julia> underlying_grid = RectilinearGrid(CPU(), topology=(Periodic, Flat, Bounded), size=(4, 4), x=(0, 1), z=(0, 1))
4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── Periodic x ∈ [0.0, 1.0) regularly spaced with Δx=0.25
├── Flat y                  
└── Bounded  z ∈ [0.0, 1.0] regularly spaced with Δz=0.25

julia> slope(x) = x
slope (generic function with 1 method)

julia> grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(slope))
4×1×4 ImmersedBoundaryGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo:
├── immersed_boundary: PartialCellBottom(mean(z)=0.5, min(z)=0.125, max(z)=0.875, ϵ=0.2)
├── underlying_grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── Periodic x ∈ [0.0, 1.0) regularly spaced with Δx=0.25
├── Flat y                  
└── Bounded  z ∈ [0.0, 1.0] regularly spaced with Δz=0.25

julia> model = HydrostaticFreeSurfaceModel(; grid)
ERROR: BoundsError: attempt to access 10×1×1 OffsetArray(::Array{Float64, 3}, -2:7, 1:1, 1:1) with eltype Float64 with indices -2:7×1:1×1:1 at index [1, 0, 1]
Stacktrace:
  [1] throw_boundserror(A::OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, I::Tuple{Int64, Int64, Int64})
    @ Base ./abstractarray.jl:737
  [2] checkbounds
    @ ./abstractarray.jl:702 [inlined]
  [3] getindex
    @ ~/.julia/packages/OffsetArrays/hwmnB/src/OffsetArrays.jl:422 [inlined]
  [4] getindex
    @ ~/atdepth/Oceananigans.jl/src/Fields/field.jl:541 [inlined]
  [5] Δzᶜᶜᶜ
    @ ~/atdepth/Oceananigans.jl/src/ImmersedBoundaries/partial_cell_bottom.jl:115 [inlined]
  [6] Δzᶜᶠᶜ
    @ ~/atdepth/Oceananigans.jl/src/ImmersedBoundaries/partial_cell_bottom.jl:141 [inlined]
  [7] getindex
    @ ~/atdepth/Oceananigans.jl/src/AbstractOperations/grid_metrics.jl:138 [inlined]
  [8] getindex
    @ ~/atdepth/Oceananigans.jl/src/AbstractOperations/conditional_operations.jl:101 [inlined]
  [9] _getindex
    @ ./abstractarray.jl:1341 [inlined]
 [10] getindex
    @ ./abstractarray.jl:1291 [inlined]
 [11] macro expansion
    @ ./reducedim.jl:317 [inlined]
 [12] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [13] _mapreducedim!(f::typeof(identity), op::typeof(Base.add_sum), R::SubArray{…}, A::Oceananigans.AbstractOperations.ConditionalOperation{…})
    @ Base ./reducedim.jl:316
 [14] mapreducedim!(f::Function, op::Function, R::SubArray{…}, A::Oceananigans.AbstractOperations.ConditionalOperation{…})
    @ Base ./reducedim.jl:324
 [15] sum!(f::Function, r::SubArray{…}, A::Oceananigans.AbstractOperations.ConditionalOperation{…})
    @ Base ./reducedim.jl:1034
 [16] sum!(r::Field{…}, a::Oceananigans.AbstractOperations.GridMetricOperation{…}; condition::Nothing, mask::Int64, kwargs::@Kwargs{})
    @ Oceananigans.Fields ~/atdepth/Oceananigans.jl/src/Fields/field.jl:696
 [17] sum!(r::Field{…}, a::Oceananigans.AbstractOperations.GridMetricOperation{…})
    @ Oceananigans.Fields ~/atdepth/Oceananigans.jl/src/Fields/field.jl:690
 [18] Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitAuxiliaryFields(grid::ImmersedBoundaryGrid{Float64, Periodic, Flat, Bounded, RectilinearGrid{…}, PartialCellBottom{…}, Nothing, Nothing, CPU})
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/atdepth/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl:237
 [19] materialize_free_surface(free_surface::SplitExplicitFreeSurface{…}, velocities::@NamedTuple{…}, grid::ImmersedBoundaryGrid{…})
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/atdepth/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl:121
 [20] HydrostaticFreeSurfaceModel(; grid::ImmersedBoundaryGrid{…}, clock::Clock{…}, momentum_advection::Centered{…}, tracer_advection::Centered{…}, buoyancy::Nothing, coriolis::Nothing, free_surface::SplitExplicitFreeSurface{…}, tracers::Nothing, forcing::@NamedTuple{}, closure::Nothing, boundary_conditions::@NamedTuple{}, particles::Nothing, biogeochemistry::Nothing, velocities::Nothing, pressure::Nothing, diffusivity_fields::Nothing, auxiliary_fields::@NamedTuple{})
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/atdepth/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl:177
 [21] top-level scope
    @ REPL[6]:1
Some type information was truncated. Use `show(err)` to see complete types.

@ali-ramadhan
Copy link
Member

Confirming that the MWE just above works on this branch on both CPU and GPU!

@ali-ramadhan
Copy link
Member

Also confirming that my original MWE from #3681 also works on this branch now!

using Oceananigans
using Oceananigans.ImmersedBoundaries: PartialCellBottom

underlying_grid = RectilinearGrid(GPU(), topology=(Periodic, Flat, Bounded), size=(4, 4), x=(0, 1), z=(0, 1))
slope(x) = x
grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(slope))

model = HydrostaticFreeSurfaceModel(; grid)
time_step!(model, 1)

@simone-silvestri simone-silvestri added the 🚨 DO NOT MERGE 🚨 IN BIG BOLD RED CAPS WITH FLASHING SIRENS label Sep 4, 2024
@glwagner
Copy link
Member

What's keeping us from merging this? It fixes something so its a step forward, I think we should merge it.

@glwagner
Copy link
Member

I implemented the clamping from #3787 and it seems to have fixed the issue:

GridFittedBottom

internal_tide_gf.mp4

PartialCellBottom

internal_tide_pc.mp4

I'm not sure I really understand why it would make a difference. I'm wondering @simone-silvestri are you sure that your test was using this code?

Either way I see no reason not to merge this. Perhaps we will continue to find improvements for PartialCellBottom but I think this is ready to be used in an example. Note that in the internal tide case the partial cells have the interesting effect of weakening the vertical velocity. That could make sense if the topography is somehow "less effectively steep" than in the full cell case.

@glwagner
Copy link
Member

This PR now supercedes #3787

@simone-silvestri simone-silvestri removed the 🚨 DO NOT MERGE 🚨 IN BIG BOLD RED CAPS WITH FLASHING SIRENS label Sep 25, 2024
@glwagner glwagner merged commit 675e0ba into main Sep 27, 2024
46 checks passed
@francispoulin
Copy link
Collaborator

Now we have merged, should I revisit #3529?

@glwagner
Copy link
Member

It's a good idea! Also @apaloczy had some tests that might be good to try to re-run.

@glwagner glwagner changed the title Fix a bug in the adapt of PartialCellBottom Fix _immered_cell and adapt for PartialCellBottom Sep 29, 2024
@glwagner
Copy link
Member

I suspect it won't change anything though because this PR fixes a bug that affects when criteria about the minimum fractional cell height is met. I don't think that criteria is met in your example (which is good)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug 🐞 Even a perfect program still has bugs GPU 👾 Where Oceananigans gets its powers from
Projects
None yet
Development

Successfully merging this pull request may close these issues.

PartialCellBottom needs to be adapted correctly for the GPU when using a HydrostaticFreeSurfaceModel
6 participants