Validation tests for immersed boundary implementations

[glwagner]

@johncmarshall54 suggested that a good validation test for an immersed boundary implementation is a topographic Rossby wave, a la:

https://journals.ametsoc.org/mwr/article/125/9/2293/104481/Representation-of-Topography-by-Shaved-Cells-in-a

which has an analytical solution

[weymouth]

That's a good reference. This approach is typically called a cut-cell method in more recent literature. Note the method does adjust the local coefficients applied to the pressure, and so the FFT solver would be out.

[glwagner]

Thanks @weymouth ! I was thinking a good approach might be to first implement a continuous forcing method (since its easy), along with a validation test or two; perhaps some with moving boundaries... ? Then we can assess error + accuracy, and use the validation tests to assess improvements in potential future discrete forcing implementations with modified Poisson solvers... ?

[navidcy]

This sounds like a good idea.

[apaloczy]

I've been looking at a coastal upwelling model with an analytical solution (Estrade et al., 2008's Equation 15, plotted below) that could be a good candidate for an immersed boundary validation test. It is an extension of Ekman's classical 1D solution to a simple 2D planar slope ($x-z$) geometry:

Estrade et al. test this with a 2D ROMS implementation (their Figure 10) that I've tried to replicate in Oceananigans:

The discrepancy is mostly in the BBL, so I'm wondering if cut cells (#3146) would improve it. But I'm also not sure if I'm prescribing bottom friction correctly at the immersed boundary (see code below run in Oceananigans v0.91.0). Any thoughts?

using Oceananigans
using Oceananigans.Units
using Printf

Lx = 200kilometers
dx = 100meters #200meters#500meters
hmin = 4meters
slope = 1e-3

D = 50meters
f = 3.8145e-05 # At ~15N.

dz = 2meters #4meters

Ti = 2π/f
Av = f*(D/π)^2/2
@info @sprintf("Ti: %.1f h", Ti/3600)
@info @sprintf("Av: %1.3e m2/s for D = %d m", Av, D) # Av = 4.831e-3 m2/s, D = 50 m, lat = 15N in Estrade et al. (2008).

te = 4Ti #10Ti
outdt = Ti/10
fout = "upwelling2Dhomog.nc"

logdt = outdt#Ti/100
maxcfl = 0.7

H = hmin + slope*Lx
Nx = Int(ceil(Lx/dx))
Nz = Int(ceil(H/dz))

underlying_grid = RectilinearGrid(CPU(),
                                  size=(Nx, Nz), halo=(3, 3),
                                  x = (-Lx, 0),
                                  z = (-H, 0),
                                  topology=(Bounded, Flat, Bounded))

h(x) = hmin + slope*x
grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(h))

τx₀ = 0    # [Pa]
τy₀ = -0.1 # [Pa]
ρ₀ = 1025  # [kg/m3]
r = 5e-3

Twind = Ti/2 # Wind ramp-up timescale.

τx₀ = τx₀/ρ₀
τy₀ = τy₀/ρ₀

# Boundary conditions (wind stress and bottom friction).
@inline wind_stress_u(x, t, p) = - p.τx₀ * tanh(t/p.Twind)
@inline wind_stress_v(x, t, p) = - p.τy₀ * tanh(t/p.Twind)
@inline drag_u(x, t, u, p) = - p.r * u
@inline drag_v(x, t, v, p) = - p.r * v
@inline immersed_drag_u(x, z, t, u, p) = - p.r * u
@inline immersed_drag_v(x, z, t, v, p) = - p.r * v

wind_stress_bc_u = FluxBoundaryCondition(wind_stress_u, parameters=(; τx₀, Twind))
wind_stress_bc_v = FluxBoundaryCondition(wind_stress_v, parameters=(; τy₀, Twind))

drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u), parameters=(; r))
drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:v), parameters=(; r))
immersed_drag_bc_u = FluxBoundaryCondition(immersed_drag_u, field_dependencies=(:u), parameters=(; r))
immersed_drag_bc_v = FluxBoundaryCondition(immersed_drag_v, field_dependencies=(:v), parameters=(; r))

