[WIP] Implementation of a modified EVP rheology

[simone-silvestri]

This PR is a tentative proposal for an implementation of a modified EVP rheology and an infrastructure to solve the momentum equations of sea ice.

This PR builds on top of #24, so a lot of changes here are duplicated. All names and implementations here can be changed.
The salient points of this implementation are:

An explicit momentum solver

encoded in the ExplicitMomentumSolver type, which contains

• the rheology
• the auxiliary fields required by the rheology (similarly to how we define the diffusivity_fields in an oceananigans model)
• the number of substeps
• the ocean-sea ice drag coefficient

This type extends the solve_momentum! function where the momentum time stepping is a substepping performed with a forward-backward scheme where velocities in the Coriolis force are taken at the m+1 substep. The ocean-sea ice stress is computed semi-implicitly.
The equations solved are (on even substeps)
$$u^{n+1} = u^n + \frac{\Delta t}{\text{substeps}} \left[ \frac{C^D \rho_o}{m_i} \lVert \boldsymbol{U}_o - \boldsymbol{U}^n \rVert (u_o - u^{n+1}) + \tau_a + f \cdot v^n + g \frac{\partial \eta}{\partial x} + \frac{\partial \sigma}{\partial x_j}\right] + \frac{G^R_x}{\text{substeps}} $$
$$v^{n+1} = v^n + \frac{\Delta t}{\text{substeps}} \left[ \frac{C^D \rho_o}{m_i} \lVert \boldsymbol{U}_o - \boldsymbol{U}^n \rVert (u_o - v^{n+1}) + \tau_a - f \cdot u^{n+1} + g \frac{\partial \eta}{\partial y} + \frac{\partial \sigma}{\partial x_j}\right] + \frac{G^R_y}{\text{substeps}} $$
On odd substeps the order reverses so we solve first for $v^{n+1}$ and use it to calculate $u^{n+1}$. $G^R$ is a "rheology specific" tendency term

The substepping requires

• a compute_stresses! kernel
• ∂ⱼ_σ₁ⱼ and ∂ⱼ_σ₂ⱼ pointwise functions
• a rheology_substep pointwise function (if these are different than the substeps specified in the ExplicitMomentumSolver as in the modified EVP formulation)
• rheology_specific_forcing_x and rheology_specific_forcing_y functions for particular tendencies a particular rheology would like to add ($G^R$ in the above equations)

While the initialization requires:

• a initialize_rheology! function
• a required_auxiliary_fields function that is called with a grid as an argument and allocates any necessary auxiliary field

All these functions must be extended by a possible new rheology to be compatible with the ExplicitMomentumSolver

A modified EVP rheology

encoded in the ExplicitViscoPlasticRheology type.
It extends the compute_stresses! function in accordance with the formulation of the rheology.
The rheology_specific_forcing_* pointwise functions are the restoring to the previous velocities, for example

@inline rheology_specific_forcing_x(i, j, k, grid, r::ExplicitViscoPlasticRheology, fields, uᵢ) = fields.uⁿ[i, j, k] - uᵢ[i, j, k]

The auxiliary fields required for this rheology are:

function required_auxiliary_fields(grid, ::ExplicitViscoPlasticRheology)
    
    # TODO: What about boundary conditions?
    σ₁₁ = Field{Center, Center, Nothing}(grid)
    σ₂₂ = Field{Center, Center, Nothing}(grid)
    σ₁₂ = Field{Face, Face, Nothing}(grid)
    
    uⁿ = Field{Face, Center, Nothing}(grid)
    vⁿ = Field{Center, Face, Nothing}(grid)
    P  = Field{Center, Center, Nothing}(grid)

    substeps = Field{Face, Face, Nothing}(grid) # Dynamic substeps a la Kimmritz et al (2016)

    # An initial (safe) educated guess
    fill!(substeps, 100)

    return (; σ₁₁, σ₂₂, σ₁₂, substeps, uⁿ, vⁿ, P)
end

while the stress divergence is simply

@inline function ∂ⱼ_σ₁ⱼ(i, j, k, grid, ::ExplicitViscoPlasticRheology, fields) 
    ∂xσ₁₁ = ∂xᶠᶜᶜ(i, j, k, grid, fields.σ₁₁)
    ∂yσ₁₂ = ∂yᶠᶜᶜ(i, j, k, grid, fields.σ₁₂)

    return ∂xσ₁₁ + ∂yσ₁₂ 
end

Validation scripts

There are two validation scripts that show the interface of the dynamics: ice_advected_by_anticyclone.jl and ice_advected_on_triangular_coastline.jl. The implementation for the moment is probably not 100% correct, because the output of the first one, that is supposed to reproduce the results in Simulating Linear Kinematic Features in Viscous-Plastic Sea Ice Models yields strange results:

sea_ice_evp.mp4

Anyways, I am curious to know what people think about this implementation and open to changes, comments or even closing the PR if we decide to go in another direction