Stress balance free drift model

docs/src/library/public.md

@@ -39,9 +39,9 @@ Modules = [ClimaSeaIce.Rheologies]
 Private = false
 ```
 
-## ClimaSeaIce.SeaIceMomentumEquations
+## ClimaSeaIce.SeaIceDynamics
 
 ```@autodocs
-Modules = [ClimaSeaIce.SeaIceMomentumEquations]
+Modules = [ClimaSeaIce.SeaIceDynamics]
 Private = false
 ```

examples/ice_advected_by_anticyclone.jl

@@ -4,7 +4,7 @@
 #
 #
 using ClimaSeaIce
-using ClimaSeaIce.SeaIceMomentumEquations
+using ClimaSeaIce.SeaIceDynamics
 using ClimaSeaIce.Rheologies
 using Oceananigans
 using Oceananigans.Units
@@ -86,18 +86,13 @@ fill_halo_regions!((Uₐ, Vₐ))
 # We use an elasto-visco-plastic rheology and WENO seventh order 
 # for advection of h and ℵ
 
-momentum_equations = SeaIceMomentumEquation(grid; 
-                                            top_momentum_stress = (u=τₐu, v=τₐv),
-                                            bottom_momentum_stress = τₒ,
-                                            coriolis = FPlane(f=1e-4),
-                                            ocean_velocities = (u = Uₒ, v = Vₒ),
-                                            rheology = ElastoViscoPlasticRheology(),
-                                            solver   = SplitExplicitSolver(substeps=150))
-
-# Define the model!
+dynamics = SeaIceMomentumEquation(grid; 
+                                  top_momentum_stress = (u=τₐu, v=τₐv),
+                                  bottom_momentum_stress = τₒ,
+                                  coriolis     = FPlane(f=1e-4))
 
 model = SeaIceModel(grid; 
-                    dynamics = momentum_equations,
+                    dynamics,
                     ice_thermodynamics = nothing, # No ice_thermodynamics here
                     advection = WENO(order=7),
                     boundary_conditions = (u=u_bcs, v=v_bcs))

examples/ice_advected_on_coastline.jl

@@ -1,5 +1,5 @@
 using ClimaSeaIce
-using ClimaSeaIce.SeaIceMomentumEquations
+using ClimaSeaIce.SeaIceDynamics
 using ClimaSeaIce.Rheologies
 using Oceananigans
 using Oceananigans.Units

src/ClimaSeaIce.jl

@@ -39,15 +39,15 @@ include("forward_euler_timestepper.jl")
 
 include("SeaIceThermodynamics/SeaIceThermodynamics.jl")
 include("Rheologies/Rheologies.jl")
-include("SeaIceMomentumEquations/SeaIceMomentumEquations.jl")
+include("SeaIceDynamics/SeaIceDynamics.jl")
 include("sea_ice_model.jl")
 include("sea_ice_advection.jl")
 include("tracer_tendency_kernel_functions.jl")
 include("sea_ice_time_stepping.jl")
 include("EnthalpyMethodSeaIceModel.jl")
 
 using .SeaIceThermodynamics
-using .SeaIceMomentumEquations
+using .SeaIceDynamics
 using .Rheologies
 
 # Advection timescale for a `SeaIceModel`. Sea ice dynamics are two-dimensional so 

src/Rheologies/Rheologies.jl

@@ -20,6 +20,11 @@ compute_stresses!(model, dynamics, rheology, Δt) = nothing
 @inline sum_of_forcing_u(i, j, k, grid, rheology, u_forcing, fields, Δt) = u_forcing(i, j, k, grid, fields)
 @inline sum_of_forcing_v(i, j, k, grid, rheology, v_forcing, fields, Δt) = v_forcing(i, j, k, grid, fields)
 
+@inline ice_stress_ux(i, j, k, grid, ::Nothing, args...) = zero(grid)
+@inline ice_stress_uy(i, j, k, grid, ::Nothing, args...) = zero(grid)
+@inline ice_stress_vx(i, j, k, grid, ::Nothing, args...) = zero(grid)
+@inline ice_stress_vy(i, j, k, grid, ::Nothing, args...) = zero(grid)
+
 include("ice_stress_divergence.jl")
 include("viscous_rheology.jl")
 include("elasto_visco_plastic_rheology.jl")

src/SeaIceMomentumEquations/SeaIceMomentumEquations.jl → src/SeaIceDynamics/SeaIceDynamics.jl

@@ -1,8 +1,8 @@
-module SeaIceMomentumEquations
+module SeaIceDynamics
 
 # The only functions provided by the module
 export compute_momentum_tendencies!, time_step_momentum!
-export SeaIceMomentumEquation, ExplicitSolver, SplitExplicitSolver, SemiImplicitStress
+export SeaIceMomentumEquation, ExplicitSolver, SplitExplicitSolver, SemiImplicitStress, StressBalanceFreeDrift
 
 using ClimaSeaIce
 
@@ -46,8 +46,9 @@ import Oceananigans: fields
 time_step_momentum!(model, dynamics, Δt) = nothing
 compute_momentum_tendencies!(model, dynamics, Δt) = nothing
 
-include("sea_ice_momentum_equations.jl")
 include("sea_ice_external_stress.jl")
+include("stress_balance_free_drift.jl")
+include("sea_ice_momentum_equations.jl")
 include("momentum_tendencies_kernel_functions.jl")
 include("explicit_momentum_equations.jl")
 include("split_explicit_momentum_equations.jl")

src/SeaIceMomentumEquations/explicit_momentum_equations.jl → src/SeaIceDynamics/explicit_momentum_equations.jl

@@ -18,7 +18,7 @@ function time_step_momentum!(model, ::ExplicitMomentumEquation, Δt)
                          ρ = model.ice_density))
 
 
-    ocean_velocities = dynamics.ocean_velocities
+    free_drift = dynamics.free_drift
     clock = model.clock
     minimum_mass = dynamics.minimum_mass
     minimum_concentration = dynamics.minimum_concentration
@@ -27,15 +27,15 @@ function time_step_momentum!(model, ::ExplicitMomentumEquation, Δt)
     bottom_stress = dynamics.external_momentum_stresses.bottom
 
     launch!(arch, grid, :xy, _step_velocities!, u, v, grid, Gⁿ, Δt, 
-            top_stress, bottom_stress, ocean_velocities, 
+            top_stress, bottom_stress, free_drift, 
             minimum_mass, minimum_concentration, clock, model_fields)
 
     return nothing
 end
 
 @kernel function _step_velocities!(u, v, grid, Gⁿ, Δt, 
                                    top_stress, bottom_stress, 
-                                   ocean_velocities, minimum_mass, minimum_concentration, clock, fields)
+                                   free_drift, minimum_mass, minimum_concentration, clock, fields)
 
     i, j = @index(Global, NTuple)
     kᴺ   = size(grid, 3)
@@ -56,8 +56,8 @@ end
         uᴰ = (u[i, j, 1] + Δt * Gⁿ.u[i, j, 1]) / (1 + Δt * τuᵢ)
         vᴰ = (v[i, j, 1] + Δt * Gⁿ.v[i, j, 1]) / (1 + Δt * τvᵢ)
 
-        uᶠ = free_drift_u(i, j, kᴺ, grid, ocean_velocities)
-        vᶠ = free_drift_v(i, j, kᴺ, grid, ocean_velocities)
+        uᶠ = free_drift_u(i, j, kᴺ, grid, free_drift, clock, fields)
+        vᶠ = free_drift_v(i, j, kᴺ, grid, free_drift, clock, fields)
 
         sea_ice = (mᶠᶜ ≥ minimum_mass) & (ℵᶠᶜ ≥ minimum_concentration)
         u[i, j, 1] = ifelse(sea_ice, uᴰ, uᶠ)

src/SeaIceMomentumEquations/sea_ice_momentum_equations.jl → src/SeaIceDynamics/sea_ice_momentum_equations.jl

@@ -1,12 +1,12 @@
 using ClimaSeaIce.Rheologies
 using Adapt
 
-struct SeaIceMomentumEquation{S, C, R, V, A, ES, FT}
+struct SeaIceMomentumEquation{S, C, R, F, A, ES, FT}
     coriolis :: C
     rheology :: R
     auxiliary_fields :: A
     solver :: S
-    ocean_velocities :: V
+    free_drift :: F
     external_momentum_stresses :: ES
     minimum_concentration :: FT
     minimum_mass :: FT
@@ -17,16 +17,19 @@ struct ExplicitSolver end
 
 """
     SeaIceMomentumEquation(grid; 
-                           coriolis=nothing,
-                           rheology=ElastoViscoPlasticRheology(eltype(grid)),
-                           auxiliary_fields=NamedTuple(),
-                           ocean_velocities=nothing,
-                           solver=ExplicitSolver(),
-                           minimum_concentration=1e-3,
-                           minimum_mass=1.0)
+                           coriolis = nothing,
+                           rheology = ElastoViscoPlasticRheology(eltype(grid)),
+                           auxiliary_fields = NamedTuple(),
+                           top_momentum_stress    = (u=nothing, v=nothing),
+                           bottom_momentum_stress = (u=nothing, v=nothing),
+                           free_drift = StressBalanceFreeDrift(top_momentum_stress, bottom_momentum_stress),
+                           solver = SplitExplicitSolver(150),
+                           minimum_concentration = 1e-3,
+                           minimum_mass = 1.0)
 
 Constructs a `SeaIceMomentumEquation` object that controls the dynamical evolution of sea-ice momentum.
 The sea-ice momentum obey the following evolution equation:
+
 ```math
     ∂u                   τₒ    τₐ
     -- + f x u = ∇ ⋅ σ + --  + -- 