u_bcs = FieldBoundaryConditions(top=wind_stress_bc_u, bottom=drag_bc_u, immersed=immersed_drag_bc_u)
v_bcs=FieldBoundaryConditions(top=wind_stress_bc_v, bottom=drag_bc_v, immersed=immersed_drag_bc_v)
bcs = (; u=u_bcs, v=v_bcs)

turbulence_closure = ScalarDiffusivity(VerticallyImplicitTimeDiscretization(), ν=Av)

model = HydrostaticFreeSurfaceModel(grid=grid, boundary_conditions=bcs,
                                    momentum_advection=nothing, # Linear to compare with analytical solution.
                                    closure=turbulence_closure,
                                    buoyancy=nothing,
                                    coriolis=FPlane(f=f))

g = 0.01
cgw = dx/sqrt(g*hmin)
Δt₀ = cgw*maxcfl/5
simulation = Simulation(model; Δt=Δt₀, stop_time=te)

start_time = time_ns()
function log_message(sim)
    prog = 100*time(sim)/sim.stop_time
    @printf("[%05.2f%%] t: %1.1f Ti, Δt: %s, max|u, v, w|: %1.2f, %1.2f, %1.1e m s⁻¹, walltime: %s\n",
            prog, time(sim)/Ti, prettytime(sim.Δt),
            maximum(abs, sim.model.velocities.u), maximum(abs, sim.model.velocities.v), maximum(abs, sim.model.velocities.w),
            prettytime((time_ns() - start_time) * 1e-9))
end

simulation.callbacks[:progress] = Callback(log_message, TimeInterval(logdt))
outputs = (; u=model.velocities.u, v=model.velocities.v)
simulation.output_writers[:fields2D] = NetCDFOutputWriter(model, outputs,
                                                      schedule=TimeInterval(outdt),
                                                      filename=fout,
                                                      overwrite_existing=true)

run!(simulation)

[glwagner]

This looks like a great case! Can you show the analytical and numerical solution side by side in a static plot?

[apaloczy]

This looks like a great case! Can you show the analytical and numerical solution side by side in a static plot?

Sure, here it is. Besides the BBL, I'm thinking that the offshore boundary may also be a bit off because I don't have a sponge there (like Estrade et al. have in their ROMS implementation).

[glwagner]

Here's a few thoughts:

• Is there a unit issue? It looks like drag is imposed in the numerical solution (the velocity goes to 0). But the drag seems weaker in the analytical solution for some reason.
• It looks like you are using a time-step criterion based on free surface gravity waves, but I think its likely you are using the split-explicit free surface which has no such restriction. So you can take longer time-steps (not that this would change the present discussion).
• it probably makes sense to test how the numerical solution depends on spatial resolution
• if you also have the plotting code with the analytical solution, it would be nice to have a single script that runs the simulation and produces the comparison. That way we can check to make sure that the plotting script doesn't also have bugs (and also have an easy way to run this ourselves)

[apaloczy]

I put the simulation and plotting code in the same script (see below). This time I ran for longer to get a steady flow all the way to the offshore boundary.

The wind and bottom stress units seem to check. I've tested a few different resolutions and drag coefficients (linear and quadratic) to approximate the no-slip bc in the analytical case (Estrade et al. note this sensitivity to drag coefficient in their ROMS implementation). I also added a sponge layer (not shown in the figure), but that didn't have a huge effect on the interior solution.

The near-bottom discrepancy is still there. I wonder if it's because it will take finer resolution to resolve the BBL without topography-following coordinates. Any thoughts? I haven't tested partial bottom cells yet, though. Could help with #3529.

using Oceananigans
using Oceananigans.Units
using Printf
using PyPlot
using NCDatasets
using Statistics: mean
import ColorSchemes.balance
import ColorSchemes.delta

Lx = 200kilometers
dx = 1kilometer #100meters
dz = 3meters
hmin = 4meters
slope = 1e-3
quadratic_drag = false

fout = "upwelling2D.nc"
D = 50meters
f = 3.8145e-05 # At ~15N.

Δt = 4minutes

if quadratic_drag
    cᴰ = 1e-3 # [unitless]
    @info @sprintf("Using quadratic drag, cᴰ = %.1e", cᴰ)
else
    r = 5e-3  # [m/s]
    @info @sprintf("Using linear drag, r = %.1e m/s", r)
end

