set! should enforce incompressibility when setting an initial condition for an Incompressible model

[navidcy]

For example, here:

Oceananigans.jl/examples/two_dimensional_turbulence.jl

Line 38 in 6b2cea3

set!(model, u=u₀, v=u₀)

clearly the initial condition is not incompressible.

set! should somehow remove the compressible part of the initial conditions and perhaps also print out a warning so that the user knows that the model is not initialized exactly with the fields they provided?

[glwagner]

There is a comment about divergent initial conditions here:

https://clima.github.io/OceananigansDocumentation/stable/model_setup/setting_initial_conditions/

A related issue is #951 .

We could somehow change the behavior of set! to project the initial condition onto an incompressible field. For example, we might use the pressure solver. While it would be surprising that the initial condition never ends up exactly as its set by the user, it is similarly surprising to observe a large change in the velocity field after a single time step. And set! is intended to be used prior to time stepping.

I'm not sure I agree about printing a warning based on some metric diagnosed from the velocity field passed to set!. However, we could print out an info message that gives the user some information about how much their velocity field was changed during the initialization process. This is related to #1013 .

[ali-ramadhan]

We could use a kwarg like set!(model, u=u₀, v=u₀, enforce_incompressibility=true) so it's optional based on what the user wants.

[glwagner]

How exactly would we enforce incompressibility? Should we perform a pressure correction step? Something like

function set!(model; enforce_incompressibility=false, kwargs...

    # everything else in `set!`

    if enforce_incompressibility
        calculate_pressure_correction!(model, 1.0)
        pressure_correct_velocities!(model, 1.0)
    end

    return nothing
end

(This syntax will work once #1057 is merged). The "1.0" is a pseudo-timestep that is technically irrelevant (but can't be Inf or 0). The two above functions come from the time-stepping routine:

Oceananigans.jl/src/TimeSteppers/runge_kutta_3.jl

Lines 73 to 74 in c35af73

	calculate_pressure_correction!(model, first_stage_Δt)
	pressure_correct_velocities!(model, first_stage_Δt)

I suppose this projects the user-defined velocity field onto an incompressible field, which is more or less what we want? It seems better than recomputing w from continuity.

[ali-ramadhan]

Yeah that seems like a nice solution! Especially as it relies on existing functionality. To recompute w from continuity we'd have to add back in some kernels.