`PerturbationAdvectionOpenBoundaryCondition` leads to nonzero mass fluxes into the domain

[tomchor]

It seems that the PerturbationAdvection matching scheme (and probably FlatExtrapolation as well) creates a nonzero net mass flux in the domain. As a demonstration, in the following MWE a flow passes through a RectilinearGrid using PerturbationAdvection in both the west and east boundary conditions. Importantly, the simulation keeps track of and prints the net mass flux through the domain calculated by integrating u at the west and east boundaries.

using Oceananigans
using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition
using Oceananigans.Diagnostics: AdvectiveCFL
using Printf

Δx = Δz = 0.05
Nx = Nz = round(Int, 2 / Δx)
grid = RectilinearGrid(size = (Nx, Nz), halo = (4, 4), extent = (1, 1),
                       topology = (Bounded, Flat, Bounded))

U₀ = 1.0
inflow_timescale = 1e-2
outflow_timescale = Inf
u_boundary_conditions = FieldBoundaryConditions(west = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale),
                                                east = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale))
boundary_conditions = (; u = u_boundary_conditions,)

model = NonhydrostaticModel(; grid, boundary_conditions,
                            timestepper = :RungeKutta3,)

uᵢ(x, z) = U₀ + 1e-4 * rand() 
set!(model, u=uᵢ)
simulation = Simulation(model; Δt=0.1Δx/U₀, stop_time=1, verbose=false)

cfl_calculator = AdvectiveCFL(simulation.Δt)
function progress(sim)
    u, v, w = model.velocities
    cfl_value = cfl_calculator(model)
    west_mass_flux = Field(Average(view(u, 1, :, :)))[]
    east_mass_flux = Field(Average(view(u, grid.Nx+1, :, :)))[]
    net_mass_flux = east_mass_flux - west_mass_flux
    @info @sprintf("time: %.3f, max|u|: %.3f, CFL: %.4f, Net mass flux: %.4e",
                   time(sim), maximum(abs, u), cfl_value, net_mass_flux)
end
add_callback!(simulation, progress, IterationInterval(20))
run!(simulation)

The script produces:

[ Info: time: 0.000, max|u|: 1.000, CFL: 0.2000, Net mass flux: 3.4537e-06
[ Info: time: 0.100, max|u|: 1.000, CFL: 0.2000, Net mass flux: 4.6397e-05
[ Info: time: 0.200, max|u|: 1.000, CFL: 0.2000, Net mass flux: 4.1665e-05
[ Info: time: 0.300, max|u|: 1.000, CFL: 0.2000, Net mass flux: 3.7411e-05
[ Info: time: 0.400, max|u|: 1.000, CFL: 0.2000, Net mass flux: 3.3592e-05
[ Info: time: 0.500, max|u|: 1.000, CFL: 0.2000, Net mass flux: 3.0163e-05
[ Info: time: 0.600, max|u|: 1.000, CFL: 0.2000, Net mass flux: 2.7083e-05
[ Info: time: 0.700, max|u|: 1.000, CFL: 0.2000, Net mass flux: 2.4319e-05
[ Info: time: 0.800, max|u|: 1.000, CFL: 0.2000, Net mass flux: 2.1836e-05
[ Info: time: 0.900, max|u|: 1.000, CFL: 0.2000, Net mass flux: 1.9607e-05
[ Info: time: 1.000, max|u|: 1.000, CFL: 0.2000, Net mass flux: 1.7605e-05

I'm not exactly sure what the consequences of this are for simulations, but I suspect this is the culprit behind the ConjugateGradientPoissonSolver blowing up for simulations with PerturbationAdvection OBCs on ImmersedBoundaryGrids (see e.g. https://github.com/CliMA/Oceananigans.jl/issues/4451).

This is was first described in https://github.com/CliMA/Oceananigans.jl/pull/4578#issuecomment-2956351418.

I'm planning on opening a PR soon which will enforce a zero mass flux when open boundaries are used.

cc @jagoosw @glwagner

[jagoosw]

What happens if you compute the mass flux at i = 2 and i = grid.Nx? I think that it will be conserved in that region since the open boundary was implemented in a way to lead to a divergence free interior. The boundaries only get stepped at before the pressure solve so that the interior center points are divergence free (and hence the whole interior is divergance free) so I don't know how the interior could have net mass flux.

[jagoosw]

Also I ran your example as is and I can not reproduce the problem, I get zero mass flux.

Update: I checked and I had the code from https://github.com/CliMA/Oceananigans.jl/pull/4539/commits/4e2d6236782f0b02dab2e0b91a8d6604d1659fe7 so this issue was introduced by getting rid of the "secret" boundary storage. Confirm that when I run on main I do get divergence between i=1 and i=grid.Nx+1 as well as i=2 and i=grid.Nx so a recent change has broken this code.

[glwagner]

We should regard any algorithm that is asymmetric between left/right boundaries as having a bug, I think.

[tomchor]

Also I ran your example as is and I can not reproduce the problem, I get zero mass flux.

Update: I checked and I had the code from 4e2d623 so this issue was introduced by getting rid of the "secret" boundary storage. Confirm that when I run on main I do get divergence between i=1 and i=grid.Nx+1 as well as i=2 and i=grid.Nx so a recent change has broken this code.

Okay, I did some more digging. You're right that https://github.com/CliMA/Oceananigans.jl/commit/4e2d6236782f0b02dab2e0b91a8d6604d1659fe7 has zero flux. However, the code on main before that PR was merged (https://github.com/CliMA/Oceananigans.jl/commit/6c2944e0b122743147ae4e768188ebb4753c954d) was producing net mass fluxes of around 1e-3 according to the MWE above. After that PR, it went down to 1e-5. So, on main at least, that's progress.

So what I think was happening is that most of the mass flux imbalance (on main before #4578) was coming from the wrong initial condition that was happening due to the "secret storage". In the commit you linked, that IC issue was fixed and the algorithm still had a secret storage for past velocities, which appears to fix the issue completely. Later I removed the secret storage, which introduced another (smaller) mass flux imbalance of about 1e-5.

But I agree with @glwagner that, ideally, the algorithm shouldn't have any asymmetries. So ideally we'd fix whatever is causing that asymmetry to be necessary in the first place.

@jagoosw I think part of the difficulty here is that @glwagner and I still don't understand where this asymmetry is coming from. You mentioned it several times, but could you please explain more precisely what is causing this? Ideally linking the relevant code?

[glwagner]

I'm arguing that it's more than just not ideal. The asymmetry cannot be explained so may be hiding some other bug. Note that even if one empirical test is working, others are not (such as simulations with the CG solver).

[jagoosw]

I spent some time going through the code and I agree now that there is no obvious aymmetry.

The only interestying thing I've noticed is that if you set the in and outflow timescales to be the same then the error goes down to ~10^-12, and if I make it 10x faster it goes down to 10^-13. Similarly if I remove the initial pertubations it goes to zero.

[tomchor]

The only interestying thing I've noticed is that if you set the in and outflow timescales to be the same then the error goes down to ~10^-12 and if I make it 10x faster it goes down to 10^-13.

By "make it faster" do you mean setting shorter timescales? If so then yeah I noticed that too. I think it makes sense since you get closer and closer to having U₀ on all points in the boundary, thus satisfying continuity.

Similarly if I remove the initial pertubations it goes to zero.

I think that's also expected, since the velocity is U₀ everywhere and it's a trivial solution.

[jagoosw]

Going through older commits need to change the callback to:

function progress(sim)
    u, v, w = model.velocities
    cfl_value = cfl_calculator(model)
    west_mass_flux = Field(Average(view(u, 1, :, :)))
    east_mass_flux = Field(Average(view(u, grid.Nx+1, :, :)))

    compute!(west_mass_flux)
    compute!(east_mass_flux)

    net_mass_flux = east_mass_flux[] - west_mass_flux[]
    @info @sprintf("time: %.3f, max|u|: %.3f, CFL: %.4f, Net mass flux: %.4e",
                   time(sim), maximum(abs, u), cfl_value, net_mass_flux)
end

because compute didn't used to get called automatically.

But I have checked and this has always been the same.

[navidcy]

Was this resolved with #4627?

[tomchor]

Yes!