Ti = 2π/f
Av = f*(D/π)^2/2
@info @sprintf("Ti: %.1f h", Ti/3600)
@info @sprintf("Av: %1.3e m2/s for D = %d m", Av, D) # Av = 4.831e-3 m2/s, D = 50 m, lat = 15N in Estrade et al. (2008).

te = 15Ti
outdt = Ti/10

logdt = Ti/100
laststeps_avg = 10

H = hmin + slope*Lx
Nx = Int(ceil(Lx/dx))
Nz = Int(ceil(H/dz))

underlying_grid = RectilinearGrid(CPU(),
                                  size=(Nx, Nz), halo=(3, 3),
                                  x = (-Lx, 0),
                                  z = (-H, 0),
                                  topology=(Bounded, Flat, Bounded))

hbot(x) = hmin + slope*x
grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(hbot))

τx₀ = 0    # [Pa]
τy₀ = -0.1 # [Pa]
ρ₀ = 1025  # [kg/m3]

Twind = Ti/10 # Wind ramp-up timescale.

# Boundary conditions (wind stress and bottom friction).
@inline wind_stress_u(x, t, p) = - (p.τx₀/p.ρ₀) * tanh(t/p.Twind)
@inline wind_stress_v(x, t, p) = - (p.τy₀/p.ρ₀) * tanh(t/p.Twind)

if quadratic_drag
    @inline drag_u(x, t, u, v, cᴰ) = - cᴰ * u * sqrt(u^2 + v^2)
    @inline drag_v(x, t, u, v, cᴰ) = - cᴰ * v * sqrt(u^2 + v^2)
    @inline immersed_drag_u(x, z, t, u, v, cᴰ) = - cᴰ * u * sqrt(u^2 + v^2)
    @inline immersed_drag_v(x, z, t, u, v, cᴰ) = - cᴰ * v * sqrt(u^2 + v^2)
    drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=cᴰ)
    drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=cᴰ)
    immersed_drag_bc_u = FluxBoundaryCondition(immersed_drag_u, field_dependencies=(:u, :v), parameters=cᴰ)
    immersed_drag_bc_v = FluxBoundaryCondition(immersed_drag_v, field_dependencies=(:u, :v), parameters=cᴰ)
else
    @inline drag_u(x, t, u, r) = - r * u
    @inline drag_v(x, t, v, r) = - r * v
    @inline immersed_drag_u(x, z, t, u, r) = - r * u
    @inline immersed_drag_v(x, z, t, v, r) = - r * v
    drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u), parameters=r)
    drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:v), parameters=r)
    immersed_drag_bc_u = FluxBoundaryCondition(immersed_drag_u, field_dependencies=(:u), parameters=r)
    immersed_drag_bc_v = FluxBoundaryCondition(immersed_drag_v, field_dependencies=(:v), parameters=r)
end

wind_stress_bc_u = FluxBoundaryCondition(wind_stress_u, parameters=(; τx₀, ρ₀, Twind))
wind_stress_bc_v = FluxBoundaryCondition(wind_stress_v, parameters=(; τy₀, ρ₀, Twind))

u_bcs = FieldBoundaryConditions(top=wind_stress_bc_u, bottom=drag_bc_u, immersed=immersed_drag_bc_u)
v_bcs=FieldBoundaryConditions(top=wind_stress_bc_v, bottom=drag_bc_v, immersed=immersed_drag_bc_v)
bcs = (; u=u_bcs, v=v_bcs)

turbulence_closure = VerticalScalarDiffusivity(ν=Av)

# Sponge layer
const sponge_width = 10kilometers
@inline sponge_mask(x, z) = exp(-(x + Lx)^2 / (2 * sponge_width^2))
uvw_sponge = Relaxation(rate=1/20minutes, mask=sponge_mask, target=0.0)
forcing_sponge = (; u=uvw_sponge, v=uvw_sponge, w=uvw_sponge)

model = HydrostaticFreeSurfaceModel(grid=grid, boundary_conditions=bcs,
                                    momentum_advection=nothing, # Linear to compare with analytical solution.
                                    closure=turbulence_closure,
                                    forcing=forcing_sponge,
                                    buoyancy=nothing,
                                    coriolis=FPlane(f=f))

