Fix test for implicit free surface with `ImmersedBoundary` and indexed `ReducedField`s #2723

navidcy · 2022-09-03T22:14:20Z

There is a test for immersed boundary inside test_implicit_free_surface.jl but it seems to be trivially passing (see, e.g., the log from tests run on main).

One problem was that set_simple_divergent_velocity! function defined at

Oceananigans.jl/test/test_implicit_free_surface_solver.jl

Lines 17 to 35 in af21542

    
           function set_simple_divergent_velocity!(model) 
        
               # Create a divergent velocity 
        
               grid = model.grid 
        
               u, v, w = model.velocities 
        
               η = model.free_surface.η 
        
               u .= 0 
        
               v .= 0 
        
               η .= 0 
        
               imid = Int(floor(grid.Nx / 2)) + 1 
        
               jmid = Int(floor(grid.Ny / 2)) + 1 
        
               CUDA.@allowscalar u[imid, jmid, 1] = 1 
        
               update_state!(model) 
        
               return nothing 
        
           end

was setting a non-zero u velocity in the center of the domain at k = 1, but that was inside the immersed boundary. I changed the set_simple_divergent_velocity! function to address that. However, still there is an issue...

I managed to pinpoint it down to the step:

Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl

Line 140 in f585ef8

    
           compute_implicit_free_surface_right_hand_side!(rhs, solver, g, Δt, ∫ᶻQ, η)

Here's a demonstration... (not a MWE). I have a model with non-zero u velocity:

