Improve EVP rheology

src/Rheologies/elasto_visco_plastic_rheology.jl

@@ -12,23 +12,32 @@ using KernelAbstractions: @kernel, @index
 # Where:
 # σᵢⱼ(u) = 2η ϵ̇ᵢⱼ + [(ζ - η) * (ϵ̇₁₁ + ϵ̇₂₂) - P / 2] δᵢⱼ
 #
-struct ElastoViscoPlasticRheology{FT}
+struct ElastoViscoPlasticRheology{FT, IP}
     ice_compressive_strength :: FT # compressive strength
     ice_compaction_hardening :: FT # compaction hardening
     yield_curve_eccentricity :: FT # elliptic yield curve eccentricity
     minimum_plastic_stress :: FT # minimum plastic parameter (transitions to viscous behaviour)
     min_relaxation_parameter :: FT # minimum number of substeps expressed as the dynamic coefficient
     max_relaxation_parameter :: FT # maximum number of substeps expressed as the dynamic coefficient
+    relaxation_strength :: FT # strength of the relaxation parameter
+    pressure_formulation :: IP # formulation of ice pressure
+    ElastoViscoPlasticRheology(P::FT, C::FT, e::FT, Δ_min::FT, α⁻::FT, α⁺::FT, c::FT, ip::IP) where {FT, IP}  = 
+        new{FT, IP}(P, C, e, Δ_min, α⁻, α⁺, c, ip)
 end
 
+struct ReplacementPressure end
+struct IceStrength end
+
 """
     ElastoViscoPlasticRheology(FT::DataType = Float64; 
                                ice_compressive_strength = 27500, 
                                ice_compaction_hardening = 20, 
                                yield_curve_eccentricity = 2, 
                                minimum_plastic_stress = 2e-9,
                                min_relaxation_parameter = 50,
-                               max_relaxation_parameter = 300)
+                               max_relaxation_parameter = 300,
+                               relaxation_strength = π^2,
+                               pressure_formulation = ReplacementPressure())
 
 Constructs an `ElastoViscoPlasticRheology` object representing a "modified" elasto-visco-plastic
 rheology for slab sea ice dynamics that follows the implementation of Kimmritz et al (2016).
@@ -44,8 +53,9 @@ parameterized as ``P★ h exp( - C ⋅ ( 1 - ℵ ))`` where ``P★`` is the `ice
 The stresses are substepped using a dynamic substepping coefficient ``α`` that is
 spatially varying and computed dynamically as in Kimmritz et al (2016)
 In particular: α = sqrt(γ²) 
-where γ² = ζ * π² * (Δt / mᵢ) / Az is a stability parameter with ``Az`` is the area of the grid cell, 
-``mᵢ`` the ice mass, and ``Δt`` the time step.
+where γ² = ζ * cα * (Δt / mᵢ) / Az is a stability parameter with ``Az`` is the area of the grid cell, 
+``mᵢ`` the ice mass, ``Δt`` the time step, and ``cα`` a numerical stability parameter which controls the 
+stregth of ``γ²``.
 
 The stresses are substepped with:
 ```math
@@ -65,28 +75,34 @@ Arguments
 Keyword Arguments
 =================
 
-- `ice_compressive_strength`: parameter expressing compressive strength (in Nm²). Default `27500`.
-- `ice_compaction_hardening`: exponent coefficient for compaction hardening. Default `20`.
-- `yield_curve_eccentricity`: eccentricity of the elliptic yield curve. Default `2`.
+- `ice_compressive_strength`: parameter expressing compressive strength (in Nm²). Default: `27500`.
+- `ice_compaction_hardening`: exponent coefficient for compaction hardening. Default: `20`.
+- `yield_curve_eccentricity`: eccentricity of the elliptic yield curve. Default: `2`.
 - `Δ_min`: Minimum value for the visco-plastic parameter. Limits the maximum viscosity of the ice, 
-           transitioning the ice from a plastic to a viscous behaviour. Default value is `1e-10`.
-- `min_relaxation_parameter`: Minimum value for the relaxation parameter `α`. Default value is `30`.
-- `max_relaxation_parameter`: Maximum value for the relaxation parameter `α`. Default value is `500`.
+           transitioning the ice from a plastic to a viscous behaviour. Default: `1e-10`.
+- `min_relaxation_parameter`: Minimum value for the relaxation parameter `α`. Default: `30`.
+- `max_relaxation_parameter`: Maximum value for the relaxation parameter `α`. Default: `500`.
+- `relaxation_strength`: parameter controlling the strength of the relaxation parameter. The maximum value is `π²`, see Kimmritz et al (2016). Default: `π² / 2`.
+- `pressure_formulation`: can use `ReplacementPressure` or `IceStrength`. The replacement pressure formulation avoids ice motion in the absence of forcing. Default: `ReplacementPressure`.
 """
 function ElastoViscoPlasticRheology(FT::DataType = Float64; 
                                     ice_compressive_strength = 27500, 
                                     ice_compaction_hardening = 20, 
                                     yield_curve_eccentricity = 2, 
                                     minimum_plastic_stress = 2e-9,
                                     min_relaxation_parameter = 50,
-                                    max_relaxation_parameter = 300)
+                                    max_relaxation_parameter = 300,
+                                    relaxation_strength = π^2,
+                                    pressure_formulation = ReplacementPressure())
 
     return ElastoViscoPlasticRheology(convert(FT, ice_compressive_strength), 
                                       convert(FT, ice_compaction_hardening), 
                                       convert(FT, yield_curve_eccentricity),
                                       convert(FT, minimum_plastic_stress),
                                       convert(FT, min_relaxation_parameter),
-                                      convert(FT, max_relaxation_parameter))
+                                      convert(FT, max_relaxation_parameter),
+                                      convert(FT, relaxation_strength),
+                                      pressure_formulation)
 end
 
 function required_auxiliary_fields(r::ElastoViscoPlasticRheology, grid)
@@ -116,7 +132,9 @@ Adapt.adapt_structure(to, r::ElastoViscoPlasticRheology) =
                                Adapt.adapt(to, r.yield_curve_eccentricity),
                                Adapt.adapt(to, r.minimum_plastic_stress),
                                Adapt.adapt(to, r.min_relaxation_parameter),
-                               Adapt.adapt(to, r.max_relaxation_parameter))
+                               Adapt.adapt(to, r.max_relaxation_parameter),
+                               Adapt.adapt(to, r.relaxation_strength),
+                               Adapt.adapt(to, r.pressure_formulation))
 
 """
     initialize_rheology!(model, rheology::ElastoViscoPlasticRheology)
@@ -212,6 +230,15 @@ end
     @inbounds fields.Δ[i, j, 1] = Δᶜᶜᶜ
 end
 
+@inline ice_pressure(i, j, k, grid, ::IceStrength, r, fields) = @inbounds fields.P[i, j, k]
+
+@inline function ice_pressure(i, j, k, grid, ::ReplacementPressure, r, fields)
+    Pᶜᶜᶜ = @inbounds fields.P[i, j, k]
+    Δᶜᶜᶜ = @inbounds fields.Δ[i, j, k]
+    Δm   = r.minimum_plastic_stress
+    return Pᶜᶜᶜ * Δᶜᶜᶜ / (Δᶜᶜᶜ + Δm)
+end
+
 # Compute the visco-plastic stresses for a slab sea ice model.
 # The function updates the internal stress variables `σ₁₁`, `σ₂₂`, and `σ₁₂` in the `rheology` object
 # following the αEVP formulation of Kimmritz et al (2016).
@@ -220,9 +247,10 @@ end
     kᴺ   = size(grid, 3)
 
     e⁻² = rheology.yield_curve_eccentricity^(-2)
-    Δm  = rheology.minimum_plastic_stress
     α⁺  = rheology.max_relaxation_parameter
     α⁻  = rheology.min_relaxation_parameter
+    cα  = rheology.relaxation_strength
+    ip  = rheology.pressure_formulation
 
     σ₁₁ = fields.σ₁₁
     σ₂₂ = fields.σ₂₂
@@ -234,13 +262,11 @@ end
     ϵ̇₂₂ = strain_rate_yy(i, j, kᴺ, grid, u, v) 
     ϵ̇₁₂ = strain_rate_xy(i, j, kᴺ, grid, u, v)
 
-    Pᶜᶜᶜ = @inbounds fields.P[i, j, 1]
     ζᶜᶜᶜ = @inbounds fields.ζ[i, j, 1]
-    Δᶜᶜᶜ = @inbounds fields.Δ[i, j, 1]
     ζᶠᶠᶜ = ℑxyᶠᶠᵃ(i, j, 1, grid, fields.ζ)
 
     # replacement pressure?
-    Pᵣ = Pᶜᶜᶜ * Δᶜᶜᶜ / (Δᶜᶜᶜ + Δm)
+    Pᵣ = ice_pressure(i, j, 1, grid, ip, rheology, fields)
 
     ηᶜᶜᶜ = ζᶜᶜᶜ * e⁻²
     ηᶠᶠᶜ = ζᶠᶠᶜ * e⁻²
@@ -256,11 +282,11 @@ end
 
     # Update coefficients for substepping using dynamic substepping
     # with spatially varying coefficients as in Kimmritz et al (2016)
-    γ²ᶜᶜᶜ = ζᶜᶜᶜ * π^2 * Δt / mᵢᶜᶜᶜ / Azᶜᶜᶜ(i, j, 1, grid)
+    γ²ᶜᶜᶜ = ζᶜᶜᶜ * cα * Δt / mᵢᶜᶜᶜ / Azᶜᶜᶜ(i, j, 1, grid)
     γ²ᶜᶜᶜ = ifelse(isnan(γ²ᶜᶜᶜ), α⁺^2, γ²ᶜᶜᶜ) # In case both ζᶜᶜᶜ and mᵢᶜᶜᶜ are zero
     γᶜᶜᶜ  = clamp(sqrt(γ²ᶜᶜᶜ), α⁻, α⁺)
 
-    γ²ᶠᶠᶜ = ζᶠᶠᶜ * π^2 * Δt / mᵢᶠᶠᶜ / Azᶠᶠᶜ(i, j, 1, grid)
+    γ²ᶠᶠᶜ = ζᶠᶠᶜ * cα * Δt / mᵢᶠᶠᶜ / Azᶠᶠᶜ(i, j, 1, grid)
     γ²ᶠᶠᶜ = ifelse(isnan(γ²ᶠᶠᶜ), α⁺^2, γ²ᶠᶠᶜ) # In case both ζᶠᶠᶜ and mᵢᶠᶠᶜ are zero
     γᶠᶠᶜ  = clamp(sqrt(γ²ᶠᶠᶜ), α⁻, α⁺)
 
@@ -306,4 +332,4 @@ end
     user_forcing = v_forcing(i, j, k, grid, fields)
     rheology_forcing = @inbounds (fields.vⁿ[i, j, k] - fields.v[i, j, k]) / Δt / ℑyᵃᶠᵃ(i, j, k, grid, fields.α)
     return user_forcing + rheology_forcing
-end
+end