Implement a few ice-ocean examples

.gitignore

@@ -24,3 +24,9 @@ docs/site/
 # Manifest.toml
 *.swp
 *.jld2
+*.mp4
+*.DS_Store
+*.nc
+*.json
+*.sha512
+data

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.9.2"
 manifest_format = "2.0"
-project_hash = "c88ef8733706dadcee8b23c8aa1dbaecea03d700"
+project_hash = "eaa6fb5c139ff8d22cf5ceaa00746809eb5e5859"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -534,9 +534,11 @@ version = "1.2.0"
 
 [[deps.Oceananigans]]
 deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"]
-git-tree-sha1 = "996eb159d0a1380d09beb7190ed0807d51ba0924"
+git-tree-sha1 = "65b49b4f1773a76d06f9b42de95626dd996a8c94"
+repo-rev = "glw/better-simulation-interface"
+repo-url = "https://github.com/CliMA/Oceananigans.jl.git"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
-version = "0.87.1"
+version = "0.88.0"
 
 [[deps.OffsetArrays]]
 deps = ["Adapt"]

Project.toml

@@ -4,10 +4,12 @@ authors = ["Climate Modeling Alliance and contributors"]
 version = "0.1.0"
 
 [deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40"
 
 [compat]
 KernelAbstractions = "0.9"

codecov.yml

@@ -0,0 +1,14 @@
+coverage:
+  range: 50..90
+  round: down
+  precision: 2
+
+  status:
+    project:
+      default:
+        target: auto
+        threshold: 5%
+    patch:
+      default:
+        target: auto
+        threshold: 20%

examples/freezing_bucket.jl

@@ -1,26 +1,55 @@
+# # A freezing bucket
+#
+# A common laboratory experiment freezes an insultated bucket of water
+# from the top down, using a metal lid to keep the top of the bucket
+# at some constant, very cold temperature. In this example, we simulate such
+# a scenario using `SlabSeaIceModel`. Here, the bucket is perfectly insulated
+# and infinitely deep, like many buckets are: if the `Simulation` is run for longer,
+# the ice will keep freezing, and freezing, and will never run out of water.
+# Also, the water in the infinite bucket is (somehow) all at the same temperature,
+# in equilibrium with the ice-water interface (and therefore fixed at the melting
+# temperature). Yep, this kind of thing happens _all the time_.
+#
+# We start by `using Oceananigans` to bring in functions for building grids
+# and `Simulation`s and the like.
+
 using Oceananigans
 using Oceananigans.Units
+
+# Next we `using ClimaSeaIce` to get some ice-specific names.
+
 using ClimaSeaIce
 
-# # Setting up a model of a bucket of water freezing from the top down
+# # An infinitely deep bucket with a single grid point
 #
-# We'll model a slab of ice in a bucket with one grid point, which
-# requires only a zero-dimensional grid.
+# Perhaps surprisingly, we need just one grid point
+# to model an possibly infinitely thick slab of ice with `SlabSeaIceModel`.
+# We would only need more than 1 grid point if our boundary conditions
+# vary in the horizontal direction.
 
 grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
 
-## Set the temperature at the top of the bucket to ``Tᵤ = -10ᵒC``.
+# Then our model of an ice slab freezing into a bucket. We set the top ice temperature
+# to `-10 ᵒC`. Note that other units besides Celsius _can_ be used, but that requires
+# setting model.phase_transitions` with appropriate parameters.
+
 model = SlabSeaIceModel(grid; top_thermal_boundary_condition=PrescribedTemperature(-10))
 
-## We'll freeze the bucket for 10 straight days
+# The default bottom thermal boundary condition for `SlabSeaIceModel` is
+# `IceWaterThermalEquilibrium` with freshwater. That's what we want!
+
+model.thermal_boundary_conditions.bottom
+
+# Ok, we're ready to freeze the bucket for 10 straight days with an initial ice
+# thickness of 1 cm,
+
 simulation = Simulation(model, Δt=10minute, stop_time=10days)
 
-## Initialize the ice thickness to 1 cm
 set!(model, h=0.01)
 
 # # Collecting data and running the simulation
 #
-# Before running the simulation, we set up a `Callback` to create
+# Before simulating the freezing bucket, we set up a `Callback` to create
 # a timeseries of the ice thickness saved at every time step.
 
 ## Container to hold the data
@@ -47,8 +76,8 @@ using CairoMakie
 
 # `timeseries` is a `Vector` of `Tuple`. So we have to do a bit of processing
 # to build `Vector`s of time `t` and thickness `h`. It's not much work though:
-t = map(ts -> ts[1], timeseries)
-h = map(ts -> ts[2], timeseries)
+t = [datum[1] for datum in timeseries]
+h = [datum[2] for datum in timeseries]
 
 # Just for fun, we also compute the velocity of the ice-water interface:
 dhdt = @. (h[2:end] - h[1:end-1]) / simulation.Δt
@@ -67,3 +96,6 @@ lines!(axd, 1e2 .* h[1:end-1], 1e6 .* dhdt)
 current_figure() # hide
 fig
 
+# If you want more ice, you can increase `simulation.stop_time` and
+# `run!(simulation)` again (or just re-run the whole script).
+

examples/melting_in_march.jl

@@ -4,11 +4,12 @@ using ClimaSeaIce
 using GLMakie
 
 # Generate a 0D grid for a single column slab model 
-grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
+grid = RectilinearGrid(size=4, x=(0, 1), topology=(Periodic, Flat, Flat))
 
 # Build a model of an ice slab that has internal conductive fluxes
 # and that emits radiation from its top surface.
-solar_insolation = -500 # W m⁻²
+solar_insolation = [-600, -800, -1000, -1200] # W m⁻²
+solar_insolation = reshape(solar_insolation, (4, 1, 1))
 outgoing_radiation = RadiativeEmission()
 
 # Define a FluxFunction representing a sensible heat flux
@@ -35,18 +36,24 @@ end
 aerodynamic_flux = FluxFunction(sensible_heat_flux; parameters)
 
 top_thermal_flux = (outgoing_radiation, solar_insolation, aerodynamic_flux)
-model = SlabSeaIceModel(grid; top_thermal_flux)
+model = SlabSeaIceModel(grid;
+                        ice_consolidation_thickness = 0.01, # m
+                        top_thermal_flux)
 set!(model, h=1)
 
-simulation = Simulation(model, Δt=10minute, stop_time=1day)
+simulation = Simulation(model, Δt=10minute, stop_time=30days)
 
 # Accumulate data
 timeseries = []
 
 function accumulate_timeseries(sim)
-    T = model.top_temperature
+    T = model.top_surface_temperature
     h = model.ice_thickness
-    push!(timeseries, (time(sim), first(h), first(T)))
+    push!(timeseries, (time(sim),
+                       h[1, 1, 1], T[1, 1, 1],
+                       h[2, 1, 1], T[2, 1, 1],
+                       h[3, 1, 1], T[3, 1, 1],
+                       h[4, 1, 1], T[4, 1, 1]))
 end
 
 simulation.callbacks[:save] = Callback(accumulate_timeseries)
@@ -56,8 +63,14 @@ run!(simulation)
 # Extract and visualize data
 
 t = [datum[1] for datum in timeseries]
-h = [datum[2] for datum in timeseries]
-T = [datum[3] for datum in timeseries]
+h1 = [datum[2] for datum in timeseries]
+T1 = [datum[3] for datum in timeseries]
+h2 = [datum[4] for datum in timeseries]
+T2 = [datum[5] for datum in timeseries]
+h3 = [datum[6] for datum in timeseries]
+T3 = [datum[7] for datum in timeseries]
+h4 = [datum[8] for datum in timeseries]
+T4 = [datum[9] for datum in timeseries]
 
 set_theme!(Theme(fontsize=24, linewidth=4))
 
@@ -66,8 +79,17 @@ fig = Figure(resolution=(1000, 800))
 axT = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Top temperature (ᵒC)")
 axh = Axis(fig[2, 1], xlabel="Time (days)", ylabel="Ice thickness (m)")
 
-lines!(axT, t / day, T)
-lines!(axh, t / day, h)
+lines!(axT, t / day, T1)
+lines!(axh, t / day, h1)
+
+lines!(axT, t / day, T2)
+lines!(axh, t / day, h2)
+
+lines!(axT, t / day, T3)
+lines!(axh, t / day, h3)
+
+lines!(axT, t / day, T4)
+lines!(axh, t / day, h4)
 
 display(fig)
 

examples/perpetual_night.jl

@@ -18,7 +18,7 @@ simulation = Simulation(model, Δt=1hour, stop_time=40days)
 timeseries = []
 
 function accumulate_timeseries(sim)
-    T = model.top_temperature
+    T = model.top_surface_temperature
     h = model.ice_thickness
     push!(timeseries, (time(sim), first(h), first(T)))
 end

src/ClimaSeaIce.jl

@@ -117,7 +117,7 @@
     cᵢ = thermo.ice_heat_capacity
     cℓ = thermo.liquid_heat_capacity
 
-    return ρᵢ * ℒ₀ + (ρℓ * cℓ - ρᵢ * cᵢ) * (T - T₀)
+    return ρℓ * ℒ₀ + (ρℓ * cℓ - ρᵢ * cᵢ) * (T - T₀)
 end
 
 struct ForwardEulerTimestepper end
@@ -132,7 +132,7 @@
     PrescribedTemperature
 
 include("EnthalpyMethodSeaIceModels.jl")
-include("SlabSeaIceModels.jl")
+include("SlabSeaIceModels/SlabSeaIceModels.jl")
 
 using .EnthalpyMethodSeaIceModels: EnthalpyMethodSeaIceModel
 using .SlabSeaIceModels: SlabSeaIceModel, ConductiveFlux