Remove the use of indices for storage in PerturbationAdvection

[tomchor]

Based on a conversation with @glwagner and @xkykai, this PR removes the use of indices for storage in PerturbationAdvection and makes a couple of fixes and improvements. On main, this index sets up a hidden/secret storage for left open boundary conditions at i=0 that saves the velocity of the previous timestep. This is not done for the right boundary, so this PR will make both left and right PerturbationAdvectionOBCs the same.

@jagoosw I'm curious to hear your thoughts here. We went through the code and couldn't find the need for a storage for the velocity at a previous timestep. In fact, if I simply get rid of it (as this PR does) the results basically don't change, with the exception of the initial velocity. Here are two identical simulations, the only difference being the use of a stored velocity by one and not the other:

_cylinder_xy_PerturbationAdvection_Nx_40.jld2_old.mp4

_cylinder_xy_PerturbationAdvection_Nx_40.jld2_new.mp4

[glwagner]

Regarding the initial velocity --- is the result on this PR correct vs the previous result wrong?

[tomchor]

Regarding the initial velocity --- is the result on this PR correct vs the previous result wrong?

The initial condition for the velocity is u=1, so yes, the previous IC is wrong.

The reason is that at t=0 the storage index (i=0) hasn't been populated yet, therefore that index contains all zeros. So, in the first iteration, do the PerturbationAdvection matching scheme ends up setting the velocity at the wall (i=1) to zero, which causes the velocity in the interior to be zero (presumably because of the pressure solver?).

[glwagner]

Regarding the initial velocity --- is the result on this PR correct vs the previous result wrong?

The initial condition for the velocity is u=1, so yes, the previous IC is wrong.

The reason is that at t=0 the storage index (i=0) hasn't been populated yet, therefore that index contains all zeros. So, in the first iteration, do the PerturbationAdvection matching scheme ends up setting the velocity at the wall (i=1) to zero, which causes the velocity in the interior to be zero (presumably because of the pressure solver?).

I guess if we want to keep the former scheme, we also need to fix this issue. The storage index should be computed in the model constructor, shouldn't it?

[tomchor]

For the record here, the issue with the initial condition in the original version of PerturbationAdvection creates a nonzero net mass flux at early times. 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)

On main, this script produces:

[ Info: time: 0.000, max|u|: 1.000, CFL: 0.1950, Net mass flux: 1.0001e+00
[ Info: time: 0.100, max|u|: 1.000, CFL: 0.2000, Net mass flux: -1.1620e-02
[ Info: time: 0.200, max|u|: 1.000, CFL: 0.2000, Net mass flux: -1.0537e-02
[ Info: time: 0.300, max|u|: 1.000, CFL: 0.2000, Net mass flux: -9.4626e-03
[ Info: time: 0.400, max|u|: 1.000, CFL: 0.2000, Net mass flux: -8.4977e-03
[ Info: time: 0.500, max|u|: 1.000, CFL: 0.2000, Net mass flux: -7.6311e-03
[ Info: time: 0.600, max|u|: 1.000, CFL: 0.2000, Net mass flux: -6.8528e-03
[ Info: time: 0.700, max|u|: 1.000, CFL: 0.2000, Net mass flux: -6.1538e-03
[ Info: time: 0.800, max|u|: 1.000, CFL: 0.2000, Net mass flux: -5.5261e-03
[ Info: time: 0.900, max|u|: 1.000, CFL: 0.2000, Net mass flux: -4.9623e-03
[ Info: time: 1.000, max|u|: 1.000, CFL: 0.2000, Net mass flux: -4.4560e-03

Indicating some non-negligible net mass flux. On this branch, 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

So the two things to note here are that:

1. Properly satisfying the initial condition is important, so even if we don't remove the storage, this needs to be dealt with
2. Even with the proper initial condition, there is a net mass flux, which we probably need to fix and which might be related to the difficulties in getting PerturbationAdvection to work with immersed boundaries. (The net mass flux is ~100 times smaller with the fix in this PR though, so I'm inclined to move forward with this PR.)

@jagoosw let me know if you have anything to add here, otherwise I'll move forward and merge this PR.