@@ -46,8 +49,8 @@ Keyword Arguments
 - `coriolis`: Parameters for the background rotation rate of the model.
 - `rheology`: The sea ice rheology model, default is `ElastoViscoPlasticRheology(eltype(grid))`.
 - `auxiliary_fields`: A named tuple of auxiliary fields, default is an empty `NamedTuple()`.
-- `ocean_velocities`: The ocean surface velocities used to limit the sea ice momentum when the mass or the concentration are
-                      below a certain threshold. default is `nothing` (indicating that the free drift velocities are zero).
+- `free_drift`: The free drift velocities used to limit sea ice momentum when the mass or the concentration are
+                below a certain threshold. default is `nothing` (indicating that the free drift velocities are zero).
 - `solver`: The momentum solver to be used.
 - `minimum_concentration`: The minimum sea ice concentration above which the sea ice velocity is dynamically calculated, default is `1e-3`.
 - `minimum_mass`: The minimum sea ice mass per area above which the sea ice velocity is dynamically calculated, default is `1.0 kg/m²`.
@@ -56,10 +59,10 @@ function SeaIceMomentumEquation(grid;
                                 coriolis = nothing,
                                 rheology = ElastoViscoPlasticRheology(eltype(grid)),
                                 auxiliary_fields = NamedTuple(),
-                                ocean_velocities = nothing,
-                                top_momentum_stress    = (u = nothing, v = nothing),
-                                bottom_momentum_stress = (u = nothing, v = nothing),
-                                solver = ExplicitSolver(),
+                                top_momentum_stress    = nothing,
+                                bottom_momentum_stress = nothing,
+                                free_drift = StressBalanceFreeDrift(top_momentum_stress, bottom_momentum_stress),
+                                solver = SplitExplicitSolver(150),
                                 minimum_concentration = 1e-3,
                                 minimum_mass = 1.0)
 