julia> model.velocities.u
129×1×5 Field{Face, Center, Center} on ImmersedBoundaryGrid on CPU
├── grid: 128×1×5 ImmersedBoundaryGrid{Float64, Bounded, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Open, east: Open, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 135×7×11 OffsetArray(::Array{Float64, 3}, -2:132, -2:4, -2:8) with eltype Float64 with indices -2:132×-2:4×-2:8
    └── max=0.1, min=0.0, mean=0.000159236

julia> model.velocities.u[imid, 1, 3]
0.1

julia> ∫ᶻQ = model.free_surface.barotropic_volume_flux
NamedTuple with 2 Fields on 128×1×5 ImmersedBoundaryGrid{Float64, Bounded, Periodic, Bounded} on CPU with 3×3×3 halo:
├── u: 129×1×1 Field{Face, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
└── v: 128×1×1 Field{Center, Face, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU

julia> ∫ᶻQ.u .= 0; ∫ᶻQ.v .= 0;

I call

julia> sum!(∫ᶻQ.u, Ax * model.velocities.u);

Now, when I print out ∫ᶻQ.u I see that all its elements are supposedly 0:

julia> ∫ᶻQ.u
129×1×1 Field{Face, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
├── grid: 128×1×5 ImmersedBoundaryGrid{Float64, Bounded, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 135×7×1 OffsetArray(::Array{Float64, 3}, -2:132, -2:4, 1:1) with eltype Float64 with indices -2:132×-2:4×1:1
    └── max=0.0, min=0.0, mean=0.0

But this returns the right answer

julia> maximum(∫ᶻQ.u.data)
8.0

The velocity set by set_simple_divergent_velocity! is 0.1 and Δz = 80 so that's correct!

Hm.... Here's a MWE now to demonstrate that somewhere something is not passed on.

using Oceananigans

using Statistics
using Oceananigans.Units
using Oceananigans.Architectures: device_event
using Oceananigans.TimeSteppers: update_state!
using LinearAlgebra: norm

using Oceananigans.Models.HydrostaticFreeSurfaceModels:
    ImplicitFreeSurface,
    FreeSurface,
    PCGImplicitFreeSurfaceSolver,
    implicit_free_surface_step!

function set_simple_divergent_velocity!(model)
    # Create a divergent velocity
    grid = model.grid

    u, v, w = model.velocities
    η = model.free_surface.η

    u .= 0
    v .= 0
    η .= 0

    imid = Int(floor(grid.Nx / 2)) + 1
    jmid = Int(floor(grid.Ny / 2)) + 1

    k_index = 1
 
    grid isa ImmersedBoundaryGrid && begin
        while grid.immersed_boundary.bottom_height[imid, jmid] > grid.underlying_grid.zᵃᵃᶜ[k_index]
            k_index += 1
        end
    end

    k_index = k_index + 1 ≤ grid.Nz ? k_index + 1 : k_index
    
    u[imid, jmid, k_index] = 0.1

    update_state!(model)

    return nothing
end


arch = CPU()

rectilinear_grid = RectilinearGrid(arch, size = (128, 1, 5),
                                    x = (-500kilometers, 500kilometers),
                                    y = (0, 1),
                                    z = (-400, 0),
                                    topology = (Bounded, Periodic, Bounded))

Lz = rectilinear_grid.Lz
width = 50kilometers
bump(x, y) = - Lz * (1 - 0.2 * exp(-x^2 / 2width^2))

grid = ImmersedBoundaryGrid(rectilinear_grid, GridFittedBottom(bump))

free_surface = ImplicitFreeSurface(solver_method=:PreconditionedConjugateGradient,
                                   abstol=1e-15, reltol=0)

model = HydrostaticFreeSurfaceModel(; grid,
                                    momentum_advection = nothing,
                                    free_surface)

set_simple_divergent_velocity!(model)

@show model.velocities.u

events = ((device_event(arch), device_event(arch)), (device_event(arch), device_event(arch)))

Δt = 900.0
implicit_free_surface_step!(model.free_surface, model, Δt, 1.5, events)

η = model.free_surface.η
@info "implicit free surface solver test, norm(η): $(norm(η)), maximum(abs, η): $(maximum(abs, η))"

gives

julia> @show model.velocities.u
model.velocities.u = 129×1×5 Field{Face, Center, Center} on ImmersedBoundaryGrid on CPU
├── grid: 128×1×5 ImmersedBoundaryGrid{Float64, Bounded, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Open, east: Open, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 135×7×11 OffsetArray(::Array{Float64, 3}, -2:132, -2:4, -2:8) with eltype Float64 with indices -2:132×-2:4×-2:8
    └── max=0.1, min=0.0, mean=0.000159236

julia> @info "implicit free surface solver test, norm(η): $(norm(η)), maximum(abs, η): $(maximum(abs, η))"
[ Info: implicit free surface solver test, norm(η): 0.0, maximum(abs, η): 0.0

cc @glwagner, @simone-silvestri

navidcy · 2022-09-03T23:56:16Z

@elise-palethorpe, fyi.

glwagner · 2022-09-06T18:48:28Z

I'll try to parse this carefully, but I do have a comment on this part:

Now, when I print out ∫ᶻQ.u I see that all its elements are supposedly 0:

julia> ∫ᶻQ.u
129×1×1 Field{Face, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
├── grid: 128×1×5 ImmersedBoundaryGrid{Float64, Bounded, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 135×7×1 OffsetArray(::Array{Float64, 3}, -2:132, -2:4, 1:1) with eltype Float64 with indices -2:132×-2:4×1:1
    └── max=0.0, min=0.0, mean=0.0

But this returns the right answer

julia> maximum(∫ᶻQ.u.data)
8.0

Crucially, maximum(∫ᶻQ.u) will ignore immersed cells and halo cells, but maximum(∫ᶻQ.u.data) does not. Does this explain what you're seeing?

simone-silvestri · 2022-09-06T18:53:29Z

One problem is that ∫ᶻQ.u is a reduced field so maximum(∫ᶻQ.u) excludes points at i, j corresponding to immersed cells at i, j, k = 1

test/test_implicit_free_surface_solver.jl

navidcy · 2022-09-06T20:13:52Z

One problem is that ∫ᶻQ.u is a reduced field so maximum(∫ᶻQ.u) excludes points at i, j corresponding to immersed cells at i, j, k = 1

I'm wasn't setting the u velocity at i, j, k = 1, 1, 1 but rather at i, j, k = 45, 1, 3 and now, after 1638f39, at i, j, k = 45, 1, 5.

Still I get

julia> η = model.free_surface.η
128×1×1 Field{Center, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
├── grid: 128×1×5 ImmersedBoundaryGrid{Float64, Bounded, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: ZeroFlux, east: ZeroFlux, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 134×7×1 OffsetArray(::Array{Float64, 3}, -2:131, -2:4, 1:1) with eltype Float64 with indices -2:131×-2:4×1:1
    └── max=0.0, min=0.0, mean=0.0

julia> @info "implicit free surface solver test, norm(η): $(norm(η)), maximum(abs, η): $(maximum(abs, η))"
[ Info: implicit free surface solver test, norm(η): 0.0, maximum(abs, η): 0.0

glwagner · 2022-09-07T16:08:39Z

I'm wasn't setting the u velocity at i, j, k = 1, 1, 1 but rather at i, j, k = 45, 1, 3

Is i, j, k = 45, 1, 3 within the immersed boundary or halo region?

navidcy · 2022-09-07T20:06:18Z

No

glwagner · 2022-09-08T15:23:22Z

test/test_implicit_free_surface_solver.jl

+    grid isa ImmersedBoundaryGrid &&
+        grid.immersed_boundary.bottom_height[i, j] ≥ grid.underlying_grid.zᵃᵃᶜ[k] &&
+            @error "The point to set u is on land!"


Would something like

inactive_cell(i, j, k, grid) && error("The nudged cell at ($i, $j, $k) is inactive.")

work?

This catches both halos and immersed cells, and doesn't require so many explicit references to grid properties that might have to be changed / maintained in the future if we change the names of those properties.

That's great suggestion. I did it!

The problem persists though: the cell I'm nudging is not inactive but still I get η=0....

navidcy · 2022-09-14T21:46:20Z

An update:

using Oceananigans

using Statistics
using Oceananigans.Units
using Oceananigans.Architectures: device_event
using Oceananigans.TimeSteppers: update_state!
using LinearAlgebra: norm

using Oceananigans.Models.HydrostaticFreeSurfaceModels:
    ImplicitFreeSurface,
    FreeSurface,
    PCGImplicitFreeSurfaceSolver,
    implicit_free_surface_step!

function set_simple_divergent_velocity!(model)
    # Create a divergent velocity
    grid = model.grid

    u, v, w = model.velocities
    η = model.free_surface.η

    u .= 0
    v .= 0
    η .= 0

    # pick a surface cell at the middle of the domain
    i, j, k = Int(floor(grid.Nx / 2)) + 1, Int(floor(grid.Ny / 2)) + 1, grid.Nz

    inactive_cell(i, j, k, grid) && error("The nudged cell at ($i, $j, $k) is inactive.")

    if grid isa RectilinearGrid
        Δy = grid.Δyᵃᶜᵃ
    end

    if grid isa LatitudeLongitudeGrid
        Δy = grid.Δyᶜᶠᵃ
    end

    if grid isa ImmersedBoundaryGrid
        if grid isa ImmersedBoundaryGrid && grid.underlying_grid isa RectilinearGrid
            Δy = grid.underlying_grid.Δyᵃᶜᵃ
        elseif grid.underlying_grid isa LatitudeLongitudeGrid
            Δy = grid.underlying_grid.Δyᶜᶠᵃ
        end
    end

    Δz = CUDA.@allowscalar grid.Δzᵃᵃᶜ

    # We prescribe the value of the zonal transport in a cell, i.e., `u * Δy * Δz`. This
    # way `norm(rhs)` of the free-surface solver does not depend on the grid extensd/resolution.
    transport = 1e5 # m³ s⁻¹
    CUDA.@allowscalar u[i, j, k] = transport / (Δy * Δz)

    update_state!(model)

    return nothing
end


arch = CPU()

rectilinear_grid = RectilinearGrid(arch, size = (128, 1, 5),
x = (-5000kilometers, 5000kilometers),
y = (0, 100kilometers),
z = (-500, 0),
topology = (Bounded, Periodic, Bounded))

Lz = rectilinear_grid.Lz
width = rectilinear_grid.Lx / 20

bump(x, y) = - Lz * (1 - 0.2 * exp(-x^2 / 2width^2))

grid = ImmersedBoundaryGrid(rectilinear_grid, GridFittedBottom(bump))

free_surface = ImplicitFreeSurface(solver_method=:PreconditionedConjugateGradient,
                                   abstol=1e-15, reltol=0)

model = HydrostaticFreeSurfaceModel(; grid,
                                    momentum_advection = nothing,
                                    free_surface)

set_simple_divergent_velocity!(model)

@show model.velocities.u

events = ((device_event(arch), device_event(arch)), (device_event(arch), device_event(arch)))

Δt = 900.0
implicit_free_surface_step!(model.free_surface, model, Δt, 1.5, events)

η = model.free_surface.η

@show norm(η)
@show maximum(abs, η)
@show maximum(abs, η.data)

gives

julia> @show maximum(abs, model.velocities.u)
maximum(abs, model.velocities.u) = 0.01
0.01

julia> @show norm(η)
norm(η) = 0.0
0.0

julia> @show maximum(abs, η)
maximum(abs, η) = 0.0
0.0

julia> @show maximum(abs, η.data)
maximum(abs, η.data) = 0.0
0.0

navidcy · 2022-09-14T21:55:39Z

The problem stems (I think) that after

Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl

Line 140 in 0aa674f

    
           compute_implicit_free_surface_right_hand_side!(rhs, solver, g, Δt, ∫ᶻQ, η)

maximum(rhs) = 0
maximum(rhs.data) = 11.330180144199204

and then when rhs is given to solve! in

Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl

Line 145 in 0aa674f

solve!(η, solver, rhs, g, Δt)

solve returns an η all zeros!

maximum(η.data) = 0.0

(cc @glwagner, @simone-silvestri)

navidcy · 2022-09-15T02:51:44Z

most tests fail now...... :(

Also I get

julia> ∫ᶻ_Axᶠᶜᶜ = Field{Face, Center, Nothing}(with_halo((3, 3, 1), grid), indices = grid.Nz)
ERROR: BoundsError: attempt to access Tuple{Int64} at index [2]
Stacktrace:
 [1] indexed_iterate
   @ ./tuple.jl:86 [inlined]
 [2] new_data(FT::DataType, grid::ImmersedBoundaryGrid{Float64, Bounded, Periodic, Bounded, RectilinearGrid{Float64, Bounded, Periodic, Bounded, Float64, Float64, Float64, OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, CPU}, GridFittedBottom{OffsetMatrix{Float64, Matrix{Float64}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, CPU}, loc::Tuple{DataType, DataType, DataType}, indices::Int64)
   @ Oceananigans.Grids ~/Research/OC2.jl/src/Grids/new_data.jl:59
 [3] (Field{Face, Center, Nothing, O, G, I, D, T, B, S, F} where {O, G, I, D, T, B, S, F})(grid::ImmersedBoundaryGrid{Float64, Bounded, Periodic, Bounded, RectilinearGrid{Float64, Bounded, Periodic, Bounded, Float64, Float64, Float64, OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, CPU}, GridFittedBottom{OffsetMatrix{Float64, Matrix{Float64}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, CPU}, T::DataType; kw::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:indices,), Tuple{Int64}}})
   @ Oceananigans.Fields ~/Research/OC2.jl/src/Fields/field.jl:158
 [4] top-level scope
   @ REPL[24]:1
 [5] top-level scope
   @ ~/.julia/packages/CUDA/DfvRa/src/initialization.jl:52

simone-silvestri · 2022-09-15T03:05:17Z

that should be

julia> ∫ᶻ_Axᶠᶜᶜ = Field{Face, Center, Nothing}(with_halo((3, 3, 1), grid), indices = (:, :, grid.Nz))
4×4×1 Field{Face, Center, Nothing} reduced over dims = (3,) on RectilinearGrid on CPU
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×1 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 10×10×1 OffsetArray(::Array{Float64, 3}, -2:7, -2:7, 4:4) with eltype Float64 with indices -2:7×-2:7×4:4
    └── max=0.0, min=0.0, mean=0.0

there are a couple of things to smooth out

simone-silvestri · 2022-09-15T03:07:33Z

This is the progress

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

julia> immersed_grid =  ImmersedBoundaryGrid(grid, GridFittedBottom((x, y) -> -0.5))
4×4×4 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo:
├── immersed_boundary: GridFittedBottom(min(h)=-5.00e-01, max(h)=-5.00e-01)
├── underlying_grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0)  regularly spaced with Δx=0.25
├── Periodic y ∈ [0.0, 1.0)  regularly spaced with Δy=0.25
└── Bounded  z ∈ [-1.0, 0.0] regularly spaced with Δz=0.25

julia> set!(Field((Center, Center, Nothing), immersed_grid, indices = (:, :, 4)), (x, y) -> x)
4×4×1 Field{Center, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
├── grid: 4×4×4 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 10×10×1 OffsetArray(::Array{Float64, 3}, -2:7, -2:7, 4:4) with eltype Float64 with indices -2:7×-2:7×4:4
    └── max=0.875, min=0.125, mean=Inf

julia> set!(Field((Center, Center, Nothing), immersed_grid, indices = (:, :, 1)), (x, y) -> x)
4×4×1 Field{Center, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
├── grid: 4×4×4 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 10×10×1 OffsetArray(::Array{Float64, 3}, -2:7, -2:7, 1:1) with eltype Float64 with indices -2:7×-2:7×1:1
    └── max=-Inf, min=Inf, mean=NaN

I am still struggling with the mean though

…into ncc/test-for-2710

simone-silvestri · 2022-09-15T04:38:38Z

the boundary conditions will be a bit tricky to handle since we pass only the data and not the field... but we can fix it

navidcy · 2022-09-15T07:29:56Z

now norm seems that it's not working for simple RectilinearGrid :)
E.g., in the pgc solver, residual_norm = norm(solver.residual) gives zero... even when q=0 but rhs $\ne$ 0.

simone-silvestri · 2022-09-15T10:45:00Z

weird, I get this

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

julia> immersed_grid =  ImmersedBoundaryGrid(grid, GridFittedBottom((x, y) -> -0.5))
4×4×4 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo:
├── immersed_boundary: GridFittedBottom(min(h)=-5.00e-01, max(h)=-5.00e-01)
├── underlying_grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0)  regularly spaced with Δx=0.25
├── Periodic y ∈ [0.0, 1.0)  regularly spaced with Δy=0.25
└── Bounded  z ∈ [-1.0, 0.0] regularly spaced with Δz=0.25

julia> f1 = set!(Field((Center, Center, Nothing), immersed_grid, indices = (:, :, 4)), (x, y) -> x)
4×4×1 Field{Center, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
├── grid: 4×4×4 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 10×10×1 OffsetArray(::Array{Float64, 3}, -2:7, -2:7, 4:4) with eltype Float64 with indices -2:7×-2:7×4:4
    └── max=0.875, min=0.125, mean=0.5

julia> f2 = set!(Field((Center, Center, Nothing), immersed_grid, indices = (:, :, 1)), (x, y) -> x)
4×4×1 Field{Center, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
├── grid: 4×4×4 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 10×10×1 OffsetArray(::Array{Float64, 3}, -2:7, -2:7, 1:1) with eltype Float64 with indices -2:7×-2:7×1:1
    └── max=-Inf, min=Inf, mean=NaN

julia> using Statistics: norm

julia> norm(f1)
2.29128784747792

julia> norm(f2)
0.0

Maybe we have not set correctly the index of the residuals?

…into ncc/test-for-2710

glwagner · 2022-09-15T13:16:48Z

weird, I get this

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

julia> immersed_grid =  ImmersedBoundaryGrid(grid, GridFittedBottom((x, y) -> -0.5))
4×4×4 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo:
├── immersed_boundary: GridFittedBottom(min(h)=-5.00e-01, max(h)=-5.00e-01)
├── underlying_grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0)  regularly spaced with Δx=0.25
├── Periodic y ∈ [0.0, 1.0)  regularly spaced with Δy=0.25
└── Bounded  z ∈ [-1.0, 0.0] regularly spaced with Δz=0.25

julia> f1 = set!(Field((Center, Center, Nothing), immersed_grid, indices = (:, :, 4)), (x, y) -> x)
4×4×1 Field{Center, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
├── grid: 4×4×4 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 10×10×1 OffsetArray(::Array{Float64, 3}, -2:7, -2:7, 4:4) with eltype Float64 with indices -2:7×-2:7×4:4
    └── max=0.875, min=0.125, mean=0.5

julia> f2 = set!(Field((Center, Center, Nothing), immersed_grid, indices = (:, :, 1)), (x, y) -> x)
4×4×1 Field{Center, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
├── grid: 4×4×4 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 10×10×1 OffsetArray(::Array{Float64, 3}, -2:7, -2:7, 1:1) with eltype Float64 with indices -2:7×-2:7×1:1
    └── max=-Inf, min=Inf, mean=NaN

julia> using Statistics: norm

julia> norm(f1)
2.29128784747792

julia> norm(f2)
0.0

Maybe we have not set correctly the index of the residuals?

What is expected? Isn't indices = (:, :, 1) below the bottom for f2?

simone-silvestri · 2022-09-15T13:51:27Z

Yeah so this should be ok.
The general idea is that ReducedField have an Integer index in field.indices at the position at which they are reduced.
Also the underlying OffsetArray will be have an offset reflecting the position of the reduced plane in the grid
In this way a reduced field will have a "position" with respect to the full 3D grid.

Reductions and reduced fields without a specified position are located at index = 1 by default

I am still bugfixing a bit though.
A problem that I am foreseeing is that the GPU has no information about either indices or reduced locations but only the data structure, maybe we should ship the indices with the data on the GPU?
Also, since the underlying data is also offset we could retrieve the information on the GPU although it will probably get a bit wonky

glwagner · 2022-09-15T15:01:12Z

Yeah so this should be ok. The general idea is that ReducedField have an Integer index in field.indices at the position at which they are reduced. Also the underlying OffsetArray will be have an offset reflecting the position of the reduced plane in the grid In this way a reduced field will have a "position" with respect to the full 3D grid.

I don't totally understand. I think we want to build out support for behavior when a field is sliced at a single index. But why would we want to build this feature for fields that are also reduced?

Reductions and reduced fields without a specified position are located at index = 1 by default

I'm opposed to this, which reduces the concepts we can express in the code. I think we retain more flexibility if we interpret a Nothing location to mean that a field has no location in the reduced direction. This is a valid mathematical concept. Otherwise, we cannot distinguish between slices and true two dimensional fields.

I think in fact that it should be invalid to specify the index of a field in the reduced direction.

glwagner · 2022-09-15T15:06:16Z

src/Fields/field.jl


 function validate_index(idx::UnitRange, loc, topo, N, H)
    all_idx = all_indices(loc, topo, N, H)
-    (first(idx) ∈ all_idx && last(idx) ∈ all_idx) || throw(ArgumentError("The indices $idx must slice $I"))
+    (first(idx) ∈ all_idx && last(idx) ∈ all_idx) || throw(ArgumentError("The indices $idx must be a slice of the domain"))


The original message was more specific and clear. Can we return to that?

but what was I?

glwagner · 2022-09-15T15:14:34Z

All of that said, I suspect that the changes in this PR are important --- it's just that we need to clarify the difference between a Field that is sliced in a certain direction, and a Field with a Nothing location. My only point is that I don't think we should conflate these concepts. I can't figure out if this PR is conflating these or not.

One consequence is: the free surface will not be a reduced field.

This design was the primary motivation for PR #2246.

simone-silvestri · 2022-09-15T17:55:32Z

Ok, I agree with your comments.
This said then, it is better to close this PR and start from a clean slate.
I'll close this PR and open a new one which will focus on the following changes:

Free surface will be: Field((Center, Center, Face), grid, indices = (:, :, grid.Nz+1))
Changes to inner workings of PCG and MG solver to accommodate this change
Immersed reductions for sliced and windowed fields (and in general conditional reductions)
Try windowed-sliced BC? (We need this for η)
Reductions on ReducedFields

glwagner · 2022-09-15T20:49:45Z

That will be so so beautiful

update for grid with immersed boundary

f585ef8

navidcy added bug 🐞 Even a perfect program still has bugs testing 🧪 Tests get priority in case of emergency evacuation labels Sep 3, 2022

navidcy requested review from glwagner and simone-silvestri September 4, 2022 05:33

simone-silvestri reviewed Sep 6, 2022

View reviewed changes

test/test_implicit_free_surface_solver.jl Outdated Show resolved Hide resolved

simplify picking a point w non-zero u velocity

1638f39

Update test_implicit_free_surface_solver.jl

ece5e79

glwagner reviewed Sep 8, 2022

View reviewed changes

glwagner approved these changes Sep 8, 2022

View reviewed changes

navidcy added 4 commits September 9, 2022 07:09

check for inactive_cell

314858b

adds debugging

0128ca2

merge main

cbada59

cleanup

0aa674f

simone-silvestri added 3 commits September 14, 2022 21:16

changed indices for reduced fields

e624b10

fixing set!

f10b415

should be ok now

d31a94f

simone-silvestri added 2 commits September 14, 2022 23:08

introduced OneFieldGridded

71031f7

fixed the mean!

be8eea7

some cleaning...

5952b4d

simone-silvestri changed the title ~~Fix test for implicit free surface with ImmersedBoundary~~ Fix test for implicit free surface with ImmersedBoundary and indexed for ReducedFields Sep 15, 2022

simone-silvestri changed the title ~~Fix test for implicit free surface with ImmersedBoundary and indexed for ReducedFields~~ Fix test for implicit free surface with ImmersedBoundary and indexed ReducedFields Sep 15, 2022

navidcy and others added 3 commits September 15, 2022 14:17

fix unitrange construct

e6b064b

only boundary conditions to fix

0398ba3

Merge branch 'ncc/test-for-2710' of github.com:CliMA/Oceananigans.jl …

da244f4

…into ncc/test-for-2710

simone-silvestri added 2 commits September 15, 2022 00:58

fixed boundary condition issue

e855183

comments...

00ee7e6

navidcy added 3 commits September 15, 2022 17:33

fix docstring

98626e2

delete stray space

b5fb600

fix indices

c889495

simone-silvestri added 2 commits September 15, 2022 07:49

fixed field_reductions

1f1cd4d

Merge branch 'ncc/test-for-2710' of github.com:CliMA/Oceananigans.jl …

77f4e20

…into ncc/test-for-2710

fixed abstract operation tests

f6adf70

glwagner reviewed Sep 15, 2022

View reviewed changes

simone-silvestri closed this Sep 15, 2022

simone-silvestri mentioned this pull request Sep 15, 2022

BC and Immersed Reductions for Sliced Fields #2740

Merged

navidcy deleted the ncc/test-for-2710 branch February 8, 2023 04:01

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix test for implicit free surface with `ImmersedBoundary` and indexed `ReducedField`s #2723

Fix test for implicit free surface with `ImmersedBoundary` and indexed `ReducedField`s #2723

navidcy commented Sep 3, 2022 •

edited

Loading

navidcy commented Sep 3, 2022

glwagner commented Sep 6, 2022

simone-silvestri commented Sep 6, 2022

navidcy commented Sep 6, 2022 •

edited

Loading

glwagner commented Sep 7, 2022

navidcy commented Sep 7, 2022

glwagner Sep 8, 2022

glwagner Sep 8, 2022

navidcy Sep 8, 2022

navidcy Sep 8, 2022

navidcy commented Sep 14, 2022

navidcy commented Sep 14, 2022

navidcy commented Sep 15, 2022

simone-silvestri commented Sep 15, 2022

simone-silvestri commented Sep 15, 2022

simone-silvestri commented Sep 15, 2022

navidcy commented Sep 15, 2022 •

edited

Loading

simone-silvestri commented Sep 15, 2022 •

edited

Loading

glwagner commented Sep 15, 2022

simone-silvestri commented Sep 15, 2022 •

edited

Loading

glwagner commented Sep 15, 2022 •

edited

Loading

glwagner Sep 15, 2022

navidcy Sep 15, 2022

glwagner commented Sep 15, 2022 •

edited

Loading

simone-silvestri commented Sep 15, 2022

glwagner commented Sep 15, 2022

	function set_simple_divergent_velocity!(model)
	# Create a divergent velocity
	grid = model.grid

	u, v, w = model.velocities
	η = model.free_surface.η

	u .= 0
	v .= 0
	η .= 0

	imid = Int(floor(grid.Nx / 2)) + 1
	jmid = Int(floor(grid.Ny / 2)) + 1
	CUDA.@allowscalar u[imid, jmid, 1] = 1

	update_state!(model)

	return nothing
	end

Fix test for implicit free surface with ImmersedBoundary and indexed ReducedFields #2723

Fix test for implicit free surface with ImmersedBoundary and indexed ReducedFields #2723

Conversation

navidcy commented Sep 3, 2022 • edited Loading

navidcy commented Sep 3, 2022

glwagner commented Sep 6, 2022

simone-silvestri commented Sep 6, 2022

navidcy commented Sep 6, 2022 • edited Loading

glwagner commented Sep 7, 2022

navidcy commented Sep 7, 2022

glwagner Sep 8, 2022

Choose a reason for hiding this comment

glwagner Sep 8, 2022

Choose a reason for hiding this comment

navidcy Sep 8, 2022

Choose a reason for hiding this comment

navidcy Sep 8, 2022

Choose a reason for hiding this comment

navidcy commented Sep 14, 2022

navidcy commented Sep 14, 2022

navidcy commented Sep 15, 2022

simone-silvestri commented Sep 15, 2022

simone-silvestri commented Sep 15, 2022

simone-silvestri commented Sep 15, 2022

navidcy commented Sep 15, 2022 • edited Loading

simone-silvestri commented Sep 15, 2022 • edited Loading

glwagner commented Sep 15, 2022

simone-silvestri commented Sep 15, 2022 • edited Loading

glwagner commented Sep 15, 2022 • edited Loading

glwagner Sep 15, 2022

Choose a reason for hiding this comment

navidcy Sep 15, 2022

Choose a reason for hiding this comment

glwagner commented Sep 15, 2022 • edited Loading

simone-silvestri commented Sep 15, 2022

glwagner commented Sep 15, 2022

Fix test for implicit free surface with `ImmersedBoundary` and indexed `ReducedField`s #2723

Fix test for implicit free surface with `ImmersedBoundary` and indexed `ReducedField`s #2723

navidcy commented Sep 3, 2022 •

edited

Loading

navidcy commented Sep 6, 2022 •

edited

Loading

navidcy commented Sep 15, 2022 •

edited

Loading

simone-silvestri commented Sep 15, 2022 •

edited

Loading

simone-silvestri commented Sep 15, 2022 •

edited

Loading

glwagner commented Sep 15, 2022 •

edited

Loading

glwagner commented Sep 15, 2022 •

edited

Loading