simulation = Simulation(model; Δt=Δt, stop_time=te)

start_time = time_ns()
function log_message(sim)
    prog = 100*time(sim)/sim.stop_time
    @printf("[%05.2f%%] t: %1.1f Ti, Δt: %s, max|u, v, w|: %1.2f, %1.2f, %1.1e m s⁻¹, walltime: %s\n",
            prog, time(sim)/Ti, prettytime(sim.Δt),
            maximum(abs, sim.model.velocities.u), maximum(abs, sim.model.velocities.v), maximum(abs, sim.model.velocities.w),
            prettytime((time_ns() - start_time) * 1e-9))
end

simulation.callbacks[:progress] = Callback(log_message, IterationInterval(1))
outputs = (; u=model.velocities.u, v=model.velocities.v)
simulation.output_writers[:fields2D] = NetCDFOutputWriter(model, outputs,
                                                      schedule=TimeInterval(outdt),
                                                      filename=fout,
                                                      overwrite_existing=true)

run!(simulation)

# Compare equilibrated numerical solution to the steady analytical solution.

# Analytical solution to the steady problem (Estrade et al., 2008, Equation 15):
# u + i*v = (u₀ + i*v₀) * sinh(c*(z + h))/cosh(c*h) + 2*i*v₀ * ((1 - S₁) / (γ*T₁ - T₂)) * (1 - cosh(c*z)/cosh(c*h))

ds = Dataset(fout, "r")
x, z = Array(ds["xC"]), Array(ds["zC"])

nl = size(ds["u"])[end] - laststeps_avg
unum = mean(ds["u"][:, 1, :, nl:end], dims=3)[:, : , 1]
vnum = mean(ds["v"][:, 1, :, nl:end], dims=3)[:, : , 1]
unum = 0.5*(unum[2:end, :] + unum[1:end-1, :]) # Move u from xF to xC.

h = @. hmin - slope*x

ɣ = f/abs(f)
c = (1 + im*ɣ) * π/D
nz, nx = length(z), length(x)

honD = h/D
hD = π*honD
honDmi = 0

α = @. (cosh(hD) * cos(hD))^2 + (sinh(hD) * sin(hD))^2
S₁ = @. cosh(hD) * cos(hD) * 1/α
S₂ = @. sinh(hD) * sin(ɣ * hD) * 1/α
T₁ = @. cosh(hD) * sinh(hD) * 1/α
T₂ = @. cos(hD) * sin(ɣ * hD) * 1/α

# Estrade et al.'s (2008) Fig. 2
darkgray = (0.5, 0.5, 0.5)
lightgray = (0.8, 0.8, 0.8)

fig, ax = plt.subplots()
ax.plot(honD, 1 .- S₁, "k", label=L"$1 - S_1$")
ax.plot(honD, S₂, "k--", label=L"$S_2$")
ax.plot(honD, T₁ .+ T₂, color=darkgray, linewidth=2, label=L"$T_1 + T_2$")
ax.plot(honD, T₁ .- T₂, color=darkgray, linestyle="dashed", linewidth=2, label=L"$T_1 - T_2$")
ax.plot(honD, (1 .- S₁)./(T₁ .- T₂), color=lightgray, linewidth=2, label=L"$(1 - S_1)/(T_1 - T_2)$")
ax.plot(honD, S₂./(T₁ .- T₂), color=lightgray, linestyle="dashed", linewidth=2, label=L"$S_2/(T_1 - T_2)$")
ax.set_xlim(0, 1.6)
ax.set_ylim(-0.1, 1.4)
ax.invert_xaxis()
ax.legend()
ax.set_xlabel(L"$h/D$", fontsize=15)

zz = [z[j] for i=1:nx, j=1:nz]
hh = [h[i] for i=1:nx, j=1:nz]

v₀ = @. π*τy₀/(ρ₀*abs(f)*D)
u₀ = ɣ*v₀

# U = (u₀ + im*v₀) * sinh(c*(zz + hh))/cosh(c*hh) + 2j * v₀ * ((1 - S₁) / (ɣ*T₁ - T₂)) * (1 - cosh(c*zz)/cosh(c*hh)) # Eq. 15 in Estrade et al. (2008).
# Ue = (u₀ + im*v₀) * sinh(c*(zz + hh))/cosh(c*hh) - 2j * v₀ * ((1 - S₁) / (ɣ*T₁ - T₂)) * cosh(c*zz)/cosh(c*hh)