@@ -73,18 +76,10 @@ function SeaIceMomentumEquation(grid;
                                   rheology, 
                                   auxiliary_fields, 
                                   solver,
-                                  ocean_velocities,
+                                  free_drift,
                                   external_momentum_stresses,
                                   convert(FT, minimum_concentration),
                                   convert(FT, minimum_mass))
 end
 
 fields(mom::SeaIceMomentumEquation) = mom.auxiliary_fields
-
-# Just passing ocean velocities without mitigation
-@inline free_drift_u(i, j, k, grid, f) = @inbounds f.u[i, j, k] 
-@inline free_drift_v(i, j, k, grid, f) = @inbounds f.v[i, j, k] 
-
-# Passing no velocities
-@inline free_drift_u(i, j, k, grid, ::Nothing) = zero(grid)
-@inline free_drift_v(i, j, k, grid, ::Nothing) = zero(grid)

src/SeaIceMomentumEquations/split_explicit_momentum_equations.jl → src/SeaIceDynamics/split_explicit_momentum_equations.jl

@@ -35,7 +35,7 @@ function time_step_momentum!(model, dynamics::SplitExplicitMomentumEquation, Δt
 
     u, v = model.velocities
 
-    ocean_velocities = dynamics.ocean_velocities
+    free_drift = dynamics.free_drift
     clock = model.clock
     coriolis = dynamics.coriolis
 
@@ -55,8 +55,10 @@ function time_step_momentum!(model, dynamics::SplitExplicitMomentumEquation, Δt
                          ℵ = model.ice_concentration, 
                          ρ = model.ice_density))
 
