User interface for specifying arbitrary boundary conditions

examples/diffusive_ice_column_model.jl

@@ -13,16 +13,48 @@ using GLMakie
 grid = RectilinearGrid(size=20, z=(-1, 0), topology=(Flat, Flat, Bounded))
 
 # Set up a simple problem and build the ice model
-atmosphere_temperature = -10  # ᵒC
-ocean_temperature      = 0.1  # ᵒC
 closure = MolecularDiffusivity(grid, κ_ice=1e-5, κ_water=1e-6)
-model = ThermodynamicSeaIceModel(; grid, closure, atmosphere_temperature, ocean_temperature)
+
+#####
+##### Create temperature boundary conditions
+#####
+
+initial_air_ice_temperature = -5
+top_T_amplitude = 5
+top_T_slope = -0.5 / day # ᵒC s⁻¹
+
+# Information about ocean cooling
+initial_ice_ocean_temperature = 1.1
+bottom_T_slope = -0.1 / day # ᵒC s⁻¹
+
+# Calculate BCs
+air_ice_temperature(x, y, t) = top_T_slope * t + top_T_amplitude * sin(2π*t/day) + initial_air_ice_temperature
+ice_ocean_temperature(x, y, t) = bottom_T_slope * t + initial_ice_ocean_temperature
+
+# Plot boundary condition functions
+dt = 10minutes
+tf = 10days
+t = 0:dt:tf
+
+set_theme!(Theme(fontsize=24, linewidth=3))
+fig = Figure()
+ax = Axis(fig[1, 1], title="Boundary Conditions", xlabel="Time (s)", ylabel="Temperature (ᵒC)")
+lines!(ax, t, air_ice_temperature.(0, 0, t), label="Air-ice surface temperature")
+lines!(ax, t, ice_ocean_temperature.(0, 0, t), label="Ice-ocean temperature")
+axislegend(ax)
+
+display(fig)
+
+top_T_bc = ValueBoundaryCondition(air_ice_temperature)
+bottom_T_bc = ValueBoundaryCondition(ice_ocean_temperature)
+T_bcs = FieldBoundaryConditions(top=top_T_bc, bottom=bottom_T_bc)
+
+model = ThermodynamicSeaIceModel(; grid, closure, boundary_conditions=(; T=T_bcs))
 
 # Initialize and run
-set!(model, T=ocean_temperature)
+set!(model, T=initial_ice_ocean_temperature)
 
 H = model.state.H
-@show H
 
 # We're using explicit time stepping. The CFL condition is
 #
@@ -33,8 +65,6 @@ H = model.state.H
 Δt = 0.1 * Δz^2 / κ
 simulation = Simulation(model; Δt)
 
-@show simulation
-
 #####
 ##### Set up diagnostics
 #####
@@ -96,13 +126,13 @@ function grab_profiles!(sim)
     return nothing
 end
 
-simulation.callbacks[:grabber] = Callback(grab_profiles!, SpecifiedTimes(10minutes, 30minutes, 1hour))
+simulation.callbacks[:grabber] = Callback(grab_profiles!, TimeInterval(1hour)) #SpecifiedTimes(10minutes, 30minutes, 1hour))
 
 #####
 ##### Run the simulation
 #####
 
-simulation.stop_time = 1hour
+simulation.stop_time = 10days
 run!(simulation)
 
 # Make a plot
@@ -114,6 +144,15 @@ axH = Axis(fig[1, 2], xlabel="Enthalpy (J m⁻³)", ylabel="z (m)")
 axϕ = Axis(fig[1, 3], xlabel="Porosity", ylabel="z (m)")
 axκ = Axis(fig[1, 4], xlabel="Diffusivity", ylabel="z (m)")
 axh = Axis(fig[2, 1:4], xlabel="Time (hours)", ylabel="Ice thickness (m)")
+axq = Axis(fig[3, 1:4], xlabel="Time (hours)", ylabel="Ice thickness (m)")
+
+xlims!(axT, -15, 1)
+xlims!(axH, -30, 10)
+xlims!(axϕ, -0.05, 1.05)
+
+Nt = length(tt)
+slider = Slider(fig[4, 1:4], range=1:Nt, startvalue=1)
+n = slider.value
 
 z = znodes(model.state.T)
 
