BBM rheology

docs/make.jl

@@ -17,7 +17,7 @@ example_scripts = [
     "melting_in_spring.jl",
     "freezing_of_a_lake.jl",
     "ice_advected_by_anticyclone.jl",
-    # "ice_advected_on_coastline.jl",
+    "ice_advected_on_coastline.jl",
     "arctic_basin_seasonal_cycle.jl"
 ]
 
@@ -31,7 +31,7 @@ example_pages = [
     "Melting in Spring" => "literated/melting_in_spring.md",
     "Freezing of a Lake" => "literated/freezing_of_a_lake.md",
     "Ice advected by anticyclone" => "literated/ice_advected_by_anticyclone.md",
-    # "Ice advected on coastline" => "literated/ice_advected_on_coastline.md",
+    "Ice advected on coastline" => "literated/ice_advected_on_coastline.md",
     "Arctic basin seasonal cycle" => "literated/arctic_basin_seasonal_cycle.md"
 ]
 

examples/freezing_bucket.jl

@@ -53,9 +53,9 @@ top_heat_boundary_condition = PrescribedTemperature(-10)
 # slab sea ice representation of ice_thermodynamics
 
 ice_thermodynamics = SlabSeaIceThermodynamics(grid;
-                                          internal_heat_flux,
-                                          phase_transitions,
-                                          top_heat_boundary_condition)
+                                              internal_heat_flux,
+                                              phase_transitions,
+                                              top_heat_boundary_condition)
 
 # We also prescribe a frazil ice heat flux that stops when the ice has reached a concentration of 1.
 # This heat flux represents the initial ice formation from a liquid bucket.

examples/ice_advected_by_anticyclone.jl

@@ -25,7 +25,7 @@ L  = 512kilometers
 
 # 2 km domain
 grid = RectilinearGrid(arch;
-                       size = (128, 128), 
+                       size = (256, 256), 
                           x = (0, L), 
                           y = (0, L), 
                        halo = (7, 7),
@@ -35,19 +35,23 @@ grid = RectilinearGrid(arch;
 ##### Value boundary conditions for velocities
 #####
 
-u_bcs = FieldBoundaryConditions(north=ValueBoundaryCondition(0),
-                                south=ValueBoundaryCondition(0))
+u_bcs = FieldBoundaryConditions(north=OpenBoundaryCondition(0),
+                                south=OpenBoundaryCondition(0),
+                                west=OpenBoundaryCondition(0),
+                                east=OpenBoundaryCondition(0))
 
-v_bcs = FieldBoundaryConditions(west=ValueBoundaryCondition(0),
-                                east=ValueBoundaryCondition(0))
+v_bcs = FieldBoundaryConditions(north=OpenBoundaryCondition(0),
+                                south=OpenBoundaryCondition(0),
+                                west=OpenBoundaryCondition(0),
+                                east=OpenBoundaryCondition(0))
 
 #####
 ##### Ocean sea-ice stress
 #####
 
 # Constant ocean velocities corresponding to a cyclonic eddy
-Uₒ = XFaceField(grid)
-Vₒ = YFaceField(grid)
+Uₒ = Field{Face, Face, Nothing}(grid)
+Vₒ = Field{Face, Face, Nothing}(grid)
 
 set!(Uₒ, (x, y) -> 𝓋ₒ * (2y - L) / L)
 set!(Vₒ, (x, y) -> 𝓋ₒ * (L - 2x) / L)
@@ -59,8 +63,8 @@ fill_halo_regions!((Uₒ, Vₒ))
 #### Atmosphere - sea ice stress 
 ####
 
-Uₐ = XFaceField(grid)
-Vₐ = YFaceField(grid)
+Uₐ = Field{Face, Face, Nothing}(grid)
+Vₐ = Field{Face, Face, Nothing}(grid)
 
 # Atmospheric velocities corresponding to an anticyclonic eddy moving north-east
 @inline center(t) = 256kilometers + 51.2kilometers * t / 86400
@@ -89,13 +93,18 @@ fill_halo_regions!((Uₐ, Vₐ))
 dynamics = SeaIceMomentumEquation(grid; 
                                   top_momentum_stress = (u=τₐu, v=τₐv),
                                   bottom_momentum_stress = τₒ,
-                                  coriolis     = FPlane(f=1e-4))
+                                  coriolis = FPlane(f=1e-4),
+                                  free_drift = StressBalanceFreeDrift((u=τₐu, v=τₐv), τₒ),
+                                  rheology = BrittleBinghamMaxwellRheology(),
+                                  solver   = SplitExplicitSolver(substeps=150))
 
+# Define the model!
 model = SeaIceModel(grid; 
                     dynamics,
                     ice_thermodynamics = nothing, # No ice_thermodynamics here
                     advection = WENO(order=7),
-                    boundary_conditions = (u=u_bcs, v=v_bcs))
+                    boundary_conditions = (u=u_bcs, v=v_bcs),
+                    tracers = :d)
 
 # We start with a concentration of ℵ = 1 and an 
 # initial height field with perturbations around 0.3 m
