diff --git a/docs/make.jl b/docs/make.jl index 85a1476d..91d58d9b 100644 --- a/docs/make.jl +++ b/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" ] diff --git a/examples/freezing_bucket.jl b/examples/freezing_bucket.jl index 692645c0..102b0ae3 100644 --- a/examples/freezing_bucket.jl +++ b/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. diff --git a/examples/ice_advected_by_anticyclone.jl b/examples/ice_advected_by_anticyclone.jl index 9298ac21..4c811fbe 100644 --- a/examples/ice_advected_by_anticyclone.jl +++ b/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)) diff --git a/examples/ice_advected_on_coastline.jl b/examples/ice_advected_on_coastline.jl index 8e238d26..f4812e75 100644 --- a/examples/ice_advected_on_coastline.jl +++ b/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) \ No newline at end of file +# ![](sea_ice_advected_on_coastline.mp4) diff --git a/src/Rheologies/Rheologies.jl b/src/Rheologies/Rheologies.jl index dc82048b..f0915a5a 100644 --- a/src/Rheologies/Rheologies.jl +++ b/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 \ No newline at end of file diff --git a/src/Rheologies/brittle_bingham_maxwell_rheology.jl b/src/Rheologies/brittle_bingham_maxwell_rheology.jl new file mode 100644 index 00000000..9cc80df2 --- /dev/null +++ b/src/Rheologies/brittle_bingham_maxwell_rheology.jl @@ -0,0 +1,258 @@ +using Oceananigans.Grids: halo_size +using Oceananigans.Units + +struct BrittleBinghamMaxwellRheology{FT} + ice_ridging_strength :: FT # compressive strength + ice_compaction_hardening :: FT # compaction hardening + ice_cohesion :: FT # cohesion + undamaged_elastic_modulus :: FT # minimum plastic parameter (transitions to viscous behaviour) + undamaged_viscous_relaxation_time :: FT # minimum number of substeps expressed as the dynamic coefficient + ridging_ice_thickness :: FT # maximum number of substeps expressed as the dynamic coefficient + poisson_ratio :: FT + friction_coefficient :: FT + maximum_compressive_strength :: FT + healing_constant :: FT + damage_parameter :: FT +end + +function BrittleBinghamMaxwellRheology(FT::DataType = Float64; + ice_ridging_strength = 1e4, + ice_compaction_hardening = 20, + ice_cohesion = 5.8e3, + undamaged_elastic_modulus = 5.96e8, + undamaged_viscous_relaxation_time = 1e7, + ridging_ice_thickness = 1, + poisson_ratio = 1 / 3, + friction_coefficient = 0.7, + maximum_compressive_strength = 2.9e7, + healing_constant = 28days * 40 / 1000, # Ks * 40 K / 1000 + damage_parameter = 5) + + return BrittleBinghamMaxwellRheology(convert(FT, ice_ridging_strength), + convert(FT, ice_compaction_hardening), + convert(FT, ice_cohesion), + convert(FT, undamaged_elastic_modulus), + convert(FT, undamaged_viscous_relaxation_time), + convert(FT, ridging_ice_thickness), + convert(FT, poisson_ratio), + convert(FT, friction_coefficient), + convert(FT, maximum_compressive_strength), + convert(FT, healing_constant), + convert(FT, damage_parameter)) +end + +rheology_prognostic_tracers(::BrittleBinghamMaxwellRheology) = (:d, :σ₁₁, :σ₁₂, :σ₂₂) + +function rheology_auxiliary_fields(::BrittleBinghamMaxwellRheology, grid) + + # TODO: What about boundary conditions? + P = Field{Center, Center, Nothing}(grid) + E = Field{Center, Center, Nothing}(grid) + λ = Field{Center, Center, Nothing}(grid) + σ¹ = Field{Center, Center, Nothing}(grid) + σ² = Field{Center, Center, Nothing}(grid) + return (; P, E, λ, σ¹, σ²) +end + +# Extend the `adapt_structure` function for the ElastoViscoPlasticRheology +Adapt.adapt_structure(to, r::BrittleBinghamMaxwellRheology) = + BrittleBinghamMaxwellRheology(Adapt.adapt(to, r.ice_ridging_strength), + Adapt.adapt(to, r.ice_compaction_hardening), + Adapt.adapt(to, r.ice_cohesion), + Adapt.adapt(to, r.undamaged_elastic_modulus), + Adapt.adapt(to, r.undamaged_viscous_relaxation_time), + Adapt.adapt(to, r.ridging_ice_thickness), + Adapt.adapt(to, r.poisson_ratio), + Adapt.adapt(to, r.friction_coefficient), + Adapt.adapt(to, r.maximum_compressive_strength), + Adapt.adapt(to, r.healing_constant), + Adapt.adapt(to, r.damage_parameter)) + +##### +##### Computation of the stresses +##### + +function initialize_rheology!(model, rheology::BrittleBinghamMaxwellRheology) + h = model.ice_thickness + ℵ = model.ice_concentration + + fields = model.dynamics.auxiliary_fields + grid = model.grid + + P★ = rheology.ice_ridging_strength + E★ = rheology.undamaged_elastic_modulus + λ★ = rheology.undamaged_viscous_relaxation_time + h★ = rheology.ridging_ice_thickness + Kh = rheology.healing_constant + α = rheology.damage_parameter + C = rheology.ice_compaction_hardening + d = model.tracers.d + Δt = model.clock.last_Δt + + ice_thermodynamics = model.ice_thermodynamics + + if !isnothing(ice_thermodynamics) + launch!(architecture(model.grid), model.grid, :xyz, _heal_outstanding_damage!, d, grid, ice_thermodynamics, Kh, Δt) + fill_halo_regions!(d) + end + + # compute on the whole grid including halos + params = KernelParameters(size(fields.P.data)[1:2], fields.P.data.offsets[1:2]) + launch!(architecture(model.grid), model.grid, params, _initialize_bbm_rheology!, fields, d, grid, P★, E★, λ★, h★, α, C, h, ℵ) + + return nothing +end + +@kernel function _heal_outstanding_damage!(d, grid, ice_thermodynamics, Kh, Δt) + i, j = @index(Global, NTuple) + + phase_transitions = ice_thermodynamics.phase_transitions + bottom_heat_bc = ice_thermodynamics.heat_boundary_conditions.bottom + liquidus = phase_transitions.liquidus + + Tu = @inbounds ice_thermodynamics.top_surface_temperature[i, j, 1] + Tb = bottom_temperature(i, j, grid, bottom_heat_bc, liquidus) + Tc = max((Tb - Tu) / Kh, 0) # restoring temperature + + @inbounds d[i, j, 1] -= Tc * Δt +end + +@kernel function _initialize_bbm_rheology!(fields, d, grid, P★, E★, λ★, h★, α, C, h, ℵ) + i, j = @index(Global, NTuple) + @inbounds begin + ecc = exp(- C * (1 - ℵ[i, j, 1])) + # Center - Center fields + fields.P[i, j, 1] = P★ * (h[i, j, 1] / h★)^(3/2) * ecc + fields.E[i, j, 1] = E★ * ecc + fields.λ[i, j, 1] = λ★ * ecc^(α - 1) + # Heal the damage with a healing constant + + # Clamp the damage between 0 and a value close to 1 + dᵢ = d[i, j, 1] + d[i, j, 1] = clamp(dᵢ, zero(dᵢ), 99999 * one(dᵢ) / 100000) + end +end + +function stress_params(grid) + TX, TY, _ = Oceananigans.Grids.topology(grid) + Fx = ifelse(TX == Flat, 1:1, 0:size(grid, 1)+1) + Fy = ifelse(TY == Flat, 1:1, 0:size(grid, 2)+1) + return KernelParameters(Fx, Fy) +end + +# Specific compute stresses for the EVP rheology +function compute_stresses!(model, dynamics, rheology::BrittleBinghamMaxwellRheology, Δt, Ns) + + grid = model.grid + arch = architecture(grid) + + ρᵢ = model.ice_density + u, v = model.velocities + fields = dynamics.auxiliary_fields + tracers = model.tracers + + params = stress_params(grid) + + # Pretty simple timestepping + Δτ = Δt / Ns + + launch!(arch, grid, params, _compute_stress_predictors!, fields, grid, rheology, tracers, u, v, ρᵢ, Δτ) + launch!(arch, grid, params, _advance_stresses_and_damage!, tracers, grid, rheology, fields, ρᵢ, Δτ) + + return nothing +end + +@inline function critical_damage(i, j, k, grid, N, c, μ, fields) + σ¹ = @inbounds fields.σ¹[i, j, k] + σ² = @inbounds fields.σ²[i, j, k] + dc = one(grid) - ifelse(σ¹ > - N, c / (σ² + μ * σ¹), - N / σ¹) + dc = ifelse(isnan(dc), zero(grid), dc) + return dc * (σ² > c - μ * σ¹) +end + +@kernel function _compute_stress_predictors!(fields, grid, rheology, tracers, u, v, ρᵢ, Δτ) + i, j = @index(Global, NTuple) + + α = rheology.damage_parameter + ν = rheology.poisson_ratio + + ϵ̇₁₁ = strain_rate_xx(i, j, 1, grid, u, v) + ϵ̇₂₂ = strain_rate_yy(i, j, 1, grid, u, v) + ϵ̇₁₂ = strain_rate_xy(i, j, 1, grid, u, v) + + Kϵ₁₁ = (ϵ̇₁₁ + ν * ϵ̇₂₂) / (1 - ν^2) + Kϵ₂₂ = (ϵ̇₂₂ + ν * ϵ̇₁₁) / (1 - ν^2) + Kϵ₁₂ = (1 - ν) * ϵ̇₁₂ / (1 - ν^2) + + dᵢ = @inbounds tracers.d[i, j, 1] + E₀ = @inbounds fields.E[i, j, 1] + λ₀ = @inbounds fields.λ[i, j, 1] + + E = E₀ * (1 - dᵢ) + λ = λ₀ * (1 - dᵢ)^(α - 1) + + # Test which isotropic stress to use + P̃ = zero(grid) + + # Implicit diagonal operator + Ω = 1 / (1 + Δτ * (1 + P̃) / λ) + + @inbounds tracers.σ₁₁[i, j, 1] = Ω * (tracers.σ₁₁[i, j, 1] + Δτ * E * Kϵ₁₁) + @inbounds tracers.σ₂₂[i, j, 1] = Ω * (tracers.σ₂₂[i, j, 1] + Δτ * E * Kϵ₂₂) + @inbounds tracers.σ₁₂[i, j, 1] = Ω * (tracers.σ₁₂[i, j, 1] + Δτ * E * Kϵ₁₂) + + σ₁₁ = @inbounds tracers.σ₁₁[i, j, 1] + σ₂₂ = @inbounds tracers.σ₂₂[i, j, 1] + σ₁₂ = @inbounds tracers.σ₁₂[i, j, 1] + + # Principal stress components + @inbounds fields.σ¹[i, j, 1] = (σ₁₁ + σ₂₂) / 2 + @inbounds fields.σ²[i, j, 1] = sqrt((σ₁₁ - σ₂₂)^2 / 4 + σ₁₂^2) +end + +@kernel function _advance_stresses_and_damage!(tracers, grid, rheology, fields, ρᵢ, Δτ) + i, j = @index(Global, NTuple) + + α = rheology.damage_parameter + ν = rheology.poisson_ratio + N = rheology.maximum_compressive_strength + c = rheology.ice_cohesion + μ = rheology.friction_coefficient + + ρ = @inbounds ρᵢ[i, j, 1] + dᵢ = @inbounds tracers.d[i, j, 1] + E₀ = @inbounds fields.E[i, j, 1] + + # The relaxation time is time-dependent + E = E₀ * (1 - dᵢ) + td = sqrt(2 * (1 + ν) * ρ / E * Azᶜᶜᶜ(i, j, 1, grid)) + + # damage tendency + Gd = reconstruction_2d(i, j, 1, grid, critical_damage, N, c, μ, fields) / td + + # Damage update + dᵢ += Gd * (1 - dᵢ) * Δτ + + # Advance stresses and clamp damage + @inbounds tracers.σ₁₁[i, j, 1] -= Gd * Δτ * tracers.σ₁₁[i, j, 1] + @inbounds tracers.σ₂₂[i, j, 1] -= Gd * Δτ * tracers.σ₂₂[i, j, 1] + @inbounds tracers.σ₁₂[i, j, 1] -= Gd * Δτ * tracers.σ₁₂[i, j, 1] + @inbounds tracers.d[i, j, 1] = clamp(dᵢ, zero(grid), 99999 * one(grid) / 100000) +end + +##### +##### Methods for the BBM rheology +##### + +# In the BBM rheology, the stresses need to be vertically integrated +@inline hσ₁₁(i, j, k, grid, fields) = @inbounds fields.σ₁₁[i, j, k] * fields.h[i, j, k] +@inline hσ₂₂(i, j, k, grid, fields) = @inbounds fields.σ₂₂[i, j, k] * fields.h[i, j, k] +@inline hσ₁₂(i, j, k, grid, fields) = @inbounds fields.σ₁₂[i, j, k] * fields.h[i, j, k] + +# Here we extend all the functions that a rheology model needs to support: +@inline ice_stress_ux(i, j, k, grid, ::BrittleBinghamMaxwellRheology, clock, fields) = ℑyᵃᶠᵃ(i, j, 1, grid, hσ₁₁, fields) +@inline ice_stress_uy(i, j, k, grid, ::BrittleBinghamMaxwellRheology, clock, fields) = ℑxᶠᵃᵃ(i, j, 1, grid, hσ₁₂, fields) +@inline ice_stress_vx(i, j, k, grid, ::BrittleBinghamMaxwellRheology, clock, fields) = ℑyᵃᶠᵃ(i, j, 1, grid, hσ₁₂, fields) +@inline ice_stress_vy(i, j, k, grid, ::BrittleBinghamMaxwellRheology, clock, fields) = ℑxᶠᵃᵃ(i, j, 1, grid, hσ₂₂, fields) + +@inline reconstruction_2d(i, j, k, grid, f, args...) = ℑxyᶜᶜᵃ(i, j, k, grid, ℑxyᶠᶠᵃ, f, args...) diff --git a/src/Rheologies/elasto_visco_plastic_rheology.jl b/src/Rheologies/elasto_visco_plastic_rheology.jl index ebcb984b..376151ba 100644 --- a/src/Rheologies/elasto_visco_plastic_rheology.jl +++ b/src/Rheologies/elasto_visco_plastic_rheology.jl @@ -105,19 +105,20 @@ function ElastoViscoPlasticRheology(FT::DataType = Float64; pressure_formulation) end -function required_auxiliary_fields(r::ElastoViscoPlasticRheology, grid) +function rheology_auxiliary_fields(r::ElastoViscoPlasticRheology, grid) # TODO: What about boundary conditions? σ₁₁ = Field{Center, Center, Nothing}(grid) σ₂₂ = Field{Center, Center, Nothing}(grid) - σ₁₂ = Field{Face, Face, Nothing}(grid) + σ₁₂ = Field{Center, Center, Nothing}(grid) - uⁿ = Field{Face, Center, Nothing}(grid) - vⁿ = Field{Center, Face, Nothing}(grid) P = Field{Center, Center, Nothing}(grid) - α = Field{Center, Center, Nothing}(grid) # Dynamic substeps a la Kimmritz et al (2016) ζ = Field{Center, Center, Nothing}(grid) Δ = Field{Center, Center, Nothing}(grid) + α = Field{Center, Center, Nothing}(grid) # Dynamic substeps a la Kimmritz et al (2016) + + uⁿ = Field{Face, Face, Nothing}(grid) + vⁿ = Field{Face, Face, Nothing}(grid) # An initial (safe) educated guess fill!(α, r.max_relaxation_parameter) @@ -170,7 +171,7 @@ end @inline ice_strength(i, j, k, grid, P★, C, h, ℵ) = @inbounds P★ * h[i, j, k] * exp(- C * (1 - ℵ[i, j, k])) # Specific compute stresses for the EVP rheology -function compute_stresses!(model, dynamics, rheology::ElastoViscoPlasticRheology, Δt) +function compute_stresses!(model, dynamics, rheology::ElastoViscoPlasticRheology, Δt, Ns) grid = model.grid arch = architecture(grid) @@ -187,17 +188,34 @@ function compute_stresses!(model, dynamics, rheology::ElastoViscoPlasticRheology parameters = KernelParameters(-Hx+2:Nx+Hx-1, -Hy+2:Ny+Hy-1) - launch!(arch, grid, parameters, _compute_evp_viscosities!, fields, grid, rheology, u, v) - launch!(arch, grid, parameters, _compute_evp_stresses!, fields, grid, rheology, u, v, h, ℵ, ρᵢ, Δt) + launch!(arch, grid, parameters, _compute_evp_viscosities!, fields, grid, rheology, u, v, h, ℵ, ρᵢ, Δt) return nothing end -@inline strain_rate_xx(i, j, k, grid, u, v) = δxᶜᵃᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u) / Azᶜᶜᶜ(i, j, k, grid) -@inline strain_rate_yy(i, j, k, grid, u, v) = δyᵃᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v) / Azᶜᶜᶜ(i, j, k, grid) -@inline strain_rate_xy(i, j, k, grid, u, v) = (δxᶠᵃᵃ(i, j, k, grid, Δy_qᶜᶠᶜ, v) + δyᵃᶠᵃ(i, j, k, grid, Δx_qᶠᶜᶜ, u)) / Azᶠᶠᶜ(i, j, k, grid) / 2 +@inline ice_pressure(i, j, k, grid, ::IceStrength, r, fields) = @inbounds fields.P[i, j, k] + +@inline function ice_pressure(i, j, k, grid, ::ReplacementPressure, r, fields) + Pᶜᶜᶜ = @inbounds fields.P[i, j, k] + Δᶜᶜᶜ = @inbounds fields.Δ[i, j, k] + Δm = r.minimum_plastic_stress + return Pᶜᶜᶜ * Δᶜᶜᶜ / (Δᶜᶜᶜ + Δm) +end + +@inline strain_rate_xx(i, j, k, grid, u, v) = ℑyᵃᶜᵃ(i, j, k, grid, δxᶜᵃᵃ, Δy_qᶠᶠᶜ, u) * Az⁻¹ᶜᶜᶜ(i, j, k, grid) +@inline strain_rate_yy(i, j, k, grid, u, v) = ℑxᶜᵃᵃ(i, j, k, grid, δyᵃᶜᵃ, Δx_qᶠᶠᶜ, v) * Az⁻¹ᶜᶜᶜ(i, j, k, grid) +@inline strain_rate_xy(i, j, k, grid, u, v) = (δxᶜᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, Δy_qᶠᶠᶜ, v) + δyᵃᶜᵃ(i, j, k, grid, ℑxᶜᵃᵃ, Δx_qᶠᶠᶜ, u)) * Az⁻¹ᶜᶜᶜ(i, j, k, grid) / 2 -@kernel function _compute_evp_viscosities!(fields, grid, rheology, u, v) +@inline ice_pressure(i, j, k, grid, ::IceStrength, r, fields) = @inbounds fields.P[i, j, k] + +@inline function ice_pressure(i, j, k, grid, ::ReplacementPressure, r, fields) + Pᶜᶜᶜ = @inbounds fields.P[i, j, k] + Δᶜᶜᶜ = @inbounds fields.Δ[i, j, k] + Δm = r.minimum_plastic_stress + return Pᶜᶜᶜ * Δᶜᶜᶜ / (Δᶜᶜᶜ + Δm) +end + +@kernel function _compute_evp_viscosities!(fields, grid, rheology, u, v, h, ℵ, ρᵢ, Δt) i, j = @index(Global, NTuple) kᴺ = size(grid, 3) @@ -210,41 +228,24 @@ end # Strain rates ϵ̇₁₁ = strain_rate_xx(i, j, kᴺ, grid, u, v) ϵ̇₂₂ = strain_rate_yy(i, j, kᴺ, grid, u, v) - - # Center - Center variables: - ϵ̇₁₂ᶜᶜᶜ = ℑxyᶜᶜᵃ(i, j, kᴺ, grid, strain_rate_xy, u, v) + ϵ̇₁₂ = strain_rate_xy(i, j, kᴺ, grid, u, v) # Ice divergence δ = ϵ̇₁₁ + ϵ̇₂₂ # Ice shear (at Centers) - s = sqrt((ϵ̇₁₁ - ϵ̇₂₂)^2 + 4ϵ̇₁₂ᶜᶜᶜ^2) + s = sqrt((ϵ̇₁₁ - ϵ̇₂₂)^2 + 4ϵ̇₁₂^2) # Visco - Plastic parameter # if Δ is very small we assume a linear viscous response # adding a minimum Δ_min (at Centers) - Δᶜᶜᶜ = max(sqrt(δ^2 + s^2 * e⁻²), Δm) - Pᶜᶜᶜ = @inbounds P[i, j, 1] + Δ = max(sqrt(δ^2 + s^2 * e⁻²), Δm) + P = @inbounds P[i, j, kᴺ] + ζ = P / 2Δ - @inbounds fields.ζ[i, j, 1] = Pᶜᶜᶜ / 2Δᶜᶜᶜ - @inbounds fields.Δ[i, j, 1] = Δᶜᶜᶜ -end - -@inline ice_pressure(i, j, k, grid, ::IceStrength, r, fields) = @inbounds fields.P[i, j, k] - -@inline function ice_pressure(i, j, k, grid, ::ReplacementPressure, r, fields) - Pᶜᶜᶜ = @inbounds fields.P[i, j, k] - Δᶜᶜᶜ = @inbounds fields.Δ[i, j, k] - Δm = r.minimum_plastic_stress - return Pᶜᶜᶜ * Δᶜᶜᶜ / (Δᶜᶜᶜ + Δm) -end - -# Compute the visco-plastic stresses for a slab sea ice model. -# The function updates the internal stress variables `σ₁₁`, `σ₂₂`, and `σ₁₂` in the `rheology` object -# following the αEVP formulation of Kimmritz et al (2016). -@kernel function _compute_evp_stresses!(fields, grid, rheology, u, v, h, ℵ, ρᵢ, Δt) - i, j = @index(Global, NTuple) - kᴺ = size(grid, 3) + # Store viscosity and the deformation for analysis purposes + @inbounds fields.ζ[i, j, kᴺ] = ζ + @inbounds fields.Δ[i, j, kᴺ] = Δ e⁻² = rheology.yield_curve_eccentricity^(-2) α⁺ = rheology.max_relaxation_parameter @@ -257,50 +258,36 @@ end σ₁₂ = fields.σ₁₂ α = fields.α - # Strain rates - ϵ̇₁₁ = strain_rate_xx(i, j, kᴺ, grid, u, v) - ϵ̇₂₂ = strain_rate_yy(i, j, kᴺ, grid, u, v) - ϵ̇₁₂ = strain_rate_xy(i, j, kᴺ, grid, u, v) - - ζᶜᶜᶜ = @inbounds fields.ζ[i, j, 1] - ζᶠᶠᶜ = ℑxyᶠᶠᵃ(i, j, 1, grid, fields.ζ) - # replacement pressure? Pᵣ = ice_pressure(i, j, 1, grid, ip, rheology, fields) - - ηᶜᶜᶜ = ζᶜᶜᶜ * e⁻² - ηᶠᶠᶜ = ζᶠᶠᶜ * e⁻² + η = ζ * e⁻² # σ(uᵖ): the tangential stress depends only shear viscosity # while the compressive stresses depend on the bulk viscosity and the ice strength - σ₁₁ᵖ⁺¹ = 2 * ηᶜᶜᶜ * ϵ̇₁₁ + ((ζᶜᶜᶜ - ηᶜᶜᶜ) * (ϵ̇₁₁ + ϵ̇₂₂) - Pᵣ / 2) - σ₂₂ᵖ⁺¹ = 2 * ηᶜᶜᶜ * ϵ̇₂₂ + ((ζᶜᶜᶜ - ηᶜᶜᶜ) * (ϵ̇₁₁ + ϵ̇₂₂) - Pᵣ / 2) - σ₁₂ᵖ⁺¹ = 2 * ηᶠᶠᶜ * ϵ̇₁₂ - - mᵢᶜᶜᶜ = ice_mass(i, j, 1, grid, h, ℵ, ρᵢ) - mᵢᶠᶠᶜ = ℑxyᶠᶠᵃ(i, j, 1, grid, ice_mass, h, ℵ, ρᵢ) + σ₁₁ᵖ⁺¹ = 2 * η * ϵ̇₁₁ + ((ζ - η) * (ϵ̇₁₁ + ϵ̇₂₂) - Pᵣ / 2) + σ₂₂ᵖ⁺¹ = 2 * η * ϵ̇₂₂ + ((ζ - η) * (ϵ̇₁₁ + ϵ̇₂₂) - Pᵣ / 2) + σ₁₂ᵖ⁺¹ = 2 * η * ϵ̇₁₂ + mᵢ = ice_mass(i, j, kᴺ, grid, h, ℵ, ρᵢ) + # Update coefficients for substepping using dynamic substepping # with spatially varying coefficients as in Kimmritz et al (2016) - γ²ᶜᶜᶜ = ζᶜᶜᶜ * cα * Δt / mᵢᶜᶜᶜ / Azᶜᶜᶜ(i, j, 1, grid) - γ²ᶜᶜᶜ = ifelse(isnan(γ²ᶜᶜᶜ), α⁺^2, γ²ᶜᶜᶜ) # In case both ζᶜᶜᶜ and mᵢᶜᶜᶜ are zero - γᶜᶜᶜ = clamp(sqrt(γ²ᶜᶜᶜ), α⁻, α⁺) - - γ²ᶠᶠᶜ = ζᶠᶠᶜ * cα * Δt / mᵢᶠᶠᶜ / Azᶠᶠᶜ(i, j, 1, grid) - γ²ᶠᶠᶜ = ifelse(isnan(γ²ᶠᶠᶜ), α⁺^2, γ²ᶠᶠᶜ) # In case both ζᶠᶠᶜ and mᵢᶠᶠᶜ are zero - γᶠᶠᶜ = clamp(sqrt(γ²ᶠᶠᶜ), α⁻, α⁺) + γ = ζ * cα * Δt / mᵢ * Az⁻¹ᶜᶜᶜ(i, j, kᴺ, grid) + α = clamp(sqrt(γ), α⁻, α⁺) + α = ifelse(isnan(α), α⁺, α) @inbounds begin # Compute the new stresses and store the value of the # dynamic substepping coefficient α - σ₁₁★ = (σ₁₁ᵖ⁺¹ - σ₁₁[i, j, 1]) / γᶜᶜᶜ - σ₂₂★ = (σ₂₂ᵖ⁺¹ - σ₂₂[i, j, 1]) / γᶜᶜᶜ - σ₁₂★ = (σ₁₂ᵖ⁺¹ - σ₁₂[i, j, 1]) / γᶠᶠᶜ - - σ₁₁[i, j, 1] += ifelse(mᵢᶜᶜᶜ > 0, σ₁₁★, zero(grid)) - σ₂₂[i, j, 1] += ifelse(mᵢᶜᶜᶜ > 0, σ₂₂★, zero(grid)) - σ₁₂[i, j, 1] += ifelse(mᵢᶠᶠᶜ > 0, σ₁₂★, zero(grid)) - α[i, j, 1] = γᶜᶜᶜ + σ₁₁★ = (σ₁₁ᵖ⁺¹ - σ₁₁[i, j, kᴺ]) / α + σ₂₂★ = (σ₂₂ᵖ⁺¹ - σ₂₂[i, j, kᴺ]) / α + σ₁₂★ = (σ₁₂ᵖ⁺¹ - σ₁₂[i, j, kᴺ]) / α + + σ₁₁[i, j, kᴺ] += ifelse(mᵢ > 0, σ₁₁★, zero(grid)) + σ₂₂[i, j, kᴺ] += ifelse(mᵢ > 0, σ₂₂★, zero(grid)) + σ₁₂[i, j, kᴺ] += ifelse(mᵢ > 0, σ₁₂★, zero(grid)) + + fields.α[i, j, kᴺ] = α end end @@ -309,27 +296,26 @@ end ##### # Here we extend all the functions that a rheology model needs to support: -@inline ice_stress_ux(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = @inbounds fields.σ₁₁[i, j, k] -@inline ice_stress_vx(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = @inbounds fields.σ₁₂[i, j, k] -@inline ice_stress_uy(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = @inbounds fields.σ₁₂[i, j, k] -@inline ice_stress_vy(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = @inbounds fields.σ₂₂[i, j, k] +@inline ice_stress_ux(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = ℑyᵃᶠᵃ(i, j, k, grid, fields.σ₁₁) +@inline ice_stress_uy(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = ℑxᶠᵃᵃ(i, j, k, grid, fields.σ₁₂) +@inline ice_stress_vx(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = ℑyᵃᶠᵃ(i, j, k, grid, fields.σ₁₂) +@inline ice_stress_vy(i, j, k, grid, ::ElastoViscoPlasticRheology, clock, fields) = ℑxᶠᵃᵃ(i, j, k, grid, fields.σ₂₂) # To help convergence to the right velocities -@inline compute_substep_Δtᶠᶜᶜ(i, j, grid, Δt, ::ElastoViscoPlasticRheology, substeps, fields) = Δt / ℑxᶠᵃᵃ(i, j, 1, grid, fields.α) -@inline compute_substep_Δtᶜᶠᶜ(i, j, grid, Δt, ::ElastoViscoPlasticRheology, substeps, fields) = Δt / ℑyᵃᶠᵃ(i, j, 1, grid, fields.α) +@inline compute_substep_Δtᶠᶠᶜ(i, j, grid, Δt, ::ElastoViscoPlasticRheology, substeps, fields) = @inbounds Δt / ℑxyᶠᶠᵃ(i, j, 1, grid, fields.α) ##### ##### Numerical forcing to help convergence ##### -@inline function sum_of_forcing_u(i, j, k, grid, ::ElastoViscoPlasticRheology, u_forcing, fields, Δt) +@inline function sum_of_forcing_u(i, j, k, grid, ::ElastoViscoPlasticRheology, u_forcing, fields, Δτ) user_forcing = u_forcing(i, j, k, grid, fields) - rheology_forcing = @inbounds (fields.uⁿ[i, j, k] - fields.u[i, j, k]) / Δt / ℑxᶠᵃᵃ(i, j, k, grid, fields.α) + rheology_forcing = @inbounds (fields.uⁿ[i, j, k] - fields.u[i, j, k]) / Δτ / ℑxyᶠᶠᵃ(i, j, k, grid, fields.α) return user_forcing + rheology_forcing end -@inline function sum_of_forcing_v(i, j, k, grid, ::ElastoViscoPlasticRheology, v_forcing, fields, Δt) +@inline function sum_of_forcing_v(i, j, k, grid, ::ElastoViscoPlasticRheology, v_forcing, fields, Δτ) user_forcing = v_forcing(i, j, k, grid, fields) - rheology_forcing = @inbounds (fields.vⁿ[i, j, k] - fields.v[i, j, k]) / Δt / ℑyᵃᶠᵃ(i, j, k, grid, fields.α) + rheology_forcing = @inbounds (fields.vⁿ[i, j, k] - fields.v[i, j, k]) / Δτ / ℑxyᶠᶠᵃ(i, j, k, grid, fields.α) return user_forcing + rheology_forcing end diff --git a/src/Rheologies/ice_stress_divergence.jl b/src/Rheologies/ice_stress_divergence.jl index 3515c1f0..9f923b11 100644 --- a/src/Rheologies/ice_stress_divergence.jl +++ b/src/Rheologies/ice_stress_divergence.jl @@ -2,6 +2,7 @@ using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, ImmersedBoundaryCondition using Oceananigans.Advection:conditional_flux_ccc, conditional_flux_ffc +using Oceananigans.Operators using Oceananigans.Operators: index_left, index_right const IBG = ImmersedBoundaryGrid @@ -19,77 +20,16 @@ const f = Face() @inline _ice_stress_vx(args...) = ice_stress_vx(args...) @inline _ice_stress_vy(args...) = ice_stress_vy(args...) -@inline _ice_stress_ux(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), ice_stress_ux(i, j, k, ibg, args...)) -@inline _ice_stress_uy(i, j, k, ibg::IBG, args...) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), ice_stress_uy(i, j, k, ibg, args...)) -@inline _ice_stress_vx(i, j, k, ibg::IBG, args...) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), ice_stress_vx(i, j, k, ibg, args...)) -@inline _ice_stress_vy(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), ice_stress_vy(i, j, k, ibg, args...)) - ##### ##### Stress divergence ##### @inline function ∂ⱼ_σ₁ⱼ(i, j, k, grid, rheology, clock, fields) - return 1 / Azᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Δy_qᶜᶜᶜ, _ice_stress_ux, rheology, clock, fields) + - δyᵃᶜᵃ(i, j, k, grid, Δx_qᶠᶠᶜ, _ice_stress_uy, rheology, clock, fields)) + return Az⁻¹ᶠᶠᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Δy_qᶜᶠᶜ, _ice_stress_ux, rheology, clock, fields) + + δyᵃᶠᵃ(i, j, k, grid, Δx_qᶠᶜᶜ, _ice_stress_uy, rheology, clock, fields)) end @inline function ∂ⱼ_σ₂ⱼ(i, j, k, grid, rheology, clock, fields) - return 1 / Azᶠᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Δy_qᶠᶠᶜ, _ice_stress_vx, rheology, clock, fields) + - δyᵃᶠᵃ(i, j, k, grid, Δx_qᶜᶜᶜ, _ice_stress_vy, rheology, clock, fields)) -end - -##### -##### Immersed Stress divergence -##### - -@inline immersed_∂ⱼ_σ₁ⱼ(i, j, k, grid, args...) = zero(grid) -@inline immersed_∂ⱼ_σ₂ⱼ(i, j, k, grid, args...) = zero(grid) - -@inline flip(::Type{Center}) = Face -@inline flip(::Type{Face}) = Center -@inline flip(::Center) = Face() -@inline flip(::Face) = Center() - -@inline function immersed_∂ⱼ_σ₁ⱼ(i, j, k, ibg::IBG, u_bc::IBC, rheology, clock, fields) - # Fetch fluxes across immersed boundary - q̃ᵂ = ib_ice_stress_ux(i, j, k, ibg, u_bc.west, rheology, clock, fields) - q̃ᴱ = ib_ice_stress_ux(i+1, j, k, ibg, u_bc.east, rheology, clock, fields) - q̃ˢ = ib_ice_stress_uy(i, j-1, k, ibg, u_bc.south, rheology, clock, fields) - q̃ᴺ = ib_ice_stress_uy(i, j, k, ibg, u_bc.north, rheology, clock, fields) - - iᵂ, jˢ, _ = map(index_left, (i, j, k), (c, c, c)) # Broadcast instead of map causes inference failure - iᴱ, jᴺ, _ = map(index_right, (i, j, k), (f, f, c)) - - # Impose i) immersed fluxes if we're on an immersed boundary or ii) zero otherwise. - qᵂ = conditional_flux_ccc(iᵂ, j, k, ibg, q̃ᵂ, zero(ibg)) * Axᶜᶜᶜ(iᵂ, j, k, grid) - qᴱ = conditional_flux_ccc(iᴱ, j, k, ibg, q̃ᴱ, zero(ibg)) * Axᶜᶜᶜ(iᴱ, j, k, grid) - qˢ = conditional_flux_ffc(i, jˢ, k, ibg, q̃ˢ, zero(ibg)) * Ayᶠᶠᶜ(i, jˢ, k, grid) - qᴺ = conditional_flux_ffc(i, jᴺ, k, ibg, q̃ᴺ, zero(ibg)) * Ayᶠᶠᶜ(i, jᴺ, k, grid) - - return (qᴱ - qᵂ + qᴺ - qˢ) / Vᶠᶜᶜ(i, j, k, grid) -end - -@inline function immersed_∂ⱼ_σ₂ⱼ(i, j, k, ibg::IBG, v_bc::IBC, rheology, clock, fields) - # Fetch fluxes across immersed boundary - q̃ᵂ = ib_ice_stress_vx(i-1, j, k, ibg, v_bc.west, rheology, clock, fields) - q̃ᴱ = ib_ice_stress_vx(i, j, k, ibg, v_bc.east, rheology, clock, fields) - q̃ˢ = ib_ice_stress_vy(i, j, k, ibg, v_bc.south, rheology, clock, fields) - q̃ᴺ = ib_ice_stress_vy(i, j+1, k, ibg, v_bc.north, rheology, clock, fields) - - iᵂ, jˢ, _ = map(index_left, (i, j, k), (f, f, c)) # Broadcast instead of map causes inference failure - iᴱ, jᴺ, _ = map(index_right, (i, j, k), (c, c, c)) - - # Impose i) immersed fluxes if we're on an immersed boundary or ii) zero otherwise. - qᵂ = conditional_flux_ffc(iᵂ, j, k, ibg, q̃ᵂ, zero(ibg)) * Axᶠᶠᶜ(iᵂ, j, k, grid) - qᴱ = conditional_flux_ffc(iᴱ, j, k, ibg, q̃ᴱ, zero(ibg)) * Axᶠᶠᶜ(iᴱ, j, k, grid) - qˢ = conditional_flux_ccc(i, jˢ, k, ibg, q̃ˢ, zero(ibg)) * Ayᶜᶜᶜ(i, jˢ, k, grid) - qᴺ = conditional_flux_ccc(i, jᴺ, k, ibg, q̃ᴺ, zero(ibg)) * Ayᶜᶜᶜ(i, jᴺ, k, grid) - - return (qᴱ - qᵂ + qᴺ - qˢ) / Vᶜᶠᶜ(i, j, k, grid) -end - -# TODO: Implement immersed fluxes (0 for the moment) -@inline ib_ice_stress_ux(i, j, k, grid, args...) = zero(grid) -@inline ib_ice_stress_vx(i, j, k, grid, args...) = zero(grid) -@inline ib_ice_stress_uy(i, j, k, grid, args...) = zero(grid) -@inline ib_ice_stress_vy(i, j, k, grid, args...) = zero(grid) \ No newline at end of file + return Az⁻¹ᶠᶠᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Δy_qᶜᶠᶜ, _ice_stress_vx, rheology, clock, fields) + + δyᵃᶠᵃ(i, j, k, grid, Δx_qᶠᶜᶜ, _ice_stress_vy, rheology, clock, fields)) +end \ No newline at end of file diff --git a/src/Rheologies/viscous_rheology.jl b/src/Rheologies/viscous_rheology.jl index c011f9eb..1d192e96 100644 --- a/src/Rheologies/viscous_rheology.jl +++ b/src/Rheologies/viscous_rheology.jl @@ -1,16 +1,31 @@ +using Base: Callable using Oceananigans.Fields: AbstractField, location -using Oceananigans.TurbulenceClosures: νᶜᶜᶜ, νᶠᶜᶠ, νᶜᶠᶠ, νᶠᶠᶜ, convert_diffusivity +using Oceananigans.TurbulenceClosures: convert_diffusivity struct ViscousRheology{N} ν :: N end +const c = Center() +const f = Face() + +@inline νᶠᶜᶜ(i, j, k, grid, loc, ν::Number, args...) = ν +@inline νᶜᶠᶜ(i, j, k, grid, loc, ν::Number, args...) = ν + +# Array / Field at `Center, Center, Center` +const Lᶜᶜᵃ = Tuple{Center, Center, <:Any} +@inline νᶠᶜᶜ(i, j, k, grid, ::Lᶜᶜᵃ, ν::AbstractArray, args...) = ℑxᵃᶠᵃ(i, j, k, grid, ν) +@inline νᶜᶠᶜ(i, j, k, grid, ::Lᶜᶜᵃ, ν::AbstractArray, args...) = ℑyᵃᶠᵃ(i, j, k, grid, ν) + +@inline νᶠᶜᶜ(i, j, k, grid, loc, ν::Callable, clock, args...) = ν(node(i, j, k, grid, f, c, c)..., clock.time) +@inline νᶜᶠᶜ(i, j, k, grid, loc, ν::Callable, clock, args...) = ν(node(i, j, k, grid, c, f, c)..., clock.time) + ViscousRheology(FT::DataType=Float64; ν = 1000.0) = ViscousRheology(convert_diffusivity(FT, ν)) @inline viscosity_location(ν) = (Center(), Center(), Center()) @inline viscosity_location(ν::AbstractField) = location(ν) -@inline ice_stress_ux(i, j, k, grid, vr::ViscousRheology, clock, fields) = νᶜᶜᶜ(i, j, k, grid, viscosity_location(vr.ν), vr.ν, clock, fields) * δxᶜᶜᶜ(i, j, k, grid, fields.u) -@inline ice_stress_vx(i, j, k, grid, vr::ViscousRheology, clock, fields) = νᶠᶠᶜ(i, j, k, grid, viscosity_location(vr.ν), vr.ν, clock, fields) * δxᶠᶠᶜ(i, j, k, grid, fields.v) -@inline ice_stress_uy(i, j, k, grid, vr::ViscousRheology, clock, fields) = νᶠᶠᶜ(i, j, k, grid, viscosity_location(vr.ν), vr.ν, clock, fields) * δyᶠᶠᶜ(i, j, k, grid, fields.u) -@inline ice_stress_vy(i, j, k, grid, vr::ViscousRheology, clock, fields) = νᶜᶜᶜ(i, j, k, grid, viscosity_location(vr.ν), vr.ν, clock, fields) * δyᶜᶜᶜ(i, j, k, grid, fields.v) +@inline ice_stress_ux(i, j, k, grid, vr::ViscousRheology, clock, fields) = νᶜᶠᶜ(i, j, k, grid, viscosity_location(vr.ν), vr.ν, clock, fields) * δxᶜᶠᶜ(i, j, k, grid, fields.u) +@inline ice_stress_vx(i, j, k, grid, vr::ViscousRheology, clock, fields) = νᶜᶠᶜ(i, j, k, grid, viscosity_location(vr.ν), vr.ν, clock, fields) * δxᶜᶠᶜ(i, j, k, grid, fields.v) +@inline ice_stress_uy(i, j, k, grid, vr::ViscousRheology, clock, fields) = νᶠᶜᶜ(i, j, k, grid, viscosity_location(vr.ν), vr.ν, clock, fields) * δyᶠᶜᶜ(i, j, k, grid, fields.u) +@inline ice_stress_vy(i, j, k, grid, vr::ViscousRheology, clock, fields) = νᶠᶜᶜ(i, j, k, grid, viscosity_location(vr.ν), vr.ν, clock, fields) * δyᶠᶜᶜ(i, j, k, grid, fields.v) diff --git a/src/SeaIceDynamics/SeaIceDynamics.jl b/src/SeaIceDynamics/SeaIceDynamics.jl index 4782e835..3fd4e736 100644 --- a/src/SeaIceDynamics/SeaIceDynamics.jl +++ b/src/SeaIceDynamics/SeaIceDynamics.jl @@ -16,17 +16,15 @@ using Oceananigans.Grids: architecture using ClimaSeaIce: ice_mass using ClimaSeaIce.Rheologies: ∂ⱼ_σ₁ⱼ, ∂ⱼ_σ₂ⱼ, - immersed_∂ⱼ_σ₁ⱼ, - immersed_∂ⱼ_σ₂ⱼ, - required_auxiliary_fields, + rheology_auxiliary_fields, compute_stresses!, initialize_rheology!, - compute_substep_Δtᶠᶜᶜ, - compute_substep_Δtᶜᶠᶜ, + compute_substep_Δtᶠᶠᶜ, sum_of_forcing_u, sum_of_forcing_v import Oceananigans: fields +import ClimaSeaIce.Rheologies: rheology_prognostic_tracers ## A Framework to solve for the ice momentum equation, in the form: ## diff --git a/src/SeaIceDynamics/explicit_momentum_equations.jl b/src/SeaIceDynamics/explicit_momentum_equations.jl index d83ab0fe..457b54ab 100644 --- a/src/SeaIceDynamics/explicit_momentum_equations.jl +++ b/src/SeaIceDynamics/explicit_momentum_equations.jl @@ -40,17 +40,15 @@ end i, j = @index(Global, NTuple) kᴺ = size(grid, 3) - ℵᶠᶜ = ℑxᶠᵃᵃ(i, j, kᴺ, grid, fields.ℵ) - ℵᶜᶠ = ℑyᵃᶠᵃ(i, j, kᴺ, grid, fields.ℵ) - mᶠᶜ = ℑxᶠᵃᵃ(i, j, kᴺ, grid, ice_mass, fields.h, fields.ℵ, fields.ρ) - mᶜᶠ = ℑyᵃᶠᵃ(i, j, kᴺ, grid, ice_mass, fields.h, fields.ℵ, fields.ρ) - + ℵᶠᶠ = ℑxyᶠᶠᵃ(i, j, 1, grid, fields.ℵ) + mᶠᶠ = ℑxyᶠᶠᵃ(i, j, 1, grid, ice_mass, fields.h, fields.ℵ, fields.ρ) + # Implicit part of the stress that depends linearly on the velocity τuᵢ = ( implicit_τx_coefficient(i, j, kᴺ, grid, bottom_stress, clock, fields) - - implicit_τx_coefficient(i, j, kᴺ, grid, top_stress, clock, fields)) / mᶠᶜ * ℵᶠᶜ + - implicit_τx_coefficient(i, j, kᴺ, grid, top_stress, clock, fields)) / mᶠᶠ * ℵᶠᶠ τvᵢ = ( implicit_τy_coefficient(i, j, kᴺ, grid, bottom_stress, clock, fields) - - implicit_τy_coefficient(i, j, kᴺ, grid, top_stress, clock, fields)) / mᶜᶠ * ℵᶜᶠ + - implicit_τy_coefficient(i, j, kᴺ, grid, top_stress, clock, fields)) / mᶠᶠ * ℵᶠᶠ @inbounds begin uᴰ = (u[i, j, 1] + Δt * Gⁿ.u[i, j, 1]) / (1 + Δt * τuᵢ) @@ -59,10 +57,9 @@ end uᶠ = free_drift_u(i, j, kᴺ, grid, free_drift, clock, fields) vᶠ = free_drift_v(i, j, kᴺ, grid, free_drift, clock, fields) - sea_ice = (mᶠᶜ ≥ minimum_mass) & (ℵᶠᶜ ≥ minimum_concentration) + sea_ice = (mᶠᶠ ≥ minimum_mass) & (ℵᶠᶠ ≥ minimum_concentration) + u[i, j, 1] = ifelse(sea_ice, uᴰ, uᶠ) - - sea_ice = (mᶜᶠ ≥ minimum_mass) & (ℵᶜᶠ ≥ minimum_concentration) v[i, j, 1] = ifelse(sea_ice, vᴰ, vᶠ) end end @@ -85,15 +82,11 @@ function compute_momentum_tendencies!(model, ::ExplicitMomentumEquation, Δt) top_stress = dynamics.external_momentum_stresses.top bottom_stress = dynamics.external_momentum_stresses.bottom - u_immersed_bc = model_fields.u.boundary_conditions.immersed - v_immersed_bc = model_fields.v.boundary_conditions.immersed - Gu = model.timestepper.Gⁿ.u Gv = model.timestepper.Gⁿ.v launch!(architecture(grid), grid, :xy, _compute_velocity_tendencies!, Gu, Gv, grid, Δt, rheology, model_fields, clock, coriolis, - u_immersed_bc, v_immersed_bc, top_stress, bottom_stress, model.forcing) return nothing @@ -101,14 +94,13 @@ end @kernel function _compute_velocity_tendencies!(Gu, Gv, grid, Δt, rheology, model_fields, clock, coriolis, - u_immersed_bc, v_immersed_bc, top_stress, bottom_stress, forcing) i, j = @index(Global, NTuple) @inbounds Gu[i, j, 1] = u_velocity_tendency(i, j, grid, Δt, rheology, model_fields, clock, coriolis, - u_immersed_bc, top_stress, bottom_stress, forcing.u) + top_stress, bottom_stress, forcing.u) @inbounds Gv[i, j, 1] = v_velocity_tendency(i, j, grid, Δt, rheology, model_fields, clock, coriolis, - v_immersed_bc, top_stress, bottom_stress, forcing.v) + top_stress, bottom_stress, forcing.v) end \ No newline at end of file diff --git a/src/SeaIceDynamics/momentum_tendencies_kernel_functions.jl b/src/SeaIceDynamics/momentum_tendencies_kernel_functions.jl index 84db232b..74a4e180 100644 --- a/src/SeaIceDynamics/momentum_tendencies_kernel_functions.jl +++ b/src/SeaIceDynamics/momentum_tendencies_kernel_functions.jl @@ -1,32 +1,29 @@ -using Oceananigans.Coriolis: y_f_cross_U, x_f_cross_U +using Oceananigans.Coriolis: fᶠᶠᵃ """compute explicit ice u-velocity tendencies""" -@inline function u_velocity_tendency(i, j, grid, Δt, +@inline function u_velocity_tendency(i, j, grid, Δτ, rheology, model_fields, clock, coriolis, - u_immersed_bc, u_top_stress, u_bottom_stress, u_forcing) kᴺ = size(grid, 3) - h = model_fields.h - ℵ = model_fields.ℵ - ρ = model_fields.ρ - U = (u = model_fields.u, v = model_fields.v) + h = model_fields.h + ℵ = model_fields.ℵ + ρ = model_fields.ρ # Ice mass (per unit area) interpolated on u points - ℵᵢ = ℑxᶠᵃᵃ(i, j, kᴺ, grid, ℵ) - mᵢ = ℑxᶠᵃᵃ(i, j, kᴺ, grid, ice_mass, h, ℵ, ρ) + ℵᵢ = ℑxyᶠᶠᵃ(i, j, 1, grid, ℵ) + mᵢ = ℑxyᶠᶠᵃ(i, j, 1, grid, ice_mass, h, ℵ, ρ) - Gᵁ = ( - x_f_cross_U(i, j, kᴺ, grid, coriolis, U) - - explicit_τx(i, j, kᴺ, grid, u_top_stress, clock, model_fields) / mᵢ * ℵᵢ - + explicit_τx(i, j, kᴺ, grid, u_bottom_stress, clock, model_fields) / mᵢ * ℵᵢ - + ∂ⱼ_σ₁ⱼ(i, j, kᴺ, grid, rheology, clock, model_fields) / mᵢ - + immersed_∂ⱼ_σ₁ⱼ(i, j, kᴺ, grid, u_immersed_bc, rheology, clock, model_fields) / mᵢ - + sum_of_forcing_u(i, j, kᴺ, grid, rheology, u_forcing, model_fields, Δt)) # sum of user defined forcing and possibly other forcing terms that are rheology-dependent + @inbounds Gᵁ = ( + fᶠᶠᵃ(i, j, 1, grid, coriolis) * model_fields.v[i, j, 1] + - explicit_τx(i, j, kᴺ, grid, u_top_stress, clock, model_fields) / mᵢ * ℵᵢ + + explicit_τx(i, j, kᴺ, grid, u_bottom_stress, clock, model_fields) / mᵢ * ℵᵢ + + ∂ⱼ_σ₁ⱼ(i, j, 1, grid, rheology, clock, model_fields) / mᵢ + + sum_of_forcing_u(i, j, kᴺ, grid, rheology, u_forcing, model_fields, Δτ)) # sum of user defined forcing and possibly other forcing terms that are rheology-dependent Gᵁ = ifelse(mᵢ ≤ 0, zero(grid), Gᵁ) @@ -34,12 +31,11 @@ using Oceananigans.Coriolis: y_f_cross_U, x_f_cross_U end """compute explicit ice v-velocity tendencies""" -@inline function v_velocity_tendency(i, j, grid, Δt, +@inline function v_velocity_tendency(i, j, grid, Δτ, rheology, model_fields, clock, coriolis, - v_immersed_bc, v_top_stress, v_bottom_stress, v_forcing) @@ -48,18 +44,16 @@ end h = model_fields.h ℵ = model_fields.ℵ ρ = model_fields.ρ - U = (u = model_fields.u, v = model_fields.v) # Ice mass (per unit area) interpolated on v points - ℵᵢ = ℑyᵃᶠᵃ(i, j, kᴺ, grid, ℵ) - mᵢ = ℑyᵃᶠᵃ(i, j, kᴺ, grid, ice_mass, h, ℵ, ρ) - - Gⱽ = ( - y_f_cross_U(i, j, kᴺ, grid, coriolis, U) - - explicit_τy(i, j, kᴺ, grid, v_top_stress, clock, model_fields) / mᵢ * ℵᵢ - + explicit_τy(i, j, kᴺ, grid, v_bottom_stress, clock, model_fields) / mᵢ * ℵᵢ - + ∂ⱼ_σ₂ⱼ(i, j, kᴺ, grid, rheology, clock, model_fields) / mᵢ - + immersed_∂ⱼ_σ₂ⱼ(i, j, kᴺ, grid, v_immersed_bc, rheology, clock, model_fields) / mᵢ - + sum_of_forcing_v(i, j, kᴺ, grid, rheology, v_forcing, model_fields, Δt)) # sum of user defined forcing and possibly other forcing terms that are rheology-dependent + ℵᵢ = ℑxyᶠᶠᵃ(i, j, 1, grid, ℵ) + mᵢ = ℑxyᶠᶠᵃ(i, j, 1, grid, ice_mass, h, ℵ, ρ) + + @inbounds Gⱽ = ( - fᶠᶠᵃ(i, j, 1, grid, coriolis) * model_fields.u[i, j, 1] + - explicit_τy(i, j, kᴺ, grid, v_top_stress, clock, model_fields) / mᵢ * ℵᵢ + + explicit_τy(i, j, kᴺ, grid, v_bottom_stress, clock, model_fields) / mᵢ * ℵᵢ + + ∂ⱼ_σ₂ⱼ(i, j, 1, grid, rheology, clock, model_fields) / mᵢ + + sum_of_forcing_v(i, j, kᴺ, grid, rheology, v_forcing, model_fields, Δτ)) # sum of user defined forcing and possibly other forcing terms that are rheology-dependent Gⱽ = ifelse(mᵢ ≤ 0, zero(grid), Gⱽ) diff --git a/src/SeaIceDynamics/sea_ice_external_stress.jl b/src/SeaIceDynamics/sea_ice_external_stress.jl index 37ca7802..5f9e43ef 100644 --- a/src/SeaIceDynamics/sea_ice_external_stress.jl +++ b/src/SeaIceDynamics/sea_ice_external_stress.jl @@ -90,26 +90,24 @@ Adapt.adapt_structure(to, τ::SemiImplicitStress) = @inline function explicit_τx(i, j, k, grid, τ::SemiImplicitStress, clock, fields) uₑ = @inbounds τ.uₑ[i, j, k] - Δu = @inbounds τ.uₑ[i, j, k] - fields.u[i, j, k] - Δv = ℑxyᶠᶜᵃ(i, j, k, grid, τ.vₑ) - ℑxyᶠᶜᵃ(i, j, k, grid, fields.v) - return τ.ρₑ * τ.Cᴰ * sqrt(Δu^2 + Δv^2) * uₑ + τi = implicit_τx_coefficient(i, j, k, grid, τ, clock, fields) + return τi * uₑ end @inline function explicit_τy(i, j, k, grid, τ::SemiImplicitStress, clock, fields) vₑ = @inbounds τ.vₑ[i, j, k] - Δv = @inbounds τ.vₑ[i, j, k] - fields.v[i, j, k] - Δu = ℑxyᶜᶠᵃ(i, j, k, grid, τ.uₑ) - ℑxyᶜᶠᵃ(i, j, k, grid, fields.u) - return τ.ρₑ * τ.Cᴰ * sqrt(Δu^2 + Δv^2) * vₑ + τi = implicit_τy_coefficient(i, j, k, grid, τ, clock, fields) + return τi * vₑ end @inline function implicit_τx_coefficient(i, j, k, grid, τ::SemiImplicitStress, clock, fields) Δu = @inbounds τ.uₑ[i, j, k] - fields.u[i, j, k] - Δv = ℑxyᶠᶜᵃ(i, j, k, grid, τ.vₑ) - ℑxyᶠᶜᵃ(i, j, k, grid, fields.v) + Δv = @inbounds τ.vₑ[i, j, k] - fields.v[i, j, k] return τ.ρₑ * τ.Cᴰ * sqrt(Δu^2 + Δv^2) end @inline function implicit_τy_coefficient(i, j, k, grid, τ::SemiImplicitStress, clock, fields) - Δu = ℑxyᶜᶠᵃ(i, j, k, grid, τ.uₑ) - ℑxyᶜᶠᵃ(i, j, k, grid, fields.u) + Δu = @inbounds τ.uₑ[i, j, k] - fields.u[i, j, k] Δv = @inbounds τ.vₑ[i, j, k] - fields.v[i, j, k] return τ.ρₑ * τ.Cᴰ * sqrt(Δu^2 + Δv^2) end diff --git a/src/SeaIceDynamics/sea_ice_momentum_equations.jl b/src/SeaIceDynamics/sea_ice_momentum_equations.jl index e30bebd9..910ab44f 100644 --- a/src/SeaIceDynamics/sea_ice_momentum_equations.jl +++ b/src/SeaIceDynamics/sea_ice_momentum_equations.jl @@ -66,7 +66,7 @@ function SeaIceMomentumEquation(grid; minimum_concentration = 1e-3, minimum_mass = 1.0) - auxiliary_fields = merge(auxiliary_fields, required_auxiliary_fields(rheology, grid)) + auxiliary_fields = merge(auxiliary_fields, rheology_auxiliary_fields(rheology, grid)) external_momentum_stresses = (top = top_momentum_stress, bottom = bottom_momentum_stress) @@ -83,3 +83,5 @@ function SeaIceMomentumEquation(grid; end fields(mom::SeaIceMomentumEquation) = mom.auxiliary_fields +rheology_prognostic_tracers(m::SeaIceMomentumEquation) = rheology_prognostic_tracers(m.rheology) + diff --git a/src/SeaIceDynamics/split_explicit_momentum_equations.jl b/src/SeaIceDynamics/split_explicit_momentum_equations.jl index 39884aeb..106c5b90 100644 --- a/src/SeaIceDynamics/split_explicit_momentum_equations.jl +++ b/src/SeaIceDynamics/split_explicit_momentum_equations.jl @@ -47,13 +47,12 @@ function time_step_momentum!(model, dynamics::SplitExplicitMomentumEquation, Δt u_forcing = model.forcing.u v_forcing = model.forcing.v - u_immersed_bc = u.boundary_conditions.immersed - v_immersed_bc = v.boundary_conditions.immersed model_fields = merge(dynamics.auxiliary_fields, model.velocities, (; h = model.ice_thickness, ℵ = model.ice_concentration, - ρ = model.ice_density)) + ρ = model.ice_density), + model.tracers) active_cells_map = Oceananigans.Grids.get_active_column_map(grid) @@ -67,7 +66,7 @@ function time_step_momentum!(model, dynamics::SplitExplicitMomentumEquation, Δt for substep in 1 : substeps # Compute stresses! depending on the particular rheology implementation - compute_stresses!(model, dynamics, rheology, Δt) + compute_stresses!(model, dynamics, rheology, Δt, substeps) # The momentum equations are solved using an alternating leap-frog algorithm # for u and v (used for the ocean - ice stresses and the coriolis term) @@ -77,24 +76,24 @@ function time_step_momentum!(model, dynamics::SplitExplicitMomentumEquation, Δt u_velocity_kernel!(u, grid, Δt, substeps, rheology, model_fields, free_drift, clock, coriolis, minimum_mass, minimum_concentration, - u_immersed_bc, top_stress, bottom_stress, u_forcing) + top_stress, bottom_stress, u_forcing) v_velocity_kernel!(v, grid, Δt, substeps, rheology, model_fields, free_drift, clock, coriolis, minimum_mass, minimum_concentration, - v_immersed_bc, top_stress, bottom_stress, v_forcing) + top_stress, bottom_stress, v_forcing) else v_velocity_kernel!(v, grid, Δt, substeps, rheology, model_fields, free_drift, clock, coriolis, minimum_mass, minimum_concentration, - v_immersed_bc, top_stress, bottom_stress, v_forcing) + top_stress, bottom_stress, v_forcing) u_velocity_kernel!(u, grid, Δt, substeps, rheology, model_fields, free_drift, clock, coriolis, minimum_mass, minimum_concentration, - u_immersed_bc, top_stress, bottom_stress, u_forcing) + top_stress, bottom_stress, u_forcing) end # TODO: This needs to be removed in some way! @@ -112,16 +111,16 @@ end model_fields, free_drift, clock, coriolis, minimum_mass, minimum_concentration, - u_immersed_bc, u_top_stress, u_bottom_stress, u_forcing) + u_top_stress, u_bottom_stress, u_forcing) i, j = @index(Global, NTuple) kᴺ = size(grid, 3) - mᵢ = ℑxᶠᵃᵃ(i, j, kᴺ, grid, ice_mass, model_fields.h, model_fields.ℵ, model_fields.ρ) - ℵᵢ = ℑxᶠᵃᵃ(i, j, kᴺ, grid, model_fields.ℵ) + mᵢ = ℑxyᶠᶠᵃ(i, j, 1, grid, ice_mass, model_fields.h, model_fields.ℵ, model_fields.ρ) + ℵᵢ = ℑxyᶠᶠᵃ(i, j, 1, grid, model_fields.ℵ) - Δτ = compute_substep_Δtᶠᶜᶜ(i, j, grid, Δt, rheology, substeps, model_fields) - Gu = u_velocity_tendency(i, j, grid, Δτ, rheology, model_fields, clock, coriolis, u_immersed_bc, u_top_stress, u_bottom_stress, u_forcing) + Δτ = compute_substep_Δtᶠᶠᶜ(i, j, grid, Δt, rheology, substeps, model_fields) + Gu = u_velocity_tendency(i, j, grid, Δτ, rheology, model_fields, clock, coriolis, u_top_stress, u_bottom_stress, u_forcing) # Implicit part of the stress that depends linearly on the velocity τuᵢ = ( implicit_τx_coefficient(i, j, kᴺ, grid, u_bottom_stress, clock, model_fields) @@ -143,16 +142,16 @@ end model_fields, free_drift, clock, coriolis, minimum_mass, minimum_concentration, - v_immersed_bc, v_top_stress, v_bottom_stress, v_forcing) + v_top_stress, v_bottom_stress, v_forcing) i, j = @index(Global, NTuple) kᴺ = size(grid, 3) - mᵢ = ℑyᵃᶠᵃ(i, j, kᴺ, grid, ice_mass, model_fields.h, model_fields.ℵ, model_fields.ρ) - ℵᵢ = ℑyᵃᶠᵃ(i, j, kᴺ, grid, model_fields.ℵ) + mᵢ = ℑxyᶠᶠᵃ(i, j, 1, grid, ice_mass, model_fields.h, model_fields.ℵ, model_fields.ρ) + ℵᵢ = ℑxyᶠᶠᵃ(i, j, 1, grid, model_fields.ℵ) - Δτ = compute_substep_Δtᶜᶠᶜ(i, j, grid, Δt, rheology, substeps, model_fields) - Gv = v_velocity_tendency(i, j, grid, Δτ, rheology, model_fields, clock, coriolis, v_immersed_bc, v_top_stress, v_bottom_stress, v_forcing) + Δτ = compute_substep_Δtᶠᶠᶜ(i, j, grid, Δt, rheology, substeps, model_fields) + Gv = v_velocity_tendency(i, j, grid, Δτ, rheology, model_fields, clock, coriolis, v_top_stress, v_bottom_stress, v_forcing) # Implicit part of the stress that depends linearly on the velocity τvᵢ = ( implicit_τy_coefficient(i, j, kᴺ, grid, v_bottom_stress, clock, model_fields) diff --git a/src/SeaIceDynamics/stress_balance_free_drift.jl b/src/SeaIceDynamics/stress_balance_free_drift.jl index 71995bf9..6522d5e3 100644 --- a/src/SeaIceDynamics/stress_balance_free_drift.jl +++ b/src/SeaIceDynamics/stress_balance_free_drift.jl @@ -51,7 +51,7 @@ const BISB = StressBalanceFreeDrift{<:SemiImplicitStress, <:Any} # Then: 𝒰ᵢ = 𝒰ᴮ - τᵀ / sqrt(Cᴮ * ||τᵀ||) @inline function free_drift_u(i, j, k, grid, f::TISB, clock, fields) τxᵀ = x_momentum_stress(i, j, k, grid, f.top_momentum_stress, clock, fields) - τyᵀ = ℑxyᶠᶜᵃ(i, j, k, grid, y_momentum_stress, f.top_momentum_stress, clock, fields) + τyᵀ = y_momentum_stress(i, j, k, grid, f.top_momentum_stress, clock, fields) τᵀ = sqrt(τxᵀ^2 + τyᵀ^2) τᴮ = f.bottom_momentum_stress @@ -62,7 +62,7 @@ const BISB = StressBalanceFreeDrift{<:SemiImplicitStress, <:Any} end @inline function free_drift_v(i, j, k, grid, f::TISB, clock, fields) - τxᵀ = ℑxyᶜᶠᵃ(i, j, k, grid, x_momentum_stress, f.top_momentum_stress, clock, fields) + τxᵀ = x_momentum_stress(i, j, k, grid, f.top_momentum_stress, clock, fields) τyᵀ = y_momentum_stress(i, j, k, grid, f.top_momentum_stress, clock, fields) τᵀ = sqrt(τxᵀ^2 + τyᵀ^2) @@ -77,7 +77,7 @@ end # Then: 𝒰ᵢ = 𝒰ᵀ - τᴮ / sqrt(Cᵀ * ||τᴮ||) @inline function free_drift_u(i, j, k, grid, f::BISB, clock, fields) τxᴮ = x_momentum_stress(i, j, k, grid, f.bottom_momentum_stress, clock, fields) - τyᴮ = ℑxyᶠᶜᵃ(i, j, k, grid, y_momentum_stress, f.bottom_momentum_stress, clock, fields) + τyᴮ = y_momentum_stress(i, j, k, grid, f.bottom_momentum_stress, clock, fields) τᴮ = sqrt(τxᴮ^2 + τyᴮ^2) τᵀ = f.top_momentum_stess @@ -88,7 +88,7 @@ end end @inline function free_drift_v(i, j, k, grid, f::BISB, clock, fields) - τxᴮ = ℑxyᶜᶠᵃ(i, j, k, grid, x_momentum_stress, f.bottom_momentum_stress, clock, fields) + τxᴮ = x_momentum_stress(i, j, k, grid, f.bottom_momentum_stress, clock, fields) τyᴮ = y_momentum_stress(i, j, k, grid, f.bottom_momentum_stress, clock, fields) τᴮ = sqrt(τxᴮ^2 + τyᴮ^2) diff --git a/src/SeaIceThermodynamics/slab_thermodynamics_tendencies.jl b/src/SeaIceThermodynamics/slab_thermodynamics_tendencies.jl index bff2482d..a4892b61 100644 --- a/src/SeaIceThermodynamics/slab_thermodynamics_tendencies.jl +++ b/src/SeaIceThermodynamics/slab_thermodynamics_tendencies.jl @@ -7,6 +7,7 @@ using Oceananigans ice_thickness, ice_concentration, ice_consolidation_thickness, + ice_salinity, top_external_heat_flux, bottom_external_heat_flux, clock, model_fields) diff --git a/src/SeaIceThermodynamics/thermodynamic_time_step.jl b/src/SeaIceThermodynamics/thermodynamic_time_step.jl index a453545f..53c23dd0 100644 --- a/src/SeaIceThermodynamics/thermodynamic_time_step.jl +++ b/src/SeaIceThermodynamics/thermodynamic_time_step.jl @@ -15,6 +15,7 @@ function thermodynamic_time_step!(model, ::SlabSeaIceThermodynamics, Δt) grid, Δt, model.clock, model.ice_consolidation_thickness, + model.tracers.S, model.ice_thermodynamics, model.external_heat_fluxes.top, model.external_heat_fluxes.bottom, @@ -40,6 +41,7 @@ end Δt, clock, ice_consolidation_thickness, + ice_salinity, ice_thermodynamics, top_external_heat_flux, bottom_external_heat_flux, @@ -58,6 +60,7 @@ end ice_thickness, ice_concentration, ice_consolidation_thickness, + ice_salinity, top_external_heat_flux, bottom_external_heat_flux, clock, model_fields) diff --git a/src/forward_euler_timestepper.jl b/src/forward_euler_timestepper.jl index 8fe7b17f..4b7e7242 100644 --- a/src/forward_euler_timestepper.jl +++ b/src/forward_euler_timestepper.jl @@ -11,7 +11,7 @@ end """ ForwardEulerTimeStepper(grid, prognostic_fields; implicit_solver = nothing, - Gⁿ = map(similar, prognostic_fields)) + Gⁿ = map(deepcopy, prognostic_fields)) Return a 1st-order Forward-Euler (FE) time stepper (`ForwardEulerTimeStepper`) on `grid`, with `tracers`. The tendency fields `Gⁿ`, usually equal to @@ -27,7 +27,7 @@ at the ``n``-th timestep. """ function ForwardEulerTimeStepper(grid, prognostic_fields; implicit_solver::IT = nothing, - Gⁿ = map(similar, prognostic_fields)) where IT + Gⁿ = map(deepcopy, prognostic_fields)) where IT FT = eltype(grid) GT = typeof(Gⁿ) diff --git a/src/sea_ice_advection.jl b/src/sea_ice_advection.jl index 5e7bd67f..6314dd83 100644 --- a/src/sea_ice_advection.jl +++ b/src/sea_ice_advection.jl @@ -1,71 +1,52 @@ using Oceananigans.Operators -using Oceananigans.ImmersedBoundaries +using Oceananigans.ImmersedBoundaries: IBG using Oceananigans.Advection: FluxFormAdvection, _advective_tracer_flux_x, _advective_tracer_flux_y, conditional_flux_fcc, - conditional_flux_cfc + conditional_flux_cfc, + UpwindScheme, + CenteredScheme -# To obtain better numerical properties, the ice thickness is advected together -# with the concentration, i.e.: -# -# A = ℵ⁻¹ ∇ ⋅ (uℵh) -# -# instead of the classical -# -# A = ∇ ⋅ (uh) +using Oceananigans.Advection: _biased_interpolate_xᶠᵃᵃ, _biased_interpolate_yᵃᶠᵃ, bias -_advective_thickness_flux_x(i, j, k, grid, scheme, U, ℵ, h) = advective_thickness_flux_x(i, j, k, grid, scheme, U, ℵ, h) -_advective_thickness_flux_y(i, j, k, grid, scheme, U, ℵ, h) = advective_thickness_flux_y(i, j, k, grid, scheme, U, ℵ, h) +@inline _advective_sea_ice_tracer_flux_x(i, j, k, grid, scheme, u, c) = advective_sea_ice_tracer_flux_x(i, j, k, grid, scheme, u, c) +@inline _advective_sea_ice_tracer_flux_y(i, j, k, grid, scheme, v, c) = advective_sea_ice_tracer_flux_y(i, j, k, grid, scheme, v, c) -_advective_thickness_flux_x(i, j, k, ibg::ImmersedBoundaryGrid, scheme, U, ℵ, h) = - conditional_flux_fcc(i, j, k, ibg, zero(ibg), advective_thickness_flux_x(i, j, k, ibg, scheme, U, ℵ, h)) +@inline _advective_sea_ice_tracer_flux_x(i, j, k, ibg::IBG, scheme, u, c) = conditional_flux_fcc(i, j, k, ibg, zero(ibg), advective_sea_ice_tracer_flux_x(i, j, k, ibg, scheme, u, c)) +@inline _advective_sea_ice_tracer_flux_y(i, j, k, ibg::IBG, scheme, v, c) = conditional_flux_cfc(i, j, k, ibg, zero(ibg), advective_sea_ice_tracer_flux_y(i, j, k, ibg, scheme, v, c)) -_advective_thickness_flux_y(i, j, k, ibg::ImmersedBoundaryGrid, scheme, U, ℵ, h) = - conditional_flux_cfc(i, j, k, ibg, zero(ibg), advective_thickness_flux_y(i, j, k, ibg, scheme, U, ℵ, h)) - -@inline function advective_thickness_flux_x(i, j, k, grid, advection, U, ℵ, h) - ϕℵ = advective_tracer_flux_x(i, j, k, grid, advection, U, ℵ) / Axᶠᶜᶜ(i, j, k, grid) - Uϕℵh = ϕℵ * advective_tracer_flux_x(i, j, k, grid, advection, U, h) - @inbounds ϕℵh = ifelse(U[i, j, k] == 0, zero(grid), Uϕℵh / U[i, j, k]) - return ϕℵh +@inline function advective_sea_ice_tracer_flux_x(i, j, k, grid, scheme::UpwindScheme, u, c) + ũ = ℑyᵃᶜᵃ(i, j, k, grid, Ax_qᶠᶠᶜ, u) + cᴿ = _biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, bias(ũ), c) + return ũ * cᴿ end -@inline function advective_thickness_flux_y(i, j, k, grid, advection, V, ℵ, h) - ϕℵ = advective_tracer_flux_y(i, j, k, grid, advection, V, ℵ) / Ayᶜᶠᶜ(i, j, k, grid) - Vϕℵh = ϕℵ * advective_tracer_flux_y(i, j, k, grid, advection, V, h) - @inbounds ϕℵh = ifelse(V[i, j, k] == 0, zero(grid), Vϕℵh / V[i, j, k]) - return ϕℵh +@inline function advective_sea_ice_tracer_flux_y(i, j, k, grid, scheme::UpwindScheme, v, c) + ṽ = ℑxᶜᵃᵃ(i, j, k, grid, Ay_qᶠᶠᶜ, v) + cᴿ = _biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, bias(ṽ), c) + return ṽ * cᴿ end -@inline div_Uℵh(i, j, k, grid, ::Nothing, U, ℵ, h) = zero(grid) - -# For thickness, we compute [ℵ⁻¹ ∇ ⋅ (uℵh)] -@inline function div_Uℵh(i, j, k, grid, advection, U, ℵ, h) - ∇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)) - - return ℵ⁻¹ * ∇Uℵh +@inline function advective_sea_ice_tracer_flux_x(i, j, k, grid, scheme::CenteredScheme, u, c) + ũ = ℑyᵃᶜᵃ(i, j, k, grid, Ax_qᶠᶠᶜ, u) + cᴿ = _symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, c) + return ũ * cᴿ end -# For thickness, we compute [ℵ⁻¹ ∇ ⋅ (uℵh)] -@inline function div_Uℵh(i, j, k, grid, advection::FluxFormAdvection, U, ℵ, h) - ∇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)) - - return ℵ⁻¹ * ∇Uℵh +@inline function advective_sea_ice_tracer_flux_y(i, j, k, grid, scheme::CenteredScheme, v, c) + ṽ = ℑxᶜᵃᵃ(i, j, k, grid, Ay_qᶠᶠᶜ, v) + cᴿ = _symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, c) + return ṽ * cᴿ end @inline horizontal_div_Uc(i, j, k, grid, ::Nothing, U, c) = zero(grid) -@inline horizontal_div_Uc(i, j, k, grid, advection, U, c) = - 1 / Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_tracer_flux_x, advection, U.u, c) + - δyᵃᶜᵃ(i, j, k, grid, _advective_tracer_flux_y, advection, U.v, c)) + +@inline horizontal_div_Uc(i, j, k, grid, scheme, U, c) = + 1 / Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_sea_ice_tracer_flux_x, scheme, U.u, c) + + δyᵃᶜᵃ(i, j, k, grid, _advective_sea_ice_tracer_flux_y, scheme, U.v, c)) -@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)) +@inline horizontal_div_Uc(i, j, k, grid, scheme::FluxFormAdvection, U, c) = + 1 / Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_sea_ice_tracer_flux_x, scheme.x, U.u, c) + + δyᵃᶜᵃ(i, j, k, grid, _advective_sea_ice_tracer_flux_y, scheme.y, U.v, c)) diff --git a/src/sea_ice_model.jl b/src/sea_ice_model.jl index 95abba47..224060f9 100644 --- a/src/sea_ice_model.jl +++ b/src/sea_ice_model.jl @@ -7,6 +7,7 @@ using Oceananigans: tupleit, tracernames using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions using Oceananigans.Forcings: model_forcing using ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions: flux_summary +using ClimaSeaIce.Rheologies: rheology_prognostic_tracers struct SeaIceModel{GR, TD, D, TS, CL, U, T, IT, IC, ID, CT, STF, A, F, Arch} <: AbstractModel{TS, Arch} architecture :: Arch @@ -31,8 +32,8 @@ struct SeaIceModel{GR, TD, D, TS, CL, U, T, IT, IC, ID, CT, STF, A, F, Arch} <: advection :: A end -assumed_sea_ice_field_location(name) = name === :u ? (Face, Center, Nothing) : - name === :v ? (Center, Face, Nothing) : +assumed_sea_ice_field_location(name) = name === :u ? (Face, Face, Nothing) : + name === :v ? (Face, Face, Nothing) : (Center, Center, Nothing) function SeaIceModel(grid; @@ -55,12 +56,14 @@ function SeaIceModel(grid; ice_consolidation_thickness = field((Center, Center, Nothing), ice_consolidation_thickness, grid) tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example) + tracers = (tracers..., rheology_prognostic_tracers(dynamics)...) # add prognostic tracers # Next, we form a list of default boundary conditions: - field_names = (:u, :v, :h, :ℵ, :S, tracernames(tracers)...) + field_names = tuple(unique((:u, :v, :h, :ℵ, :S, tracernames(tracers)...))...) bc_tuple = Tuple(FieldBoundaryConditions(grid, assumed_sea_ice_field_location(name)) for name in field_names) + default_boundary_conditions = NamedTuple{field_names}(bc_tuple) # Then we merge specified, embedded, and default boundary conditions. Specified boundary conditions @@ -69,8 +72,8 @@ function SeaIceModel(grid; boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, field_names) if isnothing(velocities) - u = Field{Face, Center, Nothing}(grid, boundary_conditions=boundary_conditions.u) - v = Field{Center, Face, Nothing}(grid, boundary_conditions=boundary_conditions.v) + u = Field{Face, Face, Nothing}(grid, boundary_conditions=boundary_conditions.u) + v = Field{Face, Face, Nothing}(grid, boundary_conditions=boundary_conditions.v) velocities = (; u, v) end @@ -82,12 +85,12 @@ function SeaIceModel(grid; ice_density = field((Center, Center, Nothing), ice_density, grid) # Construct prognostic fields if not provided - ice_thickness = Field{Center, Center, Nothing}(grid, boundary_conditions=boundary_conditions.h) + ice_thickness = Field{Center, Center, Nothing}(grid, boundary_conditions=boundary_conditions.h) ice_concentration = Field{Center, Center, Nothing}(grid, boundary_conditions=boundary_conditions.ℵ) # Adding thickness and concentration if not there prognostic_fields = merge(tracers, (; h = ice_thickness, ℵ = ice_concentration)) - prognostic_fields = if ice_salinity isa ConstantField + prognostic_fields = if ice_salinity isa ConstantField prognostic_fields else merge(prognostic_fields, (; S = ice_salinity)) @@ -97,8 +100,8 @@ 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 = ForwardEulerTimeStepper(grid, prognostic_fields) + tracers = merge(tracers, (; S = ice_salinity)) if !isnothing(ice_thermodynamics) if isnothing(top_heat_flux) @@ -140,8 +143,6 @@ end const SIM = SeaIceModel function set!(model::SIM; h=nothing, ℵ=nothing) - grid = model.grid - arch = architecture(model) !isnothing(h) && set!(model.ice_thickness, h) !isnothing(ℵ) && set!(model.ice_concentration, ℵ) diff --git a/src/sea_ice_time_stepping.jl b/src/sea_ice_time_stepping.jl index d041aa26..811c0df8 100644 --- a/src/sea_ice_time_stepping.jl +++ b/src/sea_ice_time_stepping.jl @@ -1,5 +1,5 @@ using Oceananigans.Utils: Time -using Oceananigans.Fields: flattened_unique_values +using Oceananigans.Fields: flattened_unique_values, ConstantField, ZeroField using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series! using Oceananigans.ImmersedBoundaries: mask_immersed_field_xy! @@ -44,7 +44,16 @@ function dynamic_time_step!(model::SIM, Δt) Gⁿ = model.timestepper.Gⁿ - launch!(arch, grid, :xy, _dynamic_step_tracers!, h, ℵ, tracers, Gⁿ, Δt) + launch!(arch, grid, :xy, _dynamic_step_ice_variables!, h, ℵ, tracers, Gⁿ, Δt) + + # Advance tracers + for tracer_idx in keys(tracers) + tracer = @inbounds tracers[tracer_idx] + if (tracer_idx ∈ keys(Gⁿ)) + tendency = @inbounds Gⁿ[tracer_idx] + interior(tracer) .+= Δt .* interior(tendency) + end + end return nothing end @@ -53,13 +62,13 @@ end # We compute hⁿ⁺¹ and ℵⁿ⁺¹ in the same kernel to account for ridging: # if ℵ > 1, we reset the concentration to 1 and adjust the thickness # to conserve the total ice volume in the cell. -@kernel function _dynamic_step_tracers!(h, ℵ, tracers, Gⁿ, Δt) +@kernel function _dynamic_step_ice_variables!(h, ℵ, tracers, Gⁿ, Δt) i, j = @index(Global, NTuple) k = 1 Ghⁿ = Gⁿ.h Gℵⁿ = Gⁿ.ℵ - + # Update ice thickness, clipping negative values @inbounds begin h⁺ = h[i, j, k] + Δt * Ghⁿ[i, j, k] diff --git a/src/tracer_tendency_kernel_functions.jl b/src/tracer_tendency_kernel_functions.jl index fb3fb196..47de96d3 100644 --- a/src/tracer_tendency_kernel_functions.jl +++ b/src/tracer_tendency_kernel_functions.jl @@ -12,25 +12,51 @@ function compute_tracer_tendencies!(model::SIM) arch = architecture(grid) launch!(arch, grid, :xy, - _compute_dynamic_tracer_tendencies!, + _compute_dynamic_ice_variables_tendencies!, model.timestepper.Gⁿ, grid, model.velocities, model.advection, model.ice_thickness, - model.ice_concentration, - model.tracers) + model.ice_concentration) + + # Advance tracers + for tracer_idx in keys(model.tracers) + tracer = @inbounds model.tracers[tracer_idx] + + if (tracer_idx ∈ keys(model.timestepper.Gⁿ)) + tendency = @inbounds model.timestepper.Gⁿ[tracer_idx] + launch!(arch, grid, :xy, + _compute_dynamic_tracer_tendency!, + tendency, + grid, + model.velocities, + model.advection, + tracer) + end + + end return nothing end -@kernel function _compute_dynamic_tracer_tendencies!(Gⁿ, - grid, - velocities, - advection, - ice_thickness, - ice_concentration, - tracers) +@kernel function _compute_dynamic_tracer_tendency!(Gⁿ, + grid, + velocities, + advection, + tracer) + i, j = @index(Global, NTuple) + kᴺ = size(grid, 3) # Assumption! The sea ice is located at the _top_ of the grid + + @inbounds Gⁿ[i, j, 1] = - horizontal_div_Uc(i, j, kᴺ, grid, advection, velocities, tracer) +end + +@kernel function _compute_dynamic_ice_variables_tendencies!(Gⁿ, + grid, + velocities, + advection, + ice_thickness, + ice_concentration) i, j = @index(Global, NTuple) kᴺ = size(grid, 3) # Assumption! The sea ice is located at the _top_ of the grid @@ -39,10 +65,20 @@ end Gⁿ.h[i, j, 1] = - horizontal_div_Uc(i, j, kᴺ, grid, advection, velocities, ice_thickness) Gⁿ.ℵ[i, j, 1] = - horizontal_div_Uc(i, j, kᴺ, grid, advection, velocities, ice_concentration) - # for (n, θ) in enumerate(tracers) - # @inbounds Gⁿ[n] = - horizontal_div_Uc(i, j, 1, grid, advection, velocities, θ) - # end + # TODO: BBM rheology needs this! + # compute_tracer_tendencies!(Gⁿ, i, j, grid, advection, velocities, tracers) end end +const EmptyTuples = Union{NamedTuple{(), Tuple{}}, Tuple{}} +compute_tracer_tendencies!(G, i, j, grid, advection, velocities, ::EmptyTuples) = nothing + +function compute_tracer_tendencies!(G, i, j, grid, advection, velocities, tracers) + # Assumption! The tracer tendencies are the first ones + for n in keys(tracers) + if n ∈ keys(G) + @inbounds G[n][i, j, 1] = - horizontal_div_Uc(i, j, 1, grid, advection, velocities, tracers[n]) + end + end +end diff --git a/test/test_time_stepping.jl b/test/test_time_stepping.jl index b4ff2a6c..440b734b 100644 --- a/test/test_time_stepping.jl +++ b/test/test_time_stepping.jl @@ -21,18 +21,18 @@ end rheologies = (ElastoViscoPlasticRheology(), ViscousRheology(ν=1000)) advections = (WENO(), UpwindBiased(order=5)) - ice_thermodynamics = (nothing, SlabSeaIceThermodynamics(grid)) + thermodynamics = (nothing, SlabSeaIceThermodynamics(grid)) coriolises = (nothing, FPlane(latitude=45), BetaPlane(latitude=45)) solvers = (ExplicitSolver(), SplitExplicitSolver()) - for coriolis in coriolises, advection in advections, rheology in rheologies, ice_thermodynamics in ice_thermodynamics, solver in solvers + for coriolis in coriolises, advection in advections, rheology in rheologies, ice_thermodynamics in thermodynamics, solver in solvers dynamics = SeaIceMomentumEquation(grid; coriolis, rheology, solver) @test time_step_sea_ice_model_works(grid; dynamics, - ice_thermodynamics=ice_thermodynamics, - advection=advection) + ice_thermodynamics, + advection) end end end \ No newline at end of file