@@ -125,21 +164,34 @@ c = ice_heat_capacity
 f(λ) = λ * exp(λ^2) * erf(λ) - St / sqrt(π)
 =#
 
-for n = 1:length(tt)
-    tn = tt[n]
-    Tn = Tt[n]
-    Hn = Ht[n]
-    ϕn = ϕt[n]
-    κn = κt[n]
-    label = "t = " * prettytime(tn)
-    scatterlines!(axT, Tn, z; label)
-    scatterlines!(axH, Hn, z; label)
-    scatterlines!(axϕ, ϕn, z; label)
-    scatterlines!(axκ, κn, z; label)
-end
+tn = @lift tt[$n]
+tnh = @lift tt[$n] / hour
+Tn = @lift Tt[$n]
+Hn = @lift Ht[$n]
+ϕn = @lift ϕt[$n]
+κn = @lift κt[$n]
+label = @lift "t = " * prettytime(tt[$n])
+scatterlines!(axT, Tn, z; label)
+scatterlines!(axH, Hn, z; label)
+scatterlines!(axϕ, ϕn, z; label)
+scatterlines!(axκ, κn, z; label)
 
+lines!(axh, th ./ hour, ht)
 lines!(axh, th ./ hour, ht)
 
+hn = @lift begin
+    m = searchsortedfirst(th, tt[$n])
+    ht[m]
+end
+
+scatter!(axh, tnh, hn, markersize=30, color=(:red, 0.5))
+
+lines!(axq, th ./ hour, air_ice_temperature.(0, 0, th), label="Air-ice surface temperature")
+lines!(axq, th ./ hour, ice_ocean_temperature.(0, 0, th), label="Ice-ocean temperature")
+vlines!(axq, tnh)
+axislegend(axq, position=:lb)
+
 axislegend(axT, position=:lb)
+
 display(fig)
 

src/ClimaSeaIce.jl

@@ -4,13 +4,13 @@ using Oceananigans.BoundaryConditions:
     fill_halo_regions!,
     regularize_field_boundary_conditions,
     FieldBoundaryConditions,
+    apply_z_bcs!,
     ValueBoundaryCondition
 
-using Oceananigans.Utils: prettysummary, prettytime
-using Oceananigans.Fields: CenterField, ZFaceField, Field, Center, Face, interior
+using Oceananigans.Utils: prettysummary, prettytime, launch!
+using Oceananigans.Fields: CenterField, ZFaceField, Field, Center, Face, interior, TracerFields
 using Oceananigans.Models: AbstractModel
 using Oceananigans.TimeSteppers: Clock, tick!
-using Oceananigans.Utils: launch!
 using Oceananigans.Operators
 
 using KernelAbstractions: @kernel, @index
@@ -28,7 +28,6 @@ mutable struct ThermodynamicSeaIceModel{Grid,
                                         State,
                                         Cp,
                                         Fu,
-                                        Ocean,
                                         Tend} <: AbstractModel{Nothing}
     grid :: Grid
     timestepper :: Tim # unused placeholder for now
