(0.90.12) Bugfix for SpitExplicitFreeSurface on MultiRegionGrids and make SpitExplicitFreeSurface the default + don't allow more halos than interior points

[navidcy]

This continues a bit #3503. As summarized below there seems to be a bug with SpitExplicitFreeSurface on MultiRegionGrid. This PR fixes this and then promotes SpitExplicitFreeSurface as the default free surface solver for MultiRegionGrids.

Bugfix: set!(::HydrostaticFreeSurfaceModel) now calls initialize!(model) to set the barotropic transport that corresponds to the u, v velocities.

Closes #3513.

[navidcy]

Currently the SplitExplicit is a bit problematic. I showcase the issue below.

using Oceananigans, Logging, Test
using Oceananigans.Models.HydrostaticFreeSurfaceModels: update_state!
using Oceananigans.MultiRegion
using Oceananigans.MultiRegion: reconstruct_global_field

Logging.disable_logging(Logging.Warn)

function diffusion_cosine_test(grid; P = XPartition, regions, closure, field_name)
    model = HydrostaticFreeSurfaceModel(grid = MultiRegionGrid(grid, partition = P(regions)),
                                        closure = closure,
                                        tracers = :c,
                                        coriolis = nothing, buoyancy=nothing)

    initial_condition(x, y, z) = cos(2x)

    f = fields(model)[field_name]

    @apply_regionally set!(f, initial_condition)

    update_state!(model)

    Δt = 1e-6 * grid.Lz^2 / closure.κ

    for _ = 1:10
        time_step!(model, Δt)
    end

    return f
end

Nx = Ny = 8

grid = RectilinearGrid(size = (Nx, Ny, 1), halo = (3, 3, 3),
                       topology = (Bounded, Bounded, Bounded),
                       x = (0, 1), y = (0, 1), z = (0, 1))

closure = ScalarDiffusivity(ν = 1, κ = 1)

field_names = [:u, :v, :c]

for field_name in field_names
    println("Testing for $(string(field_name)):")

    fs = diffusion_cosine_test(grid; closure, regions = 1, field_name)
    fs = Array(interior(fs))

    f = diffusion_cosine_test(grid; closure, regions = 2, field_name)
    f = Array(interior(reconstruct_global_field(f)))

    println("  Comparing diffusion of $(string(field_name)) on single and multi-region grid:")
    @show all(f .≈ fs)
end

gives

Testing for u:
  Comparing diffusion of u on single and multi-region grid:
all(f .≈ fs) = false
Testing for v:
  Comparing diffusion of v on single and multi-region grid:
all(f .≈ fs) = false
Testing for c:
  Comparing diffusion of c on single and multi-region grid:
all(f .≈ fs) = true

Thus, the test seems to be OK for diffusing tracer but not for u or v!

cc @simone-silvestri

[navidcy]

Here's how the u field from the code-snippet above looks like when evolved on a SingleRegion and MultiRegion grids:

[simone-silvestri]

Hmmm okok looks like there's a problem there. It seems to me that also the tracer does not work is it right? Does it work with the implicit solver?

[navidcy]

Why you say the tracer doesn't work? It works!

Testing for c:
  Comparing diffusion of c on single and multi-region grid:
all(f .≈ fs) = true

And yes, this works with the implicit solver.

This is essentially one of the tests from the test_multiregion_advection_diffusion.jl... and they pass on main using the UnifiedImplicitSolver with MultiRegionGrid.

[navidcy]

Perhaps there is an @apply_regionally that is missing from somewhere?
The fact that tracers are OK but velocities are not hints where the problem might be coming from?

[simone-silvestri]

Ok, sorry I didn't register that true. So it might be a halo passing problem: if the corners are not communicated properly the results would be wrong. I don't think it is an @apply_regionally problem. If that were the problem we would get an error somewhere because all the data structures would be wrong.

[navidcy]

If you want we can Zoom and try to sort it out. No rush -- whenever.

[navidcy]

So @simone-silvestri, the tests now pass but I think it's totally wrong. It just happens that whatever wrong thing happens happens both in serial and multi-region....

using Oceananigans, Logging, Test
using Oceananigans.Models.HydrostaticFreeSurfaceModels: update_state!
using Oceananigans.MultiRegion
using Oceananigans.MultiRegion: reconstruct_global_field

using GLMakie

Logging.disable_logging(Logging.Warn)

function initialize_model(grid; P = XPartition, regions, closure, field_name)

    model = HydrostaticFreeSurfaceModel(; grid = MultiRegionGrid(grid, partition = P(regions)),
                                          free_surface = SplitExplicitFreeSurface(substeps=5),
                                          closure,
                                          tracers = :c,
                                          coriolis = nothing,
                                          buoyancy=nothing)

    initial_condition(x, y, z) = cos(2x)

    f = fields(model)[field_name]

    @apply_regionally set!(f, initial_condition)

    update_state!(model)

    return model
end

Nx = Ny = 32

grid = RectilinearGrid(size = (Nx, Ny, 1), halo = (3, 3, 3),
                       topology = (Bounded, Bounded, Bounded),
                       x = (0, 1), y = (0, 1), z = (0, 1))

closure = ScalarDiffusivity(ν = 1, κ = 1)

field_name = :u

model_serial   = initialize_model(grid; P = XPartition, regions = 1, field_name, closure)
model_parallel = initialize_model(grid; P = XPartition, regions = 2, field_name, closure)

function plot_comparison(models, field_name, title)

    model_serial, model_parallel = models
    f_serial   = fields(model_serial)[field_name]
    f_parallel = reconstruct_global_field(fields(model_parallel)[field_name])

    fig = Figure()
    ax1 = Axis(fig[1, 1])
    ax2 = Axis(fig[2, 1])

    x, y, z = nodes(f_serial)
    hm1 = heatmap!(ax1, x, y, interior(f_serial, :, :, 1))
    Colorbar(fig[1, 2], hm1)

    x, y, z = nodes(f_parallel)
    hm2 = heatmap!(ax2, x, y, interior(f_parallel, :, :, 1))
    Colorbar(fig[2, 2], hm2)

    fig[0, :] = Label(fig, title, fontsize=20, tellwidth=false)

    return fig
end

fig = plot_comparison((model_serial, model_parallel), field_name, "before timestep")

save("before_timestep.png", fig)

Δt = 1e-6 * grid.Lz^2 / closure.κ

for _ = 1:1
    time_step!(model_serial, Δt)
    time_step!(model_parallel, Δt)
end

fig = plot_comparison((model_serial, model_parallel), field_name, "after 1 timestep")

save("after_1_timestep.png", fig)