vg = @. 2*im * v₀ * (1 - S₁) / (ɣ*T₁ - T₂) # Geostrophic flow, Eq. 14 in Estrade et al. (2008).
U = @. (u₀ + im*v₀) * sinh(c*(zz + hh))/cosh(c*hh) + 2*im * v₀ * ((1 - S₁) / (ɣ*T₁ - T₂)) * (1 - cosh(c*zz)/cosh(c*hh)) # Full (frictional + geostrophic) velocity field, Eq. 15 in Estrade et al. (2008).

Ue = @. U - vg # Frictional part of the flow.
vg = imag(vg)  # vg is purely imaginary.

for i=1:nx
    fbot = z .<= - hh[i]
    U[i, fbot] .= NaN
    Ue[i, fbot] .= NaN
end

u, v = real(U), imag(U)
ue, ve = real(Ue), imag(Ue)
# Note that u = ue here (no geostrophic cross-shelf velocity).

maxu = 0.2
maxv = 0.5
numlev = 20

levsu = LinRange(-maxu, maxu, numlev)
levsv = LinRange(-maxv, maxv, numlev)

honD = @. hh/D
zonD = @. zz/D
balancecm = ColorMap(balance.colors)
deltacm = ColorMap(delta.colors)

# Analytical solution.
fig, (ax1, ax2) = plt.subplots(figsize=(14, 6), ncols=2, sharex=true, sharey=true)
cs1 = ax1.pcolormesh(honD, zonD, u, vmin=-maxu, vmax=maxu, cmap=balancecm)
cs2 = ax2.pcolormesh(honD, zonD, v, vmin=-maxv, vmax=maxv, cmap=balancecm)
plt.colorbar(cs1).set_label(L"$u$", fontsize=16)
plt.colorbar(cs2).set_label(L"$v$", fontsize=16)
ax1.contour(honD, zonD, u, levels=levsu, colors="k")
ax2.contour(honD, zonD, v, levels=levsv, colors="k")
cc1 = ax1.contour(honD, zonD, u, levels=[0], colors="g")
cc2 = ax2.contour(honD, zonD, v, levels=[0], colors="g")
ax1.plot(honD, -honD, "gray", linewidth=3)
ax2.plot(honD, -honD, "gray", linewidth=3)
ax1.set_xlim(honDmi, 2.5) # Compare with Estrade et al.'s (2008) Fig. 3
ax1.set_ylim(honDmi, -2.5) # Compare with Estrade et al.'s (2008) Fig. 3
ax1.invert_xaxis()
ax1.invert_yaxis()
ax1.set_xlabel(L"$h/D$", fontsize=15)
ax2.set_xlabel(L"$h/D$", fontsize=15)
ax1.set_ylabel(L"$z/D$", fontsize=15)
ax1.set_title("Analytical solution", fontsize=18, x=1.25)
fig.savefig("estrade_analytical.png", bbox_inches="tight", dpi=150)

# Then plot the analytical solution for the total flow next to Oceananigans' (u,v) fields.
udiff = u - unum
vdiff = v - vnum
maxudiff, maxvdiff = 0.1, 0.2

fig, ax = plt.subplots(figsize=(21, 14), ncols=3, nrows=2, sharex=true, sharey=true)
ax1, ax2, ax3, ax4, ax5, ax6 = ax

cs1 = ax1.pcolormesh(honD, zonD, u, vmin=-maxu, vmax=maxu, cmap=balancecm)
cs2 = ax2.pcolormesh(honD, zonD, v, vmin=-maxv, vmax=maxv, cmap=balancecm)
cs3 = ax3.pcolormesh(honD, zonD, unum, vmin=-maxu, vmax=maxu, cmap=balancecm)
cs4 = ax4.pcolormesh(honD, zonD, vnum, vmin=-maxv, vmax=maxv, cmap=balancecm)
cs5 = ax5.pcolormesh(honD, zonD, udiff, vmin=-maxudiff, vmax=maxudiff, cmap=deltacm)
cs6 = ax6.pcolormesh(honD, zonD, vdiff, vmin=-maxvdiff, vmax=maxvdiff, cmap=deltacm)