@@ -38,7 +37,6 @@ mutable struct ThermodynamicSeaIceModel{Grid,
     ice_heat_capacity :: Cp
     water_heat_capacity :: Cp
     fusion_enthalpy :: Fu
-    ocean_temperature :: Ocean
     tendencies :: Tend
 end
 
@@ -47,8 +45,7 @@ const TSIM = ThermodynamicSeaIceModel
 function Base.show(io::IO, model::TSIM)
     clock = model.clock
     print(io, "ThermodynamicSeaIceModel(t=", prettytime(clock.time), ", iteration=", clock.iteration, ")", '\n')
-    print(io, "    grid: ", summary(model.grid), '\n')
-    print(io, "    ocean_temperature: ", model.ocean_temperature)
+    print(io, "    grid: ", summary(model.grid))
 end
 
 const reference_density = 999.8 # kg m⁻³
@@ -63,18 +60,16 @@ function ThermodynamicSeaIceModel(; grid,
                                   ice_heat_capacity = 2090.0 / reference_density,
                                   water_heat_capacity = 3991.0 / reference_density,
                                   fusion_enthalpy = 3.3e5 / reference_density,
-                                  atmosphere_temperature = -10, # ᵒC
-                                  ocean_temperature = 0) # ᵒC
-
-    # Build temperature field
-    top_T_bc = ValueBoundaryCondition(atmosphere_temperature)
-    bottom_T_bc = ValueBoundaryCondition(ocean_temperature)
-    T_location = (Center, Center, Center)
-    T_bcs = FieldBoundaryConditions(grid, T_location, top=top_T_bc, bottom=bottom_T_bc)
-    temperature = CenterField(grid, boundary_conditions=T_bcs)
-    enthalpy = CenterField(grid)
-    porosity = CenterField(grid)
-    state = (T=temperature, H=enthalpy, ϕ=porosity)
+                                  boundary_conditions = NamedTuple())
+
+    # Prognostic fields: temperature, enthalpy, porosity
+    field_names = (:T, :H, :ϕ)
+
+    # Build user-defined boundary conditions
+    user_boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, field_names)
+
+    # Initialize the `state` NamedTuple
+    state = TracerFields(field_names, grid, user_boundary_conditions)
 
     tendencies = (; H=CenterField(grid))
     clock = Clock{eltype(grid)}(0, 0, 1)
@@ -87,7 +82,6 @@ function ThermodynamicSeaIceModel(; grid,
                                     ice_heat_capacity,
                                     water_heat_capacity,
                                     fusion_enthalpy,
-                                    ocean_temperature,
                                     tendencies)
 end
 
@@ -152,7 +146,7 @@ function update_temperature!(model)
     interior(T) .= interior(H) ./ c
 
     arch = model.grid.architecture
-    fill_halo_regions!(T, arch, model.clock, fields(model))
+    fill_halo_regions!(T, model.clock, fields(model))
 
     return nothing
 end
@@ -168,7 +162,7 @@ function update_enthalpy!(model)
     interior(H) .= c .* interior(T) + ℒ * interior(ϕ)
 
     arch = model.grid.architecture
-    fill_halo_regions!(H, arch, model.clock, fields(model))
+    fill_halo_regions!(H, model.clock, fields(model))
 
     return nothing
 end
@@ -186,33 +180,37 @@ end
 function time_step!(model, Δt; callbacks=nothing)
     grid = model.grid
     arch = grid.architecture
-    U = model.state
+    Ψ = model.state
     G = model.tendencies
     closure = model.closure
 
-    launch!(arch, grid, :xyz, compute_tendencies!, G, grid, closure, U)
-    launch!(arch, grid, :xyz, step_fields!, U, G, Δt)
+    launch!(arch, grid, :xyz, compute_tendencies!, G, grid, closure, Ψ)
+
+    # Calculate fluxes
+    apply_z_bcs!(G.H, Ψ.H, arch, model.clock, model.state)
+
+    launch!(arch, grid, :xyz, step_fields!, Ψ, G, Δt)
     update_state!(model)
     tick!(model.clock, Δt)
 
     return nothing
 end
 
-@kernel function step_fields!(U, G, Δt)
+@kernel function step_fields!(Ψ, G, Δt)
     i, j, k = @index(Global, NTuple)
-    @inbounds U.H[i, j, k] += Δt * G.H[i, j, k]
+    @inbounds Ψ.H[i, j, k] += Δt * G.H[i, j, k]
 end
 
 #####
 ##### Physics
 #####
 
 """ Calculate the right-hand-side of the free surface displacement (η) equation. """
-@kernel function compute_tendencies!(G, grid, closure, U)
+@kernel function compute_tendencies!(G, grid, closure, Ψ)
     i, j, k = @index(Global, NTuple)
 
     # Temperature tendency
-    @inbounds G.H[i, j, k] = ∂z_κ_∂z_T(i, j, k, grid, closure, U.T)
+    @inbounds G.H[i, j, k] = ∂z_κ_∂z_T(i, j, k, grid, closure, Ψ.T)
 end
 
 struct MolecularDiffusivity{C}