Changes for forward/reverse compatibility

Project.toml

@@ -29,15 +29,15 @@ BenchmarkTools = "1"
 CairoMakie = "0.11"
 Infiltrator = "1"
 JLD2 = "0.4"
-Muninn = "0.2.7"
+Muninn = "0.3"
 OrdinaryDiffEq = "6"
 PlotThemes = "3"
 Plots = "1"
 ProgressMeter = "1"
 PyCall = "1"
 Reexport = "1"
 Revise = "3"
-Sleipnir = "0.5"
+Sleipnir = "0.6"
 Tullio = "0.3"
 julia = "1.7"
 

src/models/iceflow/SIA2D/SIA2D.jl

@@ -7,7 +7,7 @@ include("SIA2D_utils.jl")
 ###### SHALLOW ICE APPROXIMATION MODELS #######
 ###############################################
 
-mutable struct SIA2Dmodel{F <: AbstractFloat, I <: Integer, R <: Real} <: SIAmodel
+mutable struct SIA2Dmodel{R <: Real, I <: Integer} <: SIAmodel
     A::Union{Ref{R}, Nothing}
     n::Union{Ref{R}, Nothing}
     C::Union{Ref{R}, Matrix{R}, Nothing}
@@ -73,22 +73,22 @@ function SIA2Dmodel(params::Sleipnir.Parameters;
     ft = params.simulation.float_type
     it = params.simulation.int_type
     if !isnothing(A)
-        A = Ref{R}(A)
+        A = Ref{ft}(A)
     end
     if !isnothing(n)
-        n = Ref{R}(n)
+        n = Ref{ft}(n)
     end
     if !isnothing(C)
-        C = Ref{R}(C)
+        C = Ref{ft}(C)
     end
     if !isnothing(Γ)
-        Γ = Ref{R}(Γ)
+        Γ = Ref{ft}(Γ)
     end
     if !isnothing(glacier_idx)
         glacier_idx = Ref{I}(glacier_idx)
     end
 
-    SIA2D_model = SIA2Dmodel{ft, it, Real}(A, n, C, H₀, H, H̄, S, dSdx, dSdy, D, Dx, Dy, dSdx_edges, dSdy_edges,
+    SIA2D_model = SIA2Dmodel{ft, it}(A, n, C, H₀, H, H̄, S, dSdx, dSdy, D, Dx, Dy, dSdx_edges, dSdy_edges,
                             ∇S, ∇Sx, ∇Sy, Fx, Fy, Fxx, Fyy, V, Vx, Vy, Γ, MB, MB_mask, MB_total, glacier_idx)
 
     return SIA2D_model
@@ -117,9 +117,9 @@ function initialize_iceflow_model!(iceflow_model::IF,
                                    ) where {IF <: IceflowModel, I <: Integer, G <: Sleipnir.AbstractGlacier}
     nx, ny = glacier.nx, glacier.ny
     F = params.simulation.float_type
-    iceflow_model.A = isnothing(iceflow_model.A) ? Ref{Real}(glacier.A) : iceflow_model.A
-    iceflow_model.n = isnothing(iceflow_model.n) ? Ref{Real}(glacier.n) : iceflow_model.n
-    iceflow_model.C = isnothing(iceflow_model.C) ? Ref{Real}(glacier.C) : iceflow_model.C
+    iceflow_model.A = isnothing(iceflow_model.A) ? Ref{F}(glacier.A) : iceflow_model.A
+    iceflow_model.n = isnothing(iceflow_model.n) ? Ref{F}(glacier.n) : iceflow_model.n
+    iceflow_model.C = isnothing(iceflow_model.C) ? Ref{F}(glacier.C) : iceflow_model.C
     iceflow_model.H₀ = deepcopy(glacier.H₀)
     iceflow_model.H = deepcopy(glacier.H₀)
     iceflow_model.H̄ = zeros(F,nx-1,ny-1)
@@ -141,7 +141,7 @@ function initialize_iceflow_model!(iceflow_model::IF,
     iceflow_model.V = zeros(F,nx,ny)
     iceflow_model.Vx = zeros(F,nx,ny)
     iceflow_model.Vy = zeros(F,nx,ny)
-    iceflow_model.Γ = isnothing(iceflow_model.Γ) ? Ref{Real}(0.0) : iceflow_model.Γ
+    iceflow_model.Γ = isnothing(iceflow_model.Γ) ? Ref{F}(0.0) : iceflow_model.Γ
     iceflow_model.MB = zeros(F,nx,ny)
     iceflow_model.MB_mask= zeros(F,nx,ny)
     iceflow_model.MB_total = zeros(F,nx,ny)
@@ -172,8 +172,8 @@ function initialize_iceflow_model(iceflow_model::IF,
     # nx, ny = glacier.nx, glacier.ny
     F = params.simulation.float_type
     nx, ny = glacier.nx, glacier.ny
-    iceflow_model.A = Ref{Real}(glacier.A)
-    iceflow_model.n = Ref{Real}(glacier.n)
+    iceflow_model.A = Ref{F}(glacier.A)
+    iceflow_model.n = Ref{F}(glacier.n)
     iceflow_model.glacier_idx = Ref{I}(glacier_idx)
     iceflow_model.H₀ = deepcopy(glacier.H₀)
     iceflow_model.H  = deepcopy(glacier.H₀)

src/models/iceflow/SIA2D/SIA2D_utils.jl

@@ -5,8 +5,7 @@ SIA2D!(dH, H, SIA2Dmodel)
 
 Compute an in-place step of the Shallow Ice Approximation PDE in a forward model
 """
-function SIA2D!(dH::Matrix{<:Real}, H::Matrix{<:Real}, simulation::SIM, t::F) where {F <: AbstractFloat, SIM <: Simulation}
-
+function SIA2D!(dH::Matrix{R}, H::Matrix{R}, simulation::SIM, t::R) where {R <:Real, SIM <: Simulation}
     # Retrieve parameters
     SIA2D_model::SIA2Dmodel = simulation.model.iceflow
     glacier::Sleipnir.Glacier2D = simulation.glaciers[simulation.model.iceflow.glacier_idx[]]
@@ -160,8 +159,8 @@ function avg_surface_V!(simulation::SIM) where {SIM <: Simulation}
     Vx₀, Vy₀ = surface_V!(iceflow_model.H₀, simulation)
     Vx,  Vy  = surface_V!(iceflow_model.H,  simulation)
 
-    iceflow_model.Vx .= (Vx₀ .+ Vx)./2.0
-    iceflow_model.Vy .= (Vy₀ .+ Vy)./2.0
+    inn1(iceflow_model.Vx) .= (Vx₀ .+ Vx)./2.0
+    inn1(iceflow_model.Vy) .= (Vy₀ .+ Vy)./2.0
     iceflow_model.V  .= (iceflow_model.Vx.^2 .+ iceflow_model.Vy.^2).^(1/2) 
 end
 

src/simulations/predictions/prediction_utils.jl

@@ -70,6 +70,7 @@ function simulate_iceflow_PDE!(
 
     # Define problem to be solved
     iceflow_prob = ODEProblem{true,SciMLBase.FullSpecialize}(du, model.iceflow.H, params.simulation.tspan, tstops=params.solver.tstops, simulation)
+
     iceflow_sol = solve(iceflow_prob, 
                         params.solver.solver, 
                         callback=cb, 
@@ -148,7 +149,7 @@ function batch_iceflow_PDE(glacier_idx::I, simulation::Prediction) where {I <: I
     cb_MB = DiscreteCallback(stop_condition, action!)
 
     # Run iceflow PDE for this glacier
-    du = params.simulation.use_iceflow ? SIA2D! : noSIA2D!
+    du = params.simulation.use_iceflow ? SIA2D : noSIA2D
     results = simulate_iceflow_PDE(simulation, model, params, cb_MB; du = du)
 
     return results
@@ -187,12 +188,12 @@ function simulate_iceflow_PDE(
     map!(x -> ifelse(x>0.0,x,0.0), model.iceflow.H, model.iceflow.H)
 
     # Average surface velocity
-    Vx, Vy, V = avg_surface_V( model.iceflow.H, simulation)
+    avg_surface_V(simulation)
 
     # Since we are doing out-of-place, we need to add this to the result
-    model.iceflow.Vx = Vx
-    model.iceflow.Vy = Vy
-    model.iceflow.V  = V
+    # model.iceflow.Vx = Vx
+    # model.iceflow.Vy = Vy
+    # model.iceflow.V  = V
 
     glacier_idx = simulation.model.iceflow.glacier_idx
     glacier::Sleipnir.Glacier2D = simulation.glaciers[glacier_idx[]]

test/halfar.jl

@@ -21,6 +21,7 @@ function unit_halfar_test(; A, n, t₀, t₁, Δx, Δy, nx, ny, h₀, r₀, rtol
 
     # Get parameters for a simulation 
     parameters = Parameters(simulation=SimulationParameters(tspan=(t₀, t₁),
+                                                            multiprocessing=false,
                                                             use_MB=false,
                                                             use_iceflow=true,
                                                             working_dir=Huginn.root_dir),