Halos should be filled as part of set!

[ali-ramadhan]

This would ensure consistency in general, as currently fields are not correct in the halo regions after set!.

This would ensure that fields are ready to be used for computed fields following set!.

cc @tomchor

[glwagner]

Posted from a slack conversation...

I’ll try to clarify what we found for @francispoulin and @tomchor, if I understand correctly what was happening. I think @tomchor was

1. Initializing a velocity field using set!(model, …)
2. Constructing a computed field ∂z_u = ComputedField(∂z(model.velocities.u))
3. Evaluating the computed field via compute!(∂z_u)
4. Plotting ∂z_u and noticing that it’s wrong on boundaries.

The reason why ∂z_u is wrong on boundaries is because the vertical derivative of u is calcaluated with a second-order difference that spans the boundary, and thus depends on “halo points” that are outside the physical domain. In Oceananigans, we use halo points to enforce Periodic, Value and Gradient boundary conditions, and will eventually use them for parallelization.

If ∂z_u were calculated during a simulation, the halo points of u would be filled. However, in this case ∂z_u was evaluated before invoking run!As a result, the halo points of u were not filled, and thus ∂z_u was wrong on boundaries.
One fix for this mentioned by @alir is to fill the halo points “manually” by calling fill_halo_regions(u, CPU()) after invoking set!. This will work provided that u does not have function boundary conditions. If the boundary conditions of u depend on other model fields or the simulation time, this will fail. The code that fills halo regions is here:

Oceananigans.jl/src/Models/IncompressibleModels/update_incompressible_model_state.jl

Lines 15 to 16 in c3b688f

	fill_halo_regions!(merge(model.velocities, model.tracers), model.architecture,
	model.clock, fields(model))

The “real” fix for this (a more generic solution that doesn’t require magical incantations in user scripts), however, is to fill_halo_regions! within set!. In fact, an even better solution is to call update_state!(model) at the end of set!. This way nonlinear diffusivities and hydrostatic pressure are also correct. We have also discussed projecting the velocity field onto an incompressible field within set!. We could make all these changes in one PR perhaps…

[glwagner]

The issue that discusses an enforce_incompressibility kwarg is #1027