@@ -132,8 +141,9 @@ simulation.callbacks[:top_stress] = Callback(compute_wind_stress, IterationInter
 h = model.ice_thickness
 ℵ = model.ice_concentration
 u, v = model.velocities
+d = model.tracers.d
 
-outputs = (; h, u, v, ℵ)
+outputs = (; h, u, v, ℵ, d)
 
 simulation.output_writers[:sea_ice] = JLD2Writer(model, outputs;
                                                  filename = "sea_ice_advected_by_anticyclone.jld2", 
@@ -172,13 +182,14 @@ htimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "h")
 utimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "u")
 vtimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "v")
 ℵtimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "ℵ")
+dtimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "d")
 
 # Visualize!
 Nt = length(htimeseries)
 iter = Observable(1)
 
 hi = @lift(htimeseries[$iter])
-ℵi = @lift(ℵtimeseries[$iter])
+di = @lift(dtimeseries[$iter])
 ui = @lift(utimeseries[$iter])
 vi = @lift(vtimeseries[$iter])
 
@@ -187,7 +198,7 @@ ax = Axis(fig[1, 1], title = "sea ice thickness")
 heatmap!(ax, hi, colormap = :magma, colorrange = (0.23, 0.37))
 
 ax = Axis(fig[1, 2], title = "sea ice concentration")
-heatmap!(ax, ℵi, colormap = Reverse(:deep), colorrange = (0.9, 1))
+heatmap!(ax, di, colormap = Reverse(:deep), colorrange = (0.9, 1))
 
 ax = Axis(fig[2, 1], title = "zonal velocity")
 heatmap!(ax, ui, colorrange = (-0.1, 0.1))

examples/ice_advected_on_coastline.jl

@@ -61,18 +61,14 @@ Oceananigans.BoundaryConditions.fill_halo_regions!(τᵤ)
 dynamics = SeaIceMomentumEquation(grid; 
                                   top_momentum_stress = (u=τᵤ, v=τᵥ),
                                   bottom_momentum_stress = τₒ, 
-                                  rheology = ElastoViscoPlasticRheology(),
+                                  free_drift = StressBalanceFreeDrift((u=τᵤ, v=τᵥ), τₒ),
+                                  rheology = BrittleBinghamMaxwellRheology(),
                                   solver = SplitExplicitSolver(substeps=150))
 
-u_bcs = FieldBoundaryConditions(top = nothing, bottom = nothing,
-                                north = ValueBoundaryCondition(0),
-                                south = ValueBoundaryCondition(0))
-
 #Define the model! 
 model = SeaIceModel(grid; 
                     advection = WENO(order=7),
                     dynamics = dynamics,
-                    boundary_conditions = (; u=u_bcs),
                     ice_thermodynamics = nothing)
 
 # We start with a concentration of ℵ = 1 everywhere
@@ -84,7 +80,7 @@ set!(model, ℵ = 1)
 #####
 
 # run the model for 10 day
-simulation = Simulation(model, Δt = 2minutes, stop_time=2days) 
+simulation = Simulation(model, Δt=2minutes, stop_time=2days) 
 
 # Container to hold the data
 htimeseries = []
@@ -159,4 +155,4 @@ CairoMakie.record(fig, "sea_ice_advected_on_coastline.mp4", 1:Nt, framerate = 8)
 end
 nothing #hide
 
-# ![](sea_ice_advected_on_coastline.mp4)
+# ![](sea_ice_advected_on_coastline.mp4)

src/Rheologies/Rheologies.jl

@@ -1,25 +1,29 @@
 module Rheologies
 
-export ViscousRheology, ElastoViscoPlasticRheology
-export ∂ⱼ_σ₁ⱼ, ∂ⱼ_σ₂ⱼ, required_auxiliary_fields
+export ViscousRheology, ElastoViscoPlasticRheology, BrittleBinghamMaxwellRheology
+export ∂ⱼ_σ₁ⱼ, ∂ⱼ_σ₂ⱼ, rheology_auxiliary_fields
 
 using Oceananigans
 using Oceananigans.Operators
 using ClimaSeaIce: ice_mass
 
 # Nothing rheology
-required_auxiliary_fields(rheology, grid) = NamedTuple()
 initialize_rheology!(model, rheology) = nothing
-compute_stresses!(model, dynamics, rheology, Δt) = nothing
+compute_stresses!(model, dynamics, rheology, Δt, Ns) = nothing
 
 # Nothing rheology or viscous rheology
-@inline compute_substep_Δtᶠᶜᶜ(i, j, grid, Δt, rheology, substeps, fields) = Δt / substeps
-@inline compute_substep_Δtᶜᶠᶜ(i, j, grid, Δt, rheology, substeps, fields) = Δt / substeps
+@inline compute_substep_Δtᶠᶠᶜ(i, j, grid, Δt, rheology, substeps, fields) = Δt / substeps
 
 # Fallback
 @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)
 
+# No needed extra tracers
+rheology_prognostic_tracers(rheology) = ()
+
+# No needed auxiliary fields
+rheology_auxiliary_fields(rheology, grid) = NamedTuple()
+
 @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)
@@ -28,5 +32,6 @@ compute_stresses!(model, dynamics, rheology, Δt) = nothing
 include("ice_stress_divergence.jl")
 include("viscous_rheology.jl")
 include("elasto_visco_plastic_rheology.jl")
+include("brittle_bingham_maxwell_rheology.jl")
 
 end