ax1.contour(honD, zonD, u, levels=levsu, colors="k")
ax2.contour(honD, zonD, v, levels=levsv, colors="k")
ax3.contour(honD, zonD, unum, levels=levsu, colors="k")
ax4.contour(honD, zonD, vnum, levels=levsv, colors="k")
cc1 = ax1.contour(honD, zonD, u, levels=[0], colors="g")
cc2 = ax2.contour(honD, zonD, v, levels=[0], colors="g")
cc3 = ax3.contour(honD, zonD, unum, levels=[0], colors="g")
cc4 = ax4.contour(honD, zonD, vnum, levels=[0], colors="g")

ax1.plot(honD, -honD, "gray", linewidth=3)
ax2.plot(honD, -honD, "gray", linewidth=3)
ax3.plot(honD, -honD, "gray", linewidth=3)
ax4.plot(honD, -honD, "gray", linewidth=3)
ax5.plot(honD, -honD, "gray", linewidth=3)
ax6.plot(honD, -honD, "gray", linewidth=3)

plt.colorbar(cs1).set_label(L"$u$", fontsize=16)
plt.colorbar(cs2).set_label(L"$v$", fontsize=16)
plt.colorbar(cs3).set_label(L"$u$", fontsize=16)
plt.colorbar(cs4).set_label(L"$v$", fontsize=16)
plt.colorbar(cs5).set_label(L"$u$ difference", fontsize=16)
plt.colorbar(cs6).set_label(L"$v$ difference", fontsize=16)

ax1.set_xlim(honDmi, 2.5) # Compare with Estrade et al.'s (2008) Fig. 3
ax1.set_ylim(honDmi, -2.5) # Compare with Estrade et al.'s (2008) Fig. 3
ax1.invert_xaxis()
ax1.invert_yaxis()
ax4.set_xlabel(L"$h/D$", fontsize=15)
ax1.set_ylabel(L"$z/D$", fontsize=15, y=0)
ax1.set_title("Analytical", fontsize=14)
ax3.set_title("Numerical", fontsize=14)
ax5.set_title("Analytical - numerical", fontsize=14)
fig.savefig("numerical-EMVR08_sidecmp.png", bbox_inches="tight", dpi=150)

[glwagner]

Nice, this is looking pretty good now!

I'm not sure if full or partial cells will always exhibit a form drag that prevents the two solutions from matching. But to investigate this further I think it's warranted to

1. Check the resolution dependence of the discrepancy
2. Try partial cells.

If there's still a mismatch even with refined resolution, this is a good motivation to pursue cut cells. @simone-silvestri @siddharthabishnu

[apaloczy]

Here's a twin partial cell case (dx = 1 km as in the GridFittedBottom case above). The pattern in the differences is almost identical, although their magnitude is slightly smaller with partial cells:

Resolution seems to have a negligible effect (dx = 200 m shown below, tried a few others):

The fact that partial cells actually reduce the magnitude of the difference seems encouraging, though. This could become a good validation case to illustrate the advantages of partial and cut cells for representing the BBL over a slope.

[glwagner]

That's very nice @apaloczy, thank you! I have to think about it more but I think the hypothesis is rational: there is an unavoidable error associated with step topography, but the error can be reduced by using partial cells vs full cells. If that holds this would be an excellent case to use to motivate/develop cut cells.

Note there are a few articles that suggest the same:

• How slippery are piecewise-constant coastlines in numerical ocean models?
• Sliding or Stumbling on the Staircase: Numerics of Ocean Circulation Along Piecewise-Constant Coastlines

This could be a good bona fide example I think. We only have an internal tide example now, right? But definitely a validation test otherwise.

[francispoulin]

Thanks for doing this and sharing this @apaloczy .

Is it just me or do the grid scale oscillations seem bigger in the partial cell method? I remember I wanted to test this method more a while ago and still haven't done so. I am sorry but will try and do that sometime soon.

[simone-silvestri]

The Partial cell method is bugged at the moment. @jm-c is working on correcting the bugs. An example is the output of _immersed_cell

Oceananigans.jl/src/ImmersedBoundaries/partial_cell_bottom.jl

