From 56931b01f8f30b140d44e6e785f1755cec7e1cf4 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 29 Jun 2023 10:43:37 -0600 Subject: [PATCH 1/6] Implement model kwarg allowing arbitrary user-defined boundary conditions --- examples/diffusive_ice_column_model.jl | 11 +++++++++-- src/ClimaSeaIce.jl | 10 +++++++++- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/examples/diffusive_ice_column_model.jl b/examples/diffusive_ice_column_model.jl index 83bff8c1..0fc81434 100644 --- a/examples/diffusive_ice_column_model.jl +++ b/examples/diffusive_ice_column_model.jl @@ -13,10 +13,17 @@ using GLMakie grid = RectilinearGrid(size=20, z=(-1, 0), topology=(Flat, Flat, Bounded)) # Set up a simple problem and build the ice model +closure = MolecularDiffusivity(grid, κ_ice=1e-5, κ_water=1e-6) + +# Temperature boundary conditions 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) +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) + +model = ThermodynamicSeaIceModel(; grid, closure, boundary_conditions=(; T=T_bcs)) # Initialize and run set!(model, T=ocean_temperature) diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl index 55c179bb..748b83df 100644 --- a/src/ClimaSeaIce.jl +++ b/src/ClimaSeaIce.jl @@ -7,7 +7,7 @@ using Oceananigans.BoundaryConditions: ValueBoundaryCondition using Oceananigans.Utils: prettysummary, prettytime -using Oceananigans.Fields: CenterField, ZFaceField, Field, Center, Face, interior +using Oceananigans.Fields: CenterField, ZFaceField, Field, Center, Face, interior, TracerFields using Oceananigans.Models: AbstractModel using Oceananigans.TimeSteppers: Clock, tick! using Oceananigans.Utils: launch! @@ -63,9 +63,16 @@ function ThermodynamicSeaIceModel(; grid, ice_heat_capacity = 2090.0 / reference_density, water_heat_capacity = 3991.0 / reference_density, fusion_enthalpy = 3.3e5 / reference_density, + boundary_conditions = NamedTuple(), atmosphere_temperature = -10, # ᵒC ocean_temperature = 0) # ᵒC + # Prognostic fields: enthalpy + field_names = (:T, :H, :ϕ) + boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, field_names) + state = TracerFields(field_names, grid, user_boundary_conditions) + + #= # Build temperature field top_T_bc = ValueBoundaryCondition(atmosphere_temperature) bottom_T_bc = ValueBoundaryCondition(ocean_temperature) @@ -75,6 +82,7 @@ function ThermodynamicSeaIceModel(; grid, enthalpy = CenterField(grid) porosity = CenterField(grid) state = (T=temperature, H=enthalpy, ϕ=porosity) + =# tendencies = (; H=CenterField(grid)) clock = Clock{eltype(grid)}(0, 0, 1) From ed6512b75321ef2c711470ed39bf72fbb56e5f71 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 29 Jun 2023 10:45:40 -0600 Subject: [PATCH 2/6] Bugfixes and cleanup --- examples/diffusive_ice_column_model.jl | 3 +-- src/ClimaSeaIce.jl | 29 +++++++------------------- 2 files changed, 9 insertions(+), 23 deletions(-) diff --git a/examples/diffusive_ice_column_model.jl b/examples/diffusive_ice_column_model.jl index 0fc81434..a0509fb1 100644 --- a/examples/diffusive_ice_column_model.jl +++ b/examples/diffusive_ice_column_model.jl @@ -20,8 +20,7 @@ atmosphere_temperature = -10 # ᵒC ocean_temperature = 0.1 # ᵒC 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) +T_bcs = FieldBoundaryConditions(top=top_T_bc, bottom=bottom_T_bc) model = ThermodynamicSeaIceModel(; grid, closure, boundary_conditions=(; T=T_bcs)) diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl index 748b83df..9db7c1e7 100644 --- a/src/ClimaSeaIce.jl +++ b/src/ClimaSeaIce.jl @@ -38,7 +38,6 @@ mutable struct ThermodynamicSeaIceModel{Grid, ice_heat_capacity :: Cp water_heat_capacity :: Cp fusion_enthalpy :: Fu - ocean_temperature :: Ocean tendencies :: Tend end @@ -47,8 +46,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,26 +61,16 @@ function ThermodynamicSeaIceModel(; grid, ice_heat_capacity = 2090.0 / reference_density, water_heat_capacity = 3991.0 / reference_density, fusion_enthalpy = 3.3e5 / reference_density, - boundary_conditions = NamedTuple(), - atmosphere_temperature = -10, # ᵒC - ocean_temperature = 0) # ᵒC + boundary_conditions = NamedTuple()) - # Prognostic fields: enthalpy + # Prognostic fields: temperature, enthalpy, porosity field_names = (:T, :H, :ϕ) - boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, field_names) - state = TracerFields(field_names, grid, user_boundary_conditions) - #= - # 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) - =# + # 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) @@ -95,7 +83,6 @@ function ThermodynamicSeaIceModel(; grid, ice_heat_capacity, water_heat_capacity, fusion_enthalpy, - ocean_temperature, tendencies) end From e429e61f99218965240070e0e960634e6a95459b Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 29 Jun 2023 10:48:54 -0600 Subject: [PATCH 3/6] Remove spurious type parameter --- src/ClimaSeaIce.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl index 9db7c1e7..3d8c91bd 100644 --- a/src/ClimaSeaIce.jl +++ b/src/ClimaSeaIce.jl @@ -28,7 +28,6 @@ mutable struct ThermodynamicSeaIceModel{Grid, State, Cp, Fu, - Ocean, Tend} <: AbstractModel{Nothing} grid :: Grid timestepper :: Tim # unused placeholder for now From 626c32ee08c4fe0ecf2b223cb94a1f9271960505 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 29 Jun 2023 11:16:24 -0600 Subject: [PATCH 4/6] Add enthalply fluxes --- src/ClimaSeaIce.jl | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl index 3d8c91bd..e0f4ab5f 100644 --- a/src/ClimaSeaIce.jl +++ b/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.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 @@ -180,21 +180,25 @@ 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 ##### @@ -202,11 +206,11 @@ end ##### """ 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} From c14632369acb2730b6096fb29e1a29d980b3b5cb Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 29 Jun 2023 12:37:18 -0600 Subject: [PATCH 5/6] Incorporate jbs idealized air-ice and ice-ocean state --- examples/diffusive_ice_column_model.jl | 92 +++++++++++++++++++------- src/ClimaSeaIce.jl | 4 +- 2 files changed, 71 insertions(+), 25 deletions(-) diff --git a/examples/diffusive_ice_column_model.jl b/examples/diffusive_ice_column_model.jl index a0509fb1..2dcff785 100644 --- a/examples/diffusive_ice_column_model.jl +++ b/examples/diffusive_ice_column_model.jl @@ -15,20 +15,46 @@ grid = RectilinearGrid(size=20, z=(-1, 0), topology=(Flat, Flat, Bounded)) # Set up a simple problem and build the ice model closure = MolecularDiffusivity(grid, κ_ice=1e-5, κ_water=1e-6) -# Temperature boundary conditions -atmosphere_temperature = -10 # ᵒC -ocean_temperature = 0.1 # ᵒC -top_T_bc = ValueBoundaryCondition(atmosphere_temperature) -bottom_T_bc = ValueBoundaryCondition(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 # @@ -39,8 +65,6 @@ H = model.state.H Δt = 0.1 * Δz^2 / κ simulation = Simulation(model; Δt) -@show simulation - ##### ##### Set up diagnostics ##### @@ -102,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 @@ -120,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) @@ -131,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) diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl index e0f4ab5f..faa2fa52 100644 --- a/src/ClimaSeaIce.jl +++ b/src/ClimaSeaIce.jl @@ -146,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 @@ -162,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 From ab800805e9344f3732a39f57f331d99141821ecb Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Fri, 7 Jul 2023 16:46:05 -0600 Subject: [PATCH 6/6] Fix plot label --- examples/diffusive_ice_column_model.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/diffusive_ice_column_model.jl b/examples/diffusive_ice_column_model.jl index 2dcff785..0c568405 100644 --- a/examples/diffusive_ice_column_model.jl +++ b/examples/diffusive_ice_column_model.jl @@ -144,7 +144,7 @@ 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)") +axq = Axis(fig[3, 1:4], xlabel="Time (hours)", ylabel="Surface temperatures (ᵒC)") xlims!(axT, -15, 1) xlims!(axH, -30, 10) @@ -190,7 +190,6 @@ lines!(axq, th ./ hour, air_ice_temperature.(0, 0, th), label="Air-ice surface t 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)