-    u_velocity_kernel!, _ = configure_kernel(arch, grid, :xy, _u_velocity_step!)
-    v_velocity_kernel!, _ = configure_kernel(arch, grid, :xy, _v_velocity_step!)
+    active_cells_map = Oceananigans.Grids.get_active_column_map(grid)
+
+    u_velocity_kernel!, _ = configure_kernel(arch, grid, :xy, _u_velocity_step!; active_cells_map)
+    v_velocity_kernel!, _ = configure_kernel(arch, grid, :xy, _v_velocity_step!; active_cells_map)
 
     substeps = dynamics.solver.substeps
 
@@ -73,24 +75,24 @@ function time_step_momentum!(model, dynamics::SplitExplicitMomentumEquation, Δt
         # In odd substeps we switch and calculate vⁿ⁺¹ = f(uⁿ) and uⁿ⁺¹ = f(vⁿ⁺¹).
         if iseven(substep) 
             u_velocity_kernel!(u, grid, Δt, substeps, rheology, model_fields, 
-                               ocean_velocities, clock, coriolis,
+                               free_drift, clock, coriolis,
                                minimum_mass, minimum_concentration, 
                                u_immersed_bc, top_stress, bottom_stress, u_forcing)
 
             v_velocity_kernel!(v, grid, Δt, substeps, rheology, model_fields, 
-                               ocean_velocities, clock, coriolis, 
+                               free_drift, clock, coriolis, 
                                minimum_mass, minimum_concentration,
                                v_immersed_bc, top_stress, bottom_stress, v_forcing)
 
         else
             v_velocity_kernel!(v, grid, Δt, substeps, rheology, model_fields, 
-                               ocean_velocities, clock, coriolis, 
+                               free_drift, clock, coriolis, 
                                minimum_mass, minimum_concentration,
                                v_immersed_bc, top_stress, bottom_stress, v_forcing)
 
 
             u_velocity_kernel!(u, grid, Δt, substeps, rheology, model_fields, 
-                               ocean_velocities, clock, coriolis,
+                               free_drift, clock, coriolis,
                                minimum_mass, minimum_concentration, 
                                u_immersed_bc, top_stress, bottom_stress, u_forcing)
         end
@@ -107,7 +109,7 @@ end
 
 @kernel function _u_velocity_step!(u, grid, Δt, 
                                    substeps, rheology, 
-                                   model_fields, ocean_velocities, 
+                                   model_fields, free_drift, 
                                    clock, coriolis, 
                                    minimum_mass, minimum_concentration,
                                    u_immersed_bc, u_top_stress, u_bottom_stress, u_forcing)
@@ -127,7 +129,7 @@ end
 
     τuᵢ = ifelse(mᵢ ≤ 0, zero(grid), τuᵢ)
     uᴰ  = @inbounds (u[i, j, 1] + Δτ * Gu) / (1 + Δτ * τuᵢ) # dynamical velocity 
-    uᶠ  = free_drift_u(i, j, kᴺ, grid, ocean_velocities) # free drift velocity
+    uᶠ  = free_drift_u(i, j, kᴺ, grid, free_drift, clock, model_fields) # free drift velocity
 
     # If the ice mass or the ice concentration are below a certain threshold, 
     # the sea ice velocity is set to the free drift velocity
@@ -138,7 +140,7 @@ end
 
 @kernel function _v_velocity_step!(v, grid, Δt, 
                                    substeps, rheology, 
-                                   model_fields, ocean_velocities, 
+                                   model_fields, free_drift, 
                                    clock, coriolis, 
                                    minimum_mass, minimum_concentration,
                                    v_immersed_bc, v_top_stress, v_bottom_stress, v_forcing)
@@ -159,7 +161,7 @@ end
     τvᵢ = ifelse(mᵢ ≤ 0, zero(grid), τvᵢ)
 
     vᴰ = @inbounds (v[i, j, 1] + Δτ * Gv) / (1 + Δτ * τvᵢ)# dynamical velocity 
-    vᶠ = free_drift_v(i, j, kᴺ, grid, ocean_velocities) # free drift velocity
+    vᶠ = free_drift_v(i, j, kᴺ, grid, free_drift, clock, model_fields)  # free drift velocity
 
     # If the ice mass or the ice concentration are below a certain threshold, 
     # the sea ice velocity is set to the free drift velocity