From 18aa3b66ba9d723f61a02c7edca86a3f816a6360 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 17 Dec 2024 17:04:53 -0700 Subject: [PATCH 1/6] Some cleanup --- src/SeaIceDynamics/SeaIceDynamics.jl | 4 +- src/SeaIceDynamics/nothing_dynamics.jl | 5 -- src/sea_ice_ab2_time_stepping.jl | 60 -------------------- src/sea_ice_advection.jl | 11 ++-- src/sea_ice_rk3_time_stepping.jl | 75 ------------------------- src/sea_ice_time_stepping.jl | 72 +++++++++++++++++++++++- src/tracer_tendency_kernel_functions.jl | 7 +++ 7 files changed, 85 insertions(+), 149 deletions(-) delete mode 100644 src/SeaIceDynamics/nothing_dynamics.jl delete mode 100644 src/sea_ice_ab2_time_stepping.jl delete mode 100644 src/sea_ice_rk3_time_stepping.jl diff --git a/src/SeaIceDynamics/SeaIceDynamics.jl b/src/SeaIceDynamics/SeaIceDynamics.jl index 8c513feb..fe715bbc 100644 --- a/src/SeaIceDynamics/SeaIceDynamics.jl +++ b/src/SeaIceDynamics/SeaIceDynamics.jl @@ -23,6 +23,6 @@ using Oceananigans.Grids: architecture ## - ice-atmosphere boundary stress (provided as an external flux) ## - ocean dynamic surface -include("nothing_dynamics.jl") +compute_momentum_tendencies!(model) = nothing -end \ No newline at end of file +end diff --git a/src/SeaIceDynamics/nothing_dynamics.jl b/src/SeaIceDynamics/nothing_dynamics.jl deleted file mode 100644 index f0e69efa..00000000 --- a/src/SeaIceDynamics/nothing_dynamics.jl +++ /dev/null @@ -1,5 +0,0 @@ -#### -#### For a `Nothing` dynamics, nothing happens! -#### - -step_momentum!(model, ::Nothing, Δt) = nothing diff --git a/src/sea_ice_ab2_time_stepping.jl b/src/sea_ice_ab2_time_stepping.jl deleted file mode 100644 index 96bb05e1..00000000 --- a/src/sea_ice_ab2_time_stepping.jl +++ /dev/null @@ -1,60 +0,0 @@ -using ClimaSeaIce.SeaIceDynamics: step_momentum! - -const AB2SeaIceModel = SeaIceModel{<:Any, <:Any, <:Any, <:QuasiAdamsBashforth2TimeStepper} - -function time_step!(model::AB2SeaIceModel, Δt; euler=false, callbacks = []) - - # Be paranoid and update state at iteration 0 - model.clock.iteration == 0 && update_state!(model) - - ab2_timestepper = model.timestepper - - # Change the default χ if necessary, which occurs if: - # * We detect that the time-step size has changed. - # * We detect that this is the "first" time-step, which means we - # need to take an euler step. Note that model.clock.last_Δt is - # initialized as Inf - # * The user has passed euler=true to time_step! - euler = euler || (Δt != model.clock.last_Δt) - - # If euler, then set χ = -0.5 - minus_point_five = convert(eltype(model.grid), -0.5) - χ = ifelse(euler, minus_point_five, ab2_timestepper.χ) - - # Set time-stepper χ (this is used in ab2_step!, but may also be used elsewhere) - χ₀ = ab2_timestepper.χ # Save initial value - ab2_timestepper.χ = χ - - # Ensure zeroing out all previous tendency fields to avoid errors in - # case G⁻ includes NaNs. See https://github.com/CliMA/Oceananigans.jl/issues/2259 - if euler - @debug "Taking a forward Euler step." - for field in ab2_timestepper.G⁻ - !isnothing(field) && fill!(field, 0) - end - end - - compute_tracer_tendencies!(model) - step_tracers!(model, Δt, 1) - - # TODO: This is an implicit (or split-explicit) step to advance momentum! - step_momentum!(model, model.ice_dynamics, Δt) - - # Only the tracers are advanced through an AB2 scheme - # (velocities are stepped in the dynamics step) - # so only tracers' tendencies are stored - store_tendencies!(model) - - tick!(model.clock, Δt) - update_state!(model) - - return nothing -end - -function timestepping_coefficients(ts::QuasiAdamsBashforth2TimeStepper, args...) - χ = ts.χ - FT = eltype(χ) - α = + convert(FT, 1.5) + χ - β = - convert(FT, 0.5) + χ - return α, β -end diff --git a/src/sea_ice_advection.jl b/src/sea_ice_advection.jl index 732bef00..20101eb8 100644 --- a/src/sea_ice_advection.jl +++ b/src/sea_ice_advection.jl @@ -16,7 +16,6 @@ using Oceananigans.Advection: FluxFormAdvection, # A = ∇ ⋅ (uh) _advective_thickness_flux_x(args...) = advective_thickness_flux_x(args...) - _advective_thickness_flux_y(args...) = advective_thickness_flux_y(args...) _advective_thickness_flux_x(i, j, k, ibg::ImmersedBoundaryGrid, args...) = @@ -44,7 +43,9 @@ end ∇Uℵh = 1 / Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_thickness_flux_x, advection, U.u, ℵ, h) + δyᵃᶜᵃ(i, j, k, grid, _advective_thickness_flux_y, advection, U.v, ℵ, h)) - @inbounds ℵ⁻¹ = ifelse(ℵ[i, j, k] != 0, 1 / ℵ[i, j, k], zero(grid)) + ℵᵢ = @inbounds ℵ[i, j, k] + non_zero_area = @inbounds ℵᵢ != 0 + ℵ⁻¹ = ifelse(non_zero_area, 1 / ℵᵢ, zero(grid)) return ℵ⁻¹ * ∇Uℵh end @@ -54,7 +55,9 @@ end ∇Uℵh = 1 / Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_thickness_flux_x, advection.x, U.u, ℵ, h) + δyᵃᶜᵃ(i, j, k, grid, _advective_thickness_flux_y, advection.y, U.v, ℵ, h)) - @inbounds ℵ⁻¹ = ifelse(ℵ[i, j, k] != 0, 1 / ℵ[i, j, k], zero(grid)) + ℵᵢ = @inbounds ℵ[i, j, k] + non_zero_area = @inbounds ℵᵢ != 0 + ℵ⁻¹ = ifelse(non_zero_area, 1 / ℵᵢ, zero(grid)) return ℵ⁻¹ * ∇Uℵh end @@ -67,4 +70,4 @@ end @inline horizontal_div_Uc(i, j, k, grid, advection::FluxFormAdvection, U, c) = 1 / Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_tracer_flux_x, advection.x, U.u, c) + δyᵃᶜᵃ(i, j, k, grid, _advective_tracer_flux_y, advection.y, U.v, c)) - \ No newline at end of file + diff --git a/src/sea_ice_rk3_time_stepping.jl b/src/sea_ice_rk3_time_stepping.jl deleted file mode 100644 index 0e10bec7..00000000 --- a/src/sea_ice_rk3_time_stepping.jl +++ /dev/null @@ -1,75 +0,0 @@ -const RK3SeaIceModel = SeaIceModel{<:Any, <:Any, <:Any, <:RungeKutta3TimeStepper} - -function time_step!(model::RK3SeaIceModel, Δt; callbacks = []) - - # Be paranoid and update state at iteration 0, in case run! is not used: - model.clock.iteration == 0 && update_state!(model) - - γ¹ = model.timestepper.γ¹ - γ² = model.timestepper.γ² - γ³ = model.timestepper.γ³ - - ζ² = model.timestepper.ζ² - ζ³ = model.timestepper.ζ³ - - first_stage_Δt = γ¹ * Δt - second_stage_Δt = (γ² + ζ²) * Δt - third_stage_Δt = (γ³ + ζ³) * Δt - - # - # First stage - # - - compute_tracer_tendencies!(model) - step_tracers!(model, Δt, 1) - - # TODO: This is an implicit (or split-explicit) step to advance momentum! - # do we need to pass Δt here or first_stage_Δt? - step_momentum!(model, model.ice_dynamics, Δt) - store_tendencies!(model) - - tick!(model.clock, first_stage_Δt) - update_state!(model) - - # - # Second stage - # - - compute_tracer_tendencies!(model) - step_tracers!(model, Δt, 2) - - # TODO: This is an implicit (or split-explicit) step to advance momentum! - # do we need to pass Δt here or second_stage_Δt? - step_momentum!(model, model.ice_dynamics, Δt) - store_tendencies!(model) - - tick!(model.clock, second_stage_Δt) - update_state!(model) - - # - # Third stage - # - - compute_tracer_tendencies!(model) - step_tracers!(model, Δt, 3) - - # TODO: This is an implicit (or split-explicit) step to advance momentum! - # do we need to pass Δt here or third_stage_Δt? - step_momentum!(model, model.ice_dynamics, Δt) - store_tendencies!(model) - - tick!(model.clock, third_stage_Δt) - update_state!(model) - - return nothing -end - -function timestepping_coefficients(ts::RungeKutta3TimeStepper, substep) - if substep == 1 - return ts.γ¹, zero(ts.γ¹) - elseif substep == 2 - return ts.γ², ts.ζ² - elseif substep == 3 - return ts.γ³, ts.ζ³ - end -end \ No newline at end of file diff --git a/src/sea_ice_time_stepping.jl b/src/sea_ice_time_stepping.jl index d5a9fffe..1b068548 100644 --- a/src/sea_ice_time_stepping.jl +++ b/src/sea_ice_time_stepping.jl @@ -10,7 +10,6 @@ function step_tracers!(model::SIM, Δt, substep) G⁻ = model.timestepper.G⁻ α, β = timestepping_coefficients(model.timestepper, substep) - launch!(arch, grid, :xyz, _step_tracers!, h, ℵ, tracers, Gⁿ, G⁻, Δt, α, β) return nothing @@ -25,7 +24,6 @@ end Ghⁿ = Gⁿ.h Gℵⁿ = Gⁿ.ℵ - Gh⁻ = G⁻.h Gℵ⁻ = G⁻.ℵ @@ -65,8 +63,74 @@ end @inbounds G⁻[n][i, j, 1] = Gⁿ[n][i, j, 1] end +const RK3SeaIceModel = SeaIceModel{<:Any, <:Any, <:Any, <:RungeKutta3TimeStepper} + +function time_step!(model::RK3SeaIceModel, Δt; callbacks = []) + + # Be paranoid and update state at iteration 0, in case run! is not used: + model.clock.iteration == 0 && update_state!(model) + + γ¹ = model.timestepper.γ¹ + γ² = model.timestepper.γ² + γ³ = model.timestepper.γ³ + + ζ² = model.timestepper.ζ² + ζ³ = model.timestepper.ζ³ + + first_stage_Δt = γ¹ * Δt + second_stage_Δt = (γ² + ζ²) * Δt + third_stage_Δt = (γ³ + ζ³) * Δt + + # + # First stage + # + + compute_tendencies!(model) + step_tracers!(model, Δt, 1) + step_momentum!(model, model.ice_dynamics, Δt) + store_tendencies!(model) + + tick!(model.clock, first_stage_Δt) + update_state!(model) + + # + # Second stage + # + + compute_tracer_tendencies!(model) + step_tracers!(model, Δt, 2) + step_momentum!(model, model.ice_dynamics, Δt) + store_tendencies!(model) + + tick!(model.clock, second_stage_Δt) + update_state!(model) + + # + # Third stage + # + + compute_tracer_tendencies!(model) + step_tracers!(model, Δt, 3) + step_momentum!(model, model.ice_dynamics, Δt) + store_tendencies!(model) + + tick!(model.clock, third_stage_Δt) + update_state!(model) + + return nothing +end + +function timestepping_coefficients(ts::RungeKutta3TimeStepper, substep) + if substep == 1 + return ts.γ¹, zero(ts.γ¹) + elseif substep == 2 + return ts.γ², ts.ζ² + elseif substep == 3 + return ts.γ³, ts.ζ³ + end +end + function update_state!(model::SIM) - foreach(prognostic_fields(model)) do field mask_immersed_field!(field) end @@ -75,3 +139,5 @@ function update_state!(model::SIM) return nothing end + + diff --git a/src/tracer_tendency_kernel_functions.jl b/src/tracer_tendency_kernel_functions.jl index dd830f5a..1c833ca9 100644 --- a/src/tracer_tendency_kernel_functions.jl +++ b/src/tracer_tendency_kernel_functions.jl @@ -1,6 +1,12 @@ using Oceananigans.Advection using ClimaSeaIce.SeaIceThermodynamics: thickness_thermodynamic_tendency +function compute_tendencies!(model::SIM) + compute_tracer_tendencies!(model) + compute_momentum_tendencies!(model) + return nothing +end + function compute_tracer_tendencies!(model::SIM) grid = model.grid arch = architecture(grid) @@ -79,3 +85,4 @@ function ice_thickness_tendency(i, j, k, grid, clock, return Gh_advection + Gh_thermodynamics + Fh end + From 1e90cda6d272712698d5f43201a8f9d1289bdea9 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 17 Dec 2024 17:23:28 -0700 Subject: [PATCH 2/6] Overhaul --- examples/ice_advection.jl | 140 ++++++++++++++++++++++++ src/ClimaSeaIce.jl | 2 - src/sea_ice_model.jl | 14 +-- src/sea_ice_time_stepping.jl | 37 ++----- src/tracer_tendency_kernel_functions.jl | 2 +- 5 files changed, 157 insertions(+), 38 deletions(-) create mode 100644 examples/ice_advection.jl diff --git a/examples/ice_advection.jl b/examples/ice_advection.jl new file mode 100644 index 00000000..080f20eb --- /dev/null +++ b/examples/ice_advection.jl @@ -0,0 +1,140 @@ +using Oceananigans +using Oceananigans.Units +using ClimaSeaIce +using Printf +using ClimaSeaIce.SeaIceDynamics + +arch = CPU() +L = 512kilometers +𝓋ₒ = 0.01 # m / s maximum ocean speed +𝓋ₐ = 30.0 # m / s maximum atmospheric speed modifier +Cᴰ = 1.2e-3 # Atmosphere - sea ice drag coefficient +ρₐ = 1.3 # kg/m³ + +# 2 km domain +grid = RectilinearGrid(arch; + size = (256, 256), + x = (0, L), + y = (0, L), + topology = (Bounded, Bounded, Flat)) + +# Constant ocean velocities corresponding to a cyclonic eddy +Uₒ = XFaceField(grid) +Vₒ = YFaceField(grid) +set!(Uₒ, (x, y) -> 𝓋ₒ * (2y - L) / L) +set!(Vₒ, (x, y) -> 𝓋ₒ * (L - 2x) / L) + +Uₐ = XFaceField(grid) +Vₐ = YFaceField(grid) + +# Atmosphere - sea ice stress +τᵤ = Field(ρₐ * Cᴰ * sqrt(Uₐ^2 + Vₐ^2) * Uₐ) +τᵥ = Field(ρₐ * Cᴰ * sqrt(Uₐ^2 + Vₐ^2) * Vₐ) + +# Atmospheric velocities corresponding to an anticyclonic eddy moving north-east +@inline center(t) = 256kilometers + 51.2kilometers * t / 86400 +@inline radius(x, y, t) = sqrt((x - center(t))^2 + (y - center(t))^2) +@inline speed(x, y, t) = 1 / 100 * exp(- radius(x, y, t) / 100kilometers) + +@inline ua_time(x, y, t) = - 𝓋ₐ * speed(x, y, t) * (cosd(72) * (x - center(t)) + sind(72) * (y - center(t))) / 1000 +@inline va_time(x, y, t) = - 𝓋ₐ * speed(x, y, t) * (cosd(72) * (y - center(t)) - sind(72) * (x - center(t))) / 1000 + +# Initialize the stress at time t = 0 +set!(Uₐ, (x, y) -> ua_time(x, y, 0)) +set!(Vₐ, (x, y) -> va_time(x, y, 0)) +compute!(τᵤ) +compute!(τᵥ) + +##### +##### Numerical details +##### + +momentum = nothing +advection = WENO(; order = 7) + +u_bcs = FieldBoundaryConditions(top = nothing, bottom = nothing, + north = ValueBoundaryCondition(0), + south = ValueBoundaryCondition(0)) + +v_bcs = FieldBoundaryConditions(top = nothing, bottom = nothing, + west = ValueBoundaryCondition(0), + east = ValueBoundaryCondition(0)) + +# Define the model! +model = SeaIceModel(grid; advection, + boundary_conditions = (u = u_bcs, v = v_bcs)) + #top_u_stress = τᵤ, + #top_v_stress = τᵥ, + #ocean_velocities = (u = Uₒ, v = Vₒ), + #coriolis = FPlane(f = 1e-4)) + +# Initial height field with perturbations around 0.3 m +h₀(x, y) = 0.3 + 0.005 * (sin(60 * x / 1000kilometers) + sin(30 * y / 1000kilometers)) + +# We start with a concentration of ℵ = 1 +set!(model, h = h₀) +set!(model, ℵ = 1) + +##### +##### Setup the simulation +##### + +# run the model for 2 days +simulation = Simulation(model, Δt = 2minutes, stop_time = 2days) + +# Evolve the wind stress field in time: +function compute_wind_stress(sim) + time = sim.model.clock.time + @inline ua(x, y) = ua_time(x, y, time) + @inline va(x, y) = va_time(x, y, time) + set!(Uₐ, ua) + set!(Vₐ, va) + compute!(τᵤ) + compute!(τᵥ) + return nothing +end + +add_callback!(simulation, compute_wind_stress) + +# Container to hold the data +htimeseries = [] +ℵtimeseries = [] +utimeseries = [] +vtimeseries = [] + +## Callback function to collect the data from the `sim`ulation +function accumulate_timeseries(sim) + h = sim.model.ice_thickness + ℵ = sim.model.ice_concentration + u = sim.model.velocities.u + v = sim.model.velocities.v + push!(htimeseries, deepcopy(interior(h))) + push!(ℵtimeseries, deepcopy(interior(ℵ))) + push!(utimeseries, deepcopy(interior(u))) + push!(vtimeseries, deepcopy(interior(v))) +end + +wall_time = Ref(time_ns()) + +function progress(sim) + h = sim.model.ice_thickness + ℵ = sim.model.ice_concentration + u = sim.model.velocities.u + v = sim.model.velocities.v + + hmax = maximum(interior(h)) + ℵmin = minimum(interior(ℵ)) + step_time = 1e-9 * (time_ns() - wall_time[]) + + @info @sprintf("Time: %s, Iteration %d, Δt %s, max(h): %.2f, min(ℵ): %.2f, wtime: %s \n", + prettytime(sim), iteration(sim), prettytime(sim.Δt), + hmax, ℵmin, prettytime(step_time)) + + wall_time[] = time_ns() +end + +simulation.callbacks[:progress] = Callback(progress, IterationInterval(5)) +simulation.callbacks[:save] = Callback(accumulate_timeseries, IterationInterval(5)) + +run!(simulation) + diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl index 25d3deea..b6c957aa 100644 --- a/src/ClimaSeaIce.jl +++ b/src/ClimaSeaIce.jl @@ -45,8 +45,6 @@ include("sea_ice_model.jl") include("sea_ice_advection.jl") include("tracer_tendency_kernel_functions.jl") include("sea_ice_time_stepping.jl") -include("sea_ice_ab2_time_stepping.jl") -include("sea_ice_rk3_time_stepping.jl") include("EnthalpyMethodSeaIceModel.jl") using .SeaIceThermodynamics diff --git a/src/sea_ice_model.jl b/src/sea_ice_model.jl index c9f1a2c3..c3a6f03e 100644 --- a/src/sea_ice_model.jl +++ b/src/sea_ice_model.jl @@ -15,7 +15,7 @@ struct SeaIceModel{GR, TD, D, TS, CL, U, T, IT, IC, STF, SMS, A} <: AbstractMode ice_concentration :: IC # Thermodynamics ice_thermodynamics :: TD - ice_dynamics :: D + rheology :: D # External boundary conditions external_heat_fluxes :: STF external_momentum_stresses :: SMS @@ -32,13 +32,11 @@ function SeaIceModel(grid; top_heat_flux = nothing, bottom_heat_flux = 0, velocities = nothing, - timestepper = :RungeKutta3, advection = nothing, - top_momentum_stress = nothing, # Fix when introducing dynamics tracers = (), boundary_conditions = NamedTuple(), ice_thermodynamics = SlabSeaIceThermodynamics(grid), - ice_dynamics = nothing) + rheology = nothing) if isnothing(velocities) velocities = (u = ZeroField(), v=ZeroField(), w=ZeroField()) @@ -62,7 +60,7 @@ function SeaIceModel(grid; # TODO: should we have ice thickness and concentration as part of the tracers or # just additional fields of the sea ice model? tracers = merge(tracers, (; S = ice_salinity)) - timestepper = TimeStepper(timestepper, grid, prognostic_tracers) + timestepper = TimeStepper(:RungeKutta3, grid, prognostic_tracers) top_heat_flux = external_top_heat_flux(ice_thermodynamics, top_heat_flux) @@ -77,9 +75,9 @@ function SeaIceModel(grid; ice_thickness, ice_concentration, ice_thermodynamics, - ice_dynamics, + rheology, external_heat_fluxes, - top_momentum_stress, + nothing, timestepper, advection) end @@ -88,7 +86,7 @@ const SIM = SeaIceModel function set!(model::SIM; h=nothing, ℵ=nothing) !isnothing(h) && set!(model.ice_thickness, h) - !isnothing(ℵ) && set!(model.ice_conentration, ℵ) + !isnothing(ℵ) && set!(model.ice_concentration, ℵ) return nothing end diff --git a/src/sea_ice_time_stepping.jl b/src/sea_ice_time_stepping.jl index 1b068548..6e9341a3 100644 --- a/src/sea_ice_time_stepping.jl +++ b/src/sea_ice_time_stepping.jl @@ -1,4 +1,4 @@ -function step_tracers!(model::SIM, Δt, substep) +function rk3_step!(model::SIM, Δt, γ, ζ) grid = model.grid arch = architecture(grid) @@ -9,8 +9,8 @@ function step_tracers!(model::SIM, Δt, substep) Gⁿ = model.timestepper.Gⁿ G⁻ = model.timestepper.G⁻ - α, β = timestepping_coefficients(model.timestepper, substep) - launch!(arch, grid, :xyz, _step_tracers!, h, ℵ, tracers, Gⁿ, G⁻, Δt, α, β) + launch!(arch, grid, :xyz, _step_tracers!, h, ℵ, tracers, Gⁿ, G⁻, Δt, γ, ζ) + #launch!(arch, grid, :xyz, _rk3_step_fields!, h, ℵ, tracers, Gⁿ, G⁻, Δt, γ, ζ) return nothing end @@ -43,16 +43,15 @@ end end function store_tendencies!(model::SIM) - grid = model.grid arch = architecture(grid) Nx, Ny, _ = size(grid) Gⁿ = model.timestepper.Gⁿ G⁻ = model.timestepper.G⁻ - Nt = length(Gⁿ) + NG = length(Gⁿ) - params = KernelParameters((Nx, Ny, Nt), (0, 0, 0)) + params = KernelParameters((Nx, Ny, NG), (0, 0, 0)) launch!(arch, model.grid, params, _store_tendencies!, G⁻, Gⁿ) return nothing @@ -86,10 +85,8 @@ function time_step!(model::RK3SeaIceModel, Δt; callbacks = []) # compute_tendencies!(model) - step_tracers!(model, Δt, 1) - step_momentum!(model, model.ice_dynamics, Δt) + rk3_step!(model, Δt, γ¹, zero(γ¹)) store_tendencies!(model) - tick!(model.clock, first_stage_Δt) update_state!(model) @@ -97,11 +94,9 @@ function time_step!(model::RK3SeaIceModel, Δt; callbacks = []) # Second stage # - compute_tracer_tendencies!(model) - step_tracers!(model, Δt, 2) - step_momentum!(model, model.ice_dynamics, Δt) + compute_tendencies!(model) + rk3_step!(model, Δt, γ², ζ²) store_tendencies!(model) - tick!(model.clock, second_stage_Δt) update_state!(model) @@ -109,27 +104,15 @@ function time_step!(model::RK3SeaIceModel, Δt; callbacks = []) # Third stage # - compute_tracer_tendencies!(model) - step_tracers!(model, Δt, 3) - step_momentum!(model, model.ice_dynamics, Δt) + compute_tendencies!(model) + rk3_step!(model, Δt, γ³, ζ³) store_tendencies!(model) - tick!(model.clock, third_stage_Δt) update_state!(model) return nothing end -function timestepping_coefficients(ts::RungeKutta3TimeStepper, substep) - if substep == 1 - return ts.γ¹, zero(ts.γ¹) - elseif substep == 2 - return ts.γ², ts.ζ² - elseif substep == 3 - return ts.γ³, ts.ζ³ - end -end - function update_state!(model::SIM) foreach(prognostic_fields(model)) do field mask_immersed_field!(field) diff --git a/src/tracer_tendency_kernel_functions.jl b/src/tracer_tendency_kernel_functions.jl index 1c833ca9..89832843 100644 --- a/src/tracer_tendency_kernel_functions.jl +++ b/src/tracer_tendency_kernel_functions.jl @@ -3,7 +3,7 @@ using ClimaSeaIce.SeaIceThermodynamics: thickness_thermodynamic_tendency function compute_tendencies!(model::SIM) compute_tracer_tendencies!(model) - compute_momentum_tendencies!(model) + # compute_momentum_tendencies!(model) return nothing end From 47eabbf381cb2eaebfc9c209cc2854d960fc93ab Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 17 Dec 2024 17:32:29 -0700 Subject: [PATCH 3/6] Overhaul again --- src/{sea_ice_advection.jl => Advection.jl} | 5 +++ src/ClimaSeaIce.jl | 8 ++-- .../SeaIceDynamics.jl => Rheologies.jl} | 3 +- src/SeaIceModels/SeaIceModels.jl | 7 ++++ .../compute_tendencies.jl} | 37 ++----------------- src/{ => SeaIceModels}/sea_ice_model.jl | 0 .../sea_ice_time_stepping.jl | 0 7 files changed, 20 insertions(+), 40 deletions(-) rename src/{sea_ice_advection.jl => Advection.jl} (98%) rename src/{SeaIceDynamics/SeaIceDynamics.jl => Rheologies.jl} (97%) create mode 100644 src/SeaIceModels/SeaIceModels.jl rename src/{tracer_tendency_kernel_functions.jl => SeaIceModels/compute_tendencies.jl} (58%) rename src/{ => SeaIceModels}/sea_ice_model.jl (100%) rename src/{ => SeaIceModels}/sea_ice_time_stepping.jl (100%) diff --git a/src/sea_ice_advection.jl b/src/Advection.jl similarity index 98% rename from src/sea_ice_advection.jl rename to src/Advection.jl index 20101eb8..36672698 100644 --- a/src/sea_ice_advection.jl +++ b/src/Advection.jl @@ -1,3 +1,7 @@ +module Advection + +export div_Uℵh, horizontal_div_Uc + using Oceananigans.Operators using Oceananigans.ImmersedBoundaries using Oceananigans.Advection: FluxFormAdvection, @@ -71,3 +75,4 @@ end 1 / Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_tracer_flux_x, advection.x, U.u, c) + δyᵃᶜᵃ(i, j, k, grid, _advective_tracer_flux_y, advection.y, U.v, c)) +end # Advection diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl index b6c957aa..e023ed67 100644 --- a/src/ClimaSeaIce.jl +++ b/src/ClimaSeaIce.jl @@ -40,11 +40,9 @@ mask_immersed_field!(::ConstantField) = nothing mask_immersed_field!(::ZeroField) = nothing include("SeaIceThermodynamics/SeaIceThermodynamics.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("Rheologies.jl") +include("Advection.jl") +include("SeaIceModels/SeaIceModels.jl") include("EnthalpyMethodSeaIceModel.jl") using .SeaIceThermodynamics diff --git a/src/SeaIceDynamics/SeaIceDynamics.jl b/src/Rheologies.jl similarity index 97% rename from src/SeaIceDynamics/SeaIceDynamics.jl rename to src/Rheologies.jl index fe715bbc..3d065db8 100644 --- a/src/SeaIceDynamics/SeaIceDynamics.jl +++ b/src/Rheologies.jl @@ -1,4 +1,4 @@ -module SeaIceDynamics +module Rheologies using ClimaSeaIce @@ -26,3 +26,4 @@ using Oceananigans.Grids: architecture compute_momentum_tendencies!(model) = nothing end + diff --git a/src/SeaIceModels/SeaIceModels.jl b/src/SeaIceModels/SeaIceModels.jl new file mode 100644 index 00000000..585341df --- /dev/null +++ b/src/SeaIceModels/SeaIceModels.jl @@ -0,0 +1,7 @@ +module SeaIceModels + +include("sea_ice_model.jl") +include("compute_tendencies.jl") +include("sea_ice_time_stepping.jl") + +end # SeaIceModels diff --git a/src/tracer_tendency_kernel_functions.jl b/src/SeaIceModels/compute_tendencies.jl similarity index 58% rename from src/tracer_tendency_kernel_functions.jl rename to src/SeaIceModels/compute_tendencies.jl index 89832843..9d367885 100644 --- a/src/tracer_tendency_kernel_functions.jl +++ b/src/SeaIceModels/compute_tendencies.jl @@ -1,4 +1,4 @@ -using Oceananigans.Advection +using ClimaSeaIce.Advection: div_Uℵh, horizontal_div_Uc using ClimaSeaIce.SeaIceThermodynamics: thickness_thermodynamic_tendency function compute_tendencies!(model::SIM) @@ -23,7 +23,6 @@ function compute_tracer_tendencies!(model::SIM) model.ice_thermodynamics, model.external_heat_fluxes.top, model.external_heat_fluxes.bottom, - nothing, #model.forcing.h, fields(model)) return nothing @@ -38,36 +37,9 @@ end thermodynamics, top_external_heat_flux, bottom_external_heat_flux, - h_forcing, model_fields) i, j, k = @index(Global, NTuple) - - @inbounds Gⁿ.h[i, j, k] = ice_thickness_tendency(i, j, k, grid, clock, - velocities, - advection, - ice_thickness, - concentration, - thermodynamics, - top_external_heat_flux, - bottom_external_heat_flux, - h_forcing, - model_fields) - - @inbounds Gⁿ.ℵ[i, j, k] = - horizontal_div_Uc(i, j, k, grid, advection, velocities, concentration) -end - -# Thickness change due to accretion and melting, restricted by minimum allowable value -function ice_thickness_tendency(i, j, k, grid, clock, - velocities, - advection, - ice_thickness, - ice_concentration, - thermodynamics, - top_external_heat_flux, - bottom_external_heat_flux, - h_forcing, - model_fields) Gh_advection = - div_Uℵh(i, j, k, grid, advection, velocities, ice_concentration, ice_thickness) @@ -78,11 +50,8 @@ function ice_thickness_tendency(i, j, k, grid, clock, top_external_heat_flux, bottom_external_heat_flux, clock, model_fields) - - # Compute forcing - Fh = zero(grid) #h_forcing(i, j, grid, clock, model_fields) - - return Gh_advection + Gh_thermodynamics + Fh + @inbounds Gⁿ.h[i, j, k] = Gh_advection + Gh_thermodynamics + @inbounds Gⁿ.ℵ[i, j, k] = - horizontal_div_Uc(i, j, k, grid, advection, velocities, concentration) end diff --git a/src/sea_ice_model.jl b/src/SeaIceModels/sea_ice_model.jl similarity index 100% rename from src/sea_ice_model.jl rename to src/SeaIceModels/sea_ice_model.jl diff --git a/src/sea_ice_time_stepping.jl b/src/SeaIceModels/sea_ice_time_stepping.jl similarity index 100% rename from src/sea_ice_time_stepping.jl rename to src/SeaIceModels/sea_ice_time_stepping.jl From e2175049e33702cacdec057a517e3a23d245d4f6 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 17 Dec 2024 17:32:52 -0700 Subject: [PATCH 4/6] Delete ice_advection --- examples/ice_advection.jl | 140 -------------------------------------- 1 file changed, 140 deletions(-) delete mode 100644 examples/ice_advection.jl diff --git a/examples/ice_advection.jl b/examples/ice_advection.jl deleted file mode 100644 index 080f20eb..00000000 --- a/examples/ice_advection.jl +++ /dev/null @@ -1,140 +0,0 @@ -using Oceananigans -using Oceananigans.Units -using ClimaSeaIce -using Printf -using ClimaSeaIce.SeaIceDynamics - -arch = CPU() -L = 512kilometers -𝓋ₒ = 0.01 # m / s maximum ocean speed -𝓋ₐ = 30.0 # m / s maximum atmospheric speed modifier -Cᴰ = 1.2e-3 # Atmosphere - sea ice drag coefficient -ρₐ = 1.3 # kg/m³ - -# 2 km domain -grid = RectilinearGrid(arch; - size = (256, 256), - x = (0, L), - y = (0, L), - topology = (Bounded, Bounded, Flat)) - -# Constant ocean velocities corresponding to a cyclonic eddy -Uₒ = XFaceField(grid) -Vₒ = YFaceField(grid) -set!(Uₒ, (x, y) -> 𝓋ₒ * (2y - L) / L) -set!(Vₒ, (x, y) -> 𝓋ₒ * (L - 2x) / L) - -Uₐ = XFaceField(grid) -Vₐ = YFaceField(grid) - -# Atmosphere - sea ice stress -τᵤ = Field(ρₐ * Cᴰ * sqrt(Uₐ^2 + Vₐ^2) * Uₐ) -τᵥ = Field(ρₐ * Cᴰ * sqrt(Uₐ^2 + Vₐ^2) * Vₐ) - -# Atmospheric velocities corresponding to an anticyclonic eddy moving north-east -@inline center(t) = 256kilometers + 51.2kilometers * t / 86400 -@inline radius(x, y, t) = sqrt((x - center(t))^2 + (y - center(t))^2) -@inline speed(x, y, t) = 1 / 100 * exp(- radius(x, y, t) / 100kilometers) - -@inline ua_time(x, y, t) = - 𝓋ₐ * speed(x, y, t) * (cosd(72) * (x - center(t)) + sind(72) * (y - center(t))) / 1000 -@inline va_time(x, y, t) = - 𝓋ₐ * speed(x, y, t) * (cosd(72) * (y - center(t)) - sind(72) * (x - center(t))) / 1000 - -# Initialize the stress at time t = 0 -set!(Uₐ, (x, y) -> ua_time(x, y, 0)) -set!(Vₐ, (x, y) -> va_time(x, y, 0)) -compute!(τᵤ) -compute!(τᵥ) - -##### -##### Numerical details -##### - -momentum = nothing -advection = WENO(; order = 7) - -u_bcs = FieldBoundaryConditions(top = nothing, bottom = nothing, - north = ValueBoundaryCondition(0), - south = ValueBoundaryCondition(0)) - -v_bcs = FieldBoundaryConditions(top = nothing, bottom = nothing, - west = ValueBoundaryCondition(0), - east = ValueBoundaryCondition(0)) - -# Define the model! -model = SeaIceModel(grid; advection, - boundary_conditions = (u = u_bcs, v = v_bcs)) - #top_u_stress = τᵤ, - #top_v_stress = τᵥ, - #ocean_velocities = (u = Uₒ, v = Vₒ), - #coriolis = FPlane(f = 1e-4)) - -# Initial height field with perturbations around 0.3 m -h₀(x, y) = 0.3 + 0.005 * (sin(60 * x / 1000kilometers) + sin(30 * y / 1000kilometers)) - -# We start with a concentration of ℵ = 1 -set!(model, h = h₀) -set!(model, ℵ = 1) - -##### -##### Setup the simulation -##### - -# run the model for 2 days -simulation = Simulation(model, Δt = 2minutes, stop_time = 2days) - -# Evolve the wind stress field in time: -function compute_wind_stress(sim) - time = sim.model.clock.time - @inline ua(x, y) = ua_time(x, y, time) - @inline va(x, y) = va_time(x, y, time) - set!(Uₐ, ua) - set!(Vₐ, va) - compute!(τᵤ) - compute!(τᵥ) - return nothing -end - -add_callback!(simulation, compute_wind_stress) - -# Container to hold the data -htimeseries = [] -ℵtimeseries = [] -utimeseries = [] -vtimeseries = [] - -## Callback function to collect the data from the `sim`ulation -function accumulate_timeseries(sim) - h = sim.model.ice_thickness - ℵ = sim.model.ice_concentration - u = sim.model.velocities.u - v = sim.model.velocities.v - push!(htimeseries, deepcopy(interior(h))) - push!(ℵtimeseries, deepcopy(interior(ℵ))) - push!(utimeseries, deepcopy(interior(u))) - push!(vtimeseries, deepcopy(interior(v))) -end - -wall_time = Ref(time_ns()) - -function progress(sim) - h = sim.model.ice_thickness - ℵ = sim.model.ice_concentration - u = sim.model.velocities.u - v = sim.model.velocities.v - - hmax = maximum(interior(h)) - ℵmin = minimum(interior(ℵ)) - step_time = 1e-9 * (time_ns() - wall_time[]) - - @info @sprintf("Time: %s, Iteration %d, Δt %s, max(h): %.2f, min(ℵ): %.2f, wtime: %s \n", - prettytime(sim), iteration(sim), prettytime(sim.Δt), - hmax, ℵmin, prettytime(step_time)) - - wall_time[] = time_ns() -end - -simulation.callbacks[:progress] = Callback(progress, IterationInterval(5)) -simulation.callbacks[:save] = Callback(accumulate_timeseries, IterationInterval(5)) - -run!(simulation) - From 59b8cc539250ead584e47f3ec45dd676e4a5a797 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 17 Dec 2024 17:36:28 -0700 Subject: [PATCH 5/6] Further cleanup --- src/ClimaSeaIce.jl | 4 +- src/sea_ice_model.jl | 126 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 129 insertions(+), 1 deletion(-) create mode 100644 src/sea_ice_model.jl diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl index e023ed67..e985cf92 100644 --- a/src/ClimaSeaIce.jl +++ b/src/ClimaSeaIce.jl @@ -1,4 +1,4 @@ -""" Ocean 🌊 Sea ice component of CliMa's Earth system model. """ +""" Ocean 🌊 Sea ice component of CliMA's Earth system model. """ module ClimaSeaIce using Oceananigans @@ -46,5 +46,7 @@ include("SeaIceModels/SeaIceModels.jl") include("EnthalpyMethodSeaIceModel.jl") using .SeaIceThermodynamics +using .SeaIceModels end # module + diff --git a/src/sea_ice_model.jl b/src/sea_ice_model.jl new file mode 100644 index 00000000..c3a6f03e --- /dev/null +++ b/src/sea_ice_model.jl @@ -0,0 +1,126 @@ +using Oceananigans.Fields: TracerFields +using Oceananigans.TimeSteppers: TimeStepper +using ClimaSeaIce.SeaIceThermodynamics: external_top_heat_flux +using Oceananigans: tupleit, tracernames +using ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions: flux_summary +using Oceananigans.Fields: ConstantField + +struct SeaIceModel{GR, TD, D, TS, CL, U, T, IT, IC, STF, SMS, A} <: AbstractModel{TS} + grid :: GR + clock :: CL + # Prognostic State + velocities :: U + tracers :: T + ice_thickness :: IT + ice_concentration :: IC + # Thermodynamics + ice_thermodynamics :: TD + rheology :: D + # External boundary conditions + external_heat_fluxes :: STF + external_momentum_stresses :: SMS + # Numerics + timestepper :: TS + advection :: A +end + +function SeaIceModel(grid; + clock = Clock{eltype(grid)}(time = 0), + ice_thickness = Field{Center, Center, Nothing}(grid), + ice_concentration = Field{Center, Center, Nothing}(grid), + ice_salinity = 0, # psu + top_heat_flux = nothing, + bottom_heat_flux = 0, + velocities = nothing, + advection = nothing, + tracers = (), + boundary_conditions = NamedTuple(), + ice_thermodynamics = SlabSeaIceThermodynamics(grid), + rheology = nothing) + + if isnothing(velocities) + velocities = (u = ZeroField(), v=ZeroField(), w=ZeroField()) + end + + tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example) + tracers = TracerFields(tracers, grid, boundary_conditions) + + # TODO: pass `clock` into `field`, so functions can be time-dependent? + # Wrap ice_salinity in a field + ice_salinity = field((Center, Center, Nothing), ice_salinity, grid) + + # Adding thickness and concentration if not there + prognostic_tracers = merge(tracers, (; h = ice_thickness, ℵ = ice_concentration)) + prognostic_tracers = if ice_salinity isa ConstantField + prognostic_tracers + else + merge(prognostic_tracers, (; S = ice_salinity)) + end + + # TODO: should we have ice thickness and concentration as part of the tracers or + # just additional fields of the sea ice model? + tracers = merge(tracers, (; S = ice_salinity)) + timestepper = TimeStepper(:RungeKutta3, grid, prognostic_tracers) + + top_heat_flux = external_top_heat_flux(ice_thermodynamics, top_heat_flux) + + # Package the external fluxes and boundary conditions + external_heat_fluxes = (top = top_heat_flux, + bottom = bottom_heat_flux) + + return SeaIceModel(grid, + clock, + velocities, + tracers, + ice_thickness, + ice_concentration, + ice_thermodynamics, + rheology, + external_heat_fluxes, + nothing, + timestepper, + advection) +end + +const SIM = SeaIceModel + +function set!(model::SIM; h=nothing, ℵ=nothing) + !isnothing(h) && set!(model.ice_thickness, h) + !isnothing(ℵ) && set!(model.ice_concentration, ℵ) + return nothing +end + +Base.summary(model::SIM) = "SeaIceModel" +prettytime(model::SIM) = prettytime(model.clock.time) +iteration(model::SIM) = model.clock.iteration + +function Base.show(io::IO, model::SIM) + grid = model.grid + arch = architecture(grid) + gridname = typeof(grid).name.wrapper + timestr = string("(time = ", prettytime(model), ", iteration = ", iteration(model), ")") + + print(io, "SeaIceModel{", typeof(arch), ", ", gridname, "}", timestr, '\n') + print(io, "├── grid: ", summary(model.grid), '\n') + print(io, "├── ice thermodynamics: ", summary(model.ice_thermodynamics), '\n') + print(io, "├── advection: ", summary(model.advection), '\n') + print(io, "└── external_heat_fluxes: ", '\n') + print(io, " ├── top: ", flux_summary(model.external_heat_fluxes.top, " │"), '\n') + print(io, " └── bottom: ", flux_summary(model.external_heat_fluxes.bottom, " ")) +end + +reset!(::SIM) = nothing +initialize!(::SIM) = nothing +default_included_properties(::SIM) = tuple(:grid) + +fields(model::SIM) = merge((; h = model.ice_thickness, + ℵ = model.ice_concentration), + model.tracers, + model.velocities, + fields(model.ice_thermodynamics)) + +# TODO: make this correct +prognostic_fields(model::SIM) = merge((; h = model.ice_thickness, + ℵ = model.ice_concentration), + model.tracers, + model.velocities) From b8a3d2739ac3ced543b6d1ee24bdae49dbdd4f7a Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 17 Dec 2024 17:38:07 -0700 Subject: [PATCH 6/6] Rm extra sea_ice_model --- src/sea_ice_model.jl | 126 ------------------------------------------- 1 file changed, 126 deletions(-) delete mode 100644 src/sea_ice_model.jl diff --git a/src/sea_ice_model.jl b/src/sea_ice_model.jl deleted file mode 100644 index c3a6f03e..00000000 --- a/src/sea_ice_model.jl +++ /dev/null @@ -1,126 +0,0 @@ -using Oceananigans.Fields: TracerFields -using Oceananigans.TimeSteppers: TimeStepper -using ClimaSeaIce.SeaIceThermodynamics: external_top_heat_flux -using Oceananigans: tupleit, tracernames -using ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions: flux_summary -using Oceananigans.Fields: ConstantField - -struct SeaIceModel{GR, TD, D, TS, CL, U, T, IT, IC, STF, SMS, A} <: AbstractModel{TS} - grid :: GR - clock :: CL - # Prognostic State - velocities :: U - tracers :: T - ice_thickness :: IT - ice_concentration :: IC - # Thermodynamics - ice_thermodynamics :: TD - rheology :: D - # External boundary conditions - external_heat_fluxes :: STF - external_momentum_stresses :: SMS - # Numerics - timestepper :: TS - advection :: A -end - -function SeaIceModel(grid; - clock = Clock{eltype(grid)}(time = 0), - ice_thickness = Field{Center, Center, Nothing}(grid), - ice_concentration = Field{Center, Center, Nothing}(grid), - ice_salinity = 0, # psu - top_heat_flux = nothing, - bottom_heat_flux = 0, - velocities = nothing, - advection = nothing, - tracers = (), - boundary_conditions = NamedTuple(), - ice_thermodynamics = SlabSeaIceThermodynamics(grid), - rheology = nothing) - - if isnothing(velocities) - velocities = (u = ZeroField(), v=ZeroField(), w=ZeroField()) - end - - tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example) - tracers = TracerFields(tracers, grid, boundary_conditions) - - # TODO: pass `clock` into `field`, so functions can be time-dependent? - # Wrap ice_salinity in a field - ice_salinity = field((Center, Center, Nothing), ice_salinity, grid) - - # Adding thickness and concentration if not there - prognostic_tracers = merge(tracers, (; h = ice_thickness, ℵ = ice_concentration)) - prognostic_tracers = if ice_salinity isa ConstantField - prognostic_tracers - else - merge(prognostic_tracers, (; S = ice_salinity)) - end - - # TODO: should we have ice thickness and concentration as part of the tracers or - # just additional fields of the sea ice model? - tracers = merge(tracers, (; S = ice_salinity)) - timestepper = TimeStepper(:RungeKutta3, grid, prognostic_tracers) - - top_heat_flux = external_top_heat_flux(ice_thermodynamics, top_heat_flux) - - # Package the external fluxes and boundary conditions - external_heat_fluxes = (top = top_heat_flux, - bottom = bottom_heat_flux) - - return SeaIceModel(grid, - clock, - velocities, - tracers, - ice_thickness, - ice_concentration, - ice_thermodynamics, - rheology, - external_heat_fluxes, - nothing, - timestepper, - advection) -end - -const SIM = SeaIceModel - -function set!(model::SIM; h=nothing, ℵ=nothing) - !isnothing(h) && set!(model.ice_thickness, h) - !isnothing(ℵ) && set!(model.ice_concentration, ℵ) - return nothing -end - -Base.summary(model::SIM) = "SeaIceModel" -prettytime(model::SIM) = prettytime(model.clock.time) -iteration(model::SIM) = model.clock.iteration - -function Base.show(io::IO, model::SIM) - grid = model.grid - arch = architecture(grid) - gridname = typeof(grid).name.wrapper - timestr = string("(time = ", prettytime(model), ", iteration = ", iteration(model), ")") - - print(io, "SeaIceModel{", typeof(arch), ", ", gridname, "}", timestr, '\n') - print(io, "├── grid: ", summary(model.grid), '\n') - print(io, "├── ice thermodynamics: ", summary(model.ice_thermodynamics), '\n') - print(io, "├── advection: ", summary(model.advection), '\n') - print(io, "└── external_heat_fluxes: ", '\n') - print(io, " ├── top: ", flux_summary(model.external_heat_fluxes.top, " │"), '\n') - print(io, " └── bottom: ", flux_summary(model.external_heat_fluxes.bottom, " ")) -end - -reset!(::SIM) = nothing -initialize!(::SIM) = nothing -default_included_properties(::SIM) = tuple(:grid) - -fields(model::SIM) = merge((; h = model.ice_thickness, - ℵ = model.ice_concentration), - model.tracers, - model.velocities, - fields(model.ice_thermodynamics)) - -# TODO: make this correct -prognostic_fields(model::SIM) = merge((; h = model.ice_thickness, - ℵ = model.ice_concentration), - model.tracers, - model.velocities)