(0.96.33) Correct diffusion with a non-linear free surface

[simone-silvestri]

PR #4546 changed the implementation of zstar and a nonlinear free surface. In particular, before #4546 we were

• 1. updating σT
• 2. performing implicit diffusion on σT
• 3. advancing σ
• 4. diagnosing T at the new timestep

PR#4546 changed this in

• 1. updating σ
• 2. computing directly T from updating σT and dividing by the new σ
• 3. implicit diffusion

After #4546, the diffusion is performed on T and not σT, but #4546 did not update the vertically implicit diffusion solver (I also think there was a bug before), for this reason, the non linear free surface solution is way more diffusive now than the same solution with a static grid. This leads to SST increasing drastically using zstar compared to using a static vertical coordinate in ClimaOcean.

This PR corrects the formulation and adds a very simple test, where if the free surface does not change, the results should be the same for static and moving grid (H+η)_moving == H_static.
If we run same test I added in this PR on main:

z_static = [i + rand() for i in -15:0]
z_static[1] = -15
z_static[end] = 0
z_moving = MutableVerticalDiscretization(z_static ./ 1.5)
    
c₀ = rand(15)

grid_static = RectilinearGrid(size=15, z=z_static, topology=(Flat, Flat, Bounded))
grid_moving = RectilinearGrid(size=15, z=z_moving, topology=(Flat, Flat, Bounded))

fill!(grid_moving.z.ηⁿ,   5)
fill!(grid_moving.z.σᶜᶜ⁻, 1.5)
fill!(grid_moving.z.σᶜᶜⁿ, 1.5)
fill!(grid_moving.z.σᶜᶠⁿ, 1.5)
fill!(grid_moving.z.σᶠᶠⁿ, 1.5)
fill!(grid_moving.z.σᶠᶜⁿ, 1.5)

model_static = HydrostaticFreeSurfaceModel(; grid = grid_static, 
                                            tracers = :c,
                                            closure = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(), κ=0.1))
                                                
model_moving = HydrostaticFreeSurfaceModel(; grid = grid_moving, 
                                            tracers = :c,
                                            closure = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(), κ=0.1))
                                        
set!(model_static, c = c₀)
set!(model_moving, c = c₀)

for _ in 1:1000
    time_step!(model_static, 1.0)
    time_step!(model_moving, 1.0)
end

lines(interior(model_static.tracers.c, 1, 1, :), label="Static Grid")
lines!(interior(model_moving.tracers.c, 1, 1, :), label="Moving Grid", linestyle=:dash)
axislegend()

we get

on main

on this PR

The detailed formulation of the equations we solve is begin added in #4588

[francispoulin]

Very nice PR @simone-silvestri . Why is it that the curve changes from increasing to decreasing on this PR? I thought that the static grid before was correct, but now we are getting different tendencies.

Also, are the equations that describe this written down anywhere?

[simone-silvestri]

Sorry, it's a consequence of the rand() at the beginning. I am writing down the equations in #4588, I ll start writing everything down tomorrow, with the objective of finishing this week.

[francispoulin]

Thanks @simone-silvestri !

[navidcy]

Clarifying question: what is σ and σT? (or are these going to be explained in #4588; shall I look there?)

[glwagner]

is there any way we can get away from this *_params design? It doesn't scale because we have to add *_params for every single addition to the model, getting messy and requiring a lot of remembering the names and places of functions

[simone-silvestri]

You prefer having a type that controls the parameters rather than a function?

[glwagner]

👍