Lines 97 to 102 in 3c4785c

	@inline function _immersed_cell(i, j, k, underlying_grid, ib::PartialCellBottom)
	# Face node above current cell
	z = znode(i, j, k+1, underlying_grid, c, c, f)
	h = @inbounds ib.bottom_height[i, j, 1]
	return z ≤ h
	end

which is incoherent with the notion of a minimum_fractional_cell_height. I.e., the test that checks whether a cell is immersed or not should include the minimal fractional cell height.
Something more along these lines

@inline function _immersed_cell(i, j, k, underlying_grid, ib::PartialCellBottom)
    # 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 
end

There are other inconsistencies as we find out that simple examples still crash also with this fix.
A simple case for which partial cells lead to the simulation crashing (independently on the time step size) is the internal tide example. If you want to contribute we are trying to correct the implementation over in #3682.

[glwagner]

For this case you can try minimum_fractional_cell_height = 0 to see if anything changes. If it doesn't then this above issue is not playing a role.

[glwagner]

@apaloczy I suspect the "bugged" that @simone-silvestri was talking about may not have affected your case in the end... so we could have tried partial cells earlier... but anyways, it could be worth trying partial cells now since we merged #3682

[apaloczy]

I've just rerun the code above with Oceananigans v0.92 using pbcs (with dx = 1 km and 200 m), and see only small differences compared to the older plots in both cases. The $v$ field is slightly less jagged now.

I also tried a run with minimum_fractional_cell_height=0 (below), and that improves the fields only very slightly compared to the default case (minimum_fractional_cell_height=0.2). If the discrepancy near the bottom is not due to the pbcs, should we expect cut cells to substantially improve things?

[glwagner]

Nice, thank you!

Minimal indeed! How is the difference in u as large as 0.1? Comparing the two solutions, it looks like they are quite close in magnitude. Is the difference computed correctly? Or perhaps this is a relative (normalized) difference?

Two things I would try are:

1. Increase vertical viscosity and
2. Increase resolution.

Honestly though, the solution looks decent to me. I hope that cut cells would improve it! In this case the boundary is linear? So cut cells should provide an "exact" representation (in the sense that as resolution is decreased, the numerical problem limits exactly to the analytical one). I'm less sure that this property holds for the grid fitted bottom, but a small resolution study would help us figure that out...

[apaloczy]

Minimal indeed! How is the difference in u as large as 0.1? Comparing the two solutions, it looks like they are quite close in magnitude. Is the difference computed correctly? Or perhaps this is a relative (normalized) difference?

So the differences are not normalized. I just evaluate the analytical expressions for $u$ and $v$ at the same locations as the numerical fields (u at xF and v at xC) and take their differences.

Resolution (from 1 km down to 10 m) doesn't improve things much, but increasing the vertical viscosity does (ignoring the sponge layer on the left boundary). Note that I'm actually changing $A_v = f(D/\pi)^2/2$ indirectly through $D$, so the horizontal and vertical coordinates also change. This is because the analytical solution is in terms of $D$.

1. Changing viscosity:

2. Changing resolution:

Honestly though, the solution looks decent to me. I hope that cut cells would improve it! In this case the boundary is linear? So cut cells should provide an "exact" representation (in the sense that as resolution is decreased, the numerical problem limits exactly to the analytical one). I'm less sure that this property holds for the grid fitted bottom, but a small resolution study would help us figure that out...

Right, the boundary is a planar slope. I think the solution looks qualitatively pretty good too, though it does not improve by choosing pbcs over a grid fitted bottom. And it doesn't look like the numerical solution converges to the analytical one with resolution, since changes are minimal even after refining by 100x.

From looking at the MITgcm docs and the references in #3123, it does seem that cut cells improve representation of flow close to topography, so this could perhaps be a good validation case for that as well. What would be needed to merge #3146? cc @siddharthabishnu @simone-silvestri

[glwagner]

What would be needed to merge #3146? cc @siddharthabishnu @simone-silvestri

I don't think there's been any work on #3146 !

It's a big effort. It would make for a pretty impactful piece of work to complete cut cells I think. You might even get a few papers out of it. But it's not something to do casually.

[glwagner]

I'm going to convert this to a discussion because there's no source code changes to make here.