Merge pull request #5 from facusapienza21/main

Changes in code API and structs to be compatible with Halfar tests

Project.toml

@@ -14,9 +14,17 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Sleipnir = "f5e6c550-199f-11ee-3608-394420200519"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
+
 [compat]
 OrdinaryDiffEq = "6"
 ProgressMeter = "1"

src/Huginn.jl

@@ -17,7 +17,6 @@ import Pkg
 using Distributed
 using ProgressMeter
 using PyCall
-using TimerOutputs
 using Reexport
 
 ### ODINN.jl dependencies  ###

src/models/iceflow/IceflowModel.jl

@@ -1,9 +1,26 @@
-
+export Model
 
 # Abstract type as a parent type for ice flow models
-abstract type IceflowModel end
+abstract type IceflowModel <: AbstractModel end
 
 # Subtype structure for Shallow Ice Approximation models
 abstract type SIAmodel <: IceflowModel end
 
-include("SIA2D/SIA2D.jl")
+include("iceflow_utils.jl")
+include("SIA2D/SIA2D.jl")
+
+"""
+function Model(;
+    iceflow::IceflowModel
+    )
+Initialize Huginn flow model
+
+"""
+function Model(;
+    iceflow::IceflowModel
+    )
+
+    model = Sleipnir.Model(iceflow, nothing, nothing)
+
+    return model
+end

src/models/iceflow/SIA2D/SIA2D.jl

@@ -1,6 +1,8 @@
 
 export SIA2Dmodel
 
+include("SIA2D_utils.jl")
+
 ###############################################
 ###### SHALLOW ICE APPROXIMATION MODELS #######
 ###############################################
@@ -34,7 +36,7 @@ mutable struct SIA2Dmodel{F <: AbstractFloat, I <: Integer} <: SIAmodel
     glacier_idx::Union{Ref{I}, Nothing}
 end
 
-function SIA2Dmodel(params::Parameters; 
+function SIA2Dmodel(params::Sleipnir.Parameters;
                     A::Union{Ref{F}, Nothing} = nothing,
                     H::Union{Matrix{F}, Nothing} = nothing,
                     H̄::Union{Matrix{F}, Nothing} = nothing,
@@ -71,22 +73,26 @@ function SIA2Dmodel(params::Parameters;
 end
 
 """
-    initialize_iceflow_model!(; glacier::Glacier,
-                                params::Parameters
-                                ) where F <: AbstractFloat
+function initialize_iceflow_model!(iceflow_model::IF,  
+    glacier_idx::I,
+    glacier::AbstractGlacier,
+    params::Sleipnir.Parameters
+    ) where {IF <: IceflowModel, I <: Int}
 
 Initialize iceflow model data structures to enable in-place mutation.
 
 Keyword arguments
 =================
+    - `iceflow_model`: Iceflow model used for simulation. 
+    - `glacier_idx`: Index of glacier.
     - `glacier`: `Glacier` to provide basic initial state of the ice flow model.
-    - `parameters`: `Parameters` to configure some physical variables
+    - `parameters`: `Parameters` to configure some physical variables.
 """
 function initialize_iceflow_model!(iceflow_model::IF,  
-                                    glacier_idx::I,
-                                    glacier::Glacier,
-                                     params::Parameters
-                                     ) where {IF <: IceflowModel, I <: Int}
+                                   glacier_idx::I,
+                                   glacier::Sleipnir.AbstractGlacier,
+                                   params::Sleipnir.Parameters
+                                   ) where {IF <: IceflowModel, I <: Int}
     nx, ny = glacier.nx, glacier.ny
     F = params.simulation.float_type
     iceflow_model.A = Ref{F}(params.physical.A)

src/models/iceflow/SIA2D/SIA2D_utils.jl

@@ -8,11 +8,10 @@ Compute an in-place step of the Shallow Ice Approximation PDE in a forward model
 function SIA2D!(dH::Matrix{F}, H::Matrix{F}, simulation::SIM, t::F) where {F <: AbstractFloat, SIM <: Simulation}
 
     # Retrieve parameters
-    @timeit get_timer("ODINN") "Variable assignment" begin
     SIA2D_model::SIA2Dmodel = simulation.model.iceflow
-    glacier::Glacier = simulation.glaciers[simulation.model.iceflow.glacier_idx[]]
-    params::Parameters = simulation.parameters
-    int_type = simulation.parameters.simulation.int_type
+    glacier::Sleipnir.Glacier2D = simulation.glaciers[simulation.model.iceflow.glacier_idx[]]
+    params::Sleipnir.Parameters = simulation.parameters
+        int_type = simulation.parameters.simulation.int_type
     H̄::Matrix{F} = SIA2D_model.H̄
     A::Ref{F} = SIA2D_model.A
     B::Matrix{F} = glacier.B
@@ -37,60 +36,46 @@ function SIA2D!(dH::Matrix{F}, H::Matrix{F}, simulation::SIM, t::F) where {F <:
     n::int_type = simulation.parameters.physical.n
     ρ::F = simulation.parameters.physical.ρ
     g::F = simulation.parameters.physical.g
-    end
 
-    @timeit get_timer("ODINN") "H to zero" begin
     # First, enforce values to be positive
     map!(x -> ifelse(x>0.0,x,0.0), H, H)
     # Update glacier surface altimetry
     S .= B .+ H
-    end
 
     # All grid variables computed in a staggered grid
     # Compute surface gradients on edges
-    @timeit get_timer("ODINN") "Surface gradients" begin
     diff_x!(dSdx, S, Δx)  
     diff_y!(dSdy, S, Δy) 
     avg_y!(∇Sx, dSdx)
     avg_x!(∇Sy, dSdy)
     ∇S .= (∇Sx.^2 .+ ∇Sy.^2).^((n - 1)/2) 
-    end
 
-    @timeit get_timer("ODINN") "Diffusivity" begin
     # @infiltrate
     avg!(H̄, H)
     Γ[] = 2.0 * A[] * (ρ * g)^n / (n+2) # 1 / m^3 s 
     D .= Γ[] .* H̄.^(n + 2) .* ∇S
-    end
 
     # Compute flux components
-    @timeit get_timer("ODINN") "Gradients edges" begin
     @views diff_x!(dSdx_edges, S[:,2:end - 1], Δx)
     @views diff_y!(dSdy_edges, S[2:end - 1,:], Δy)
-    end
+
     # Cap surface elevaton differences with the upstream ice thickness to 
     # imporse boundary condition of the SIA equation
-    @timeit get_timer("ODINN") "Capping flux" begin
     η₀ = params.physical.η₀
-    dSdx_edges .= @views @. min(dSdx_edges,  η₀ * H[1:end-1, 2:end-1]/Δx,  η₀ * H[2:end, 2:end-1]/Δx)
-    dSdy_edges .= @views @. min(dSdy_edges,  η₀ * H[2:end-1, 1:end-1]/Δy,  η₀ * H[2:end-1, 2:end]/Δy) 
-    dSdx_edges .= @views @. max(dSdx_edges, -η₀ * H[1:end-1, 2:end-1]/Δx, -η₀ * H[2:end, 2:end-1]/Δx)
-    dSdy_edges .= @views @. max(dSdy_edges, -η₀ * H[2:end-1, 1:end-1]/Δy, -η₀ * H[2:end-1, 2:end]/Δy)
-    end
+    dSdx_edges .= @views @. min(dSdx_edges,  η₀ * H[2:end, 2:end-1]/Δx)
+    dSdx_edges .= @views @. max(dSdx_edges, -η₀ * H[1:end-1, 2:end-1]/Δx)
+    dSdy_edges .= @views @. min(dSdy_edges,  η₀ * H[2:end-1, 2:end]/Δy)
+    dSdy_edges .= @views @. max(dSdy_edges, -η₀ * H[2:end-1, 1:end-1]/Δy)
 
-    @timeit get_timer("ODINN") "Flux" begin
     avg_y!(Dx, D)
     avg_x!(Dy, D)
     Fx .= .-Dx .* dSdx_edges
     Fy .= .-Dy .* dSdy_edges 
-    end
 
     #  Flux divergence
-    @timeit get_timer("ODINN") "dH" begin
     diff_x!(Fxx, Fx, Δx)
     diff_y!(Fyy, Fy, Δy)
     inn(dH) .= .-(Fxx .+ Fyy) 
-    end
 end
 
 # Dummy function to bypass ice flow
@@ -130,10 +115,14 @@ function SIA2D(H, SIA2Dmodel)
     # imporse boundary condition of the SIA equation
     # We need to do this with Tullio or something else that allow us to set indices.
     η₀ = 1.0
-    dSdx_edges .= min.(dSdx_edges,  η₀ * H[1:end-1, 2:end-1]./Δx,  η₀ * H[2:end, 2:end-1]./Δx)
-    dSdy_edges .= min.(dSdy_edges,  η₀ * H[2:end-1, 1:end-1]./Δy,  η₀ * H[2:end-1, 2:end]./Δy) 
-    dSdx_edges .= max.(dSdx_edges, -η₀ * H[1:end-1, 2:end-1]./Δx, -η₀ * H[2:end, 2:end-1]./Δx)
-    dSdy_edges .= max.(dSdy_edges, -η₀ * H[2:end-1, 1:end-1]./Δy, -η₀ * H[2:end-1, 2:end]./Δy)
+    dSdx_edges .= @views @. min(dSdx_edges,  η₀ * H[2:end, 2:end-1]/Δx)
+    dSdx_edges .= @views @. max(dSdx_edges, -η₀ * H[1:end-1, 2:end-1]/Δx)
+    dSdy_edges .= @views @. min(dSdy_edges,  η₀ * H[2:end-1, 2:end]/Δy)
+    dSdy_edges .= @views @. max(dSdy_edges, -η₀ * H[2:end-1, 1:end-1]/Δy)
+    # dSdx_edges .= min.(dSdx_edges,  η₀ * H[1:end-1, 2:end-1]./Δx,  η₀ * H[2:end, 2:end-1]./Δx)
+    # dSdy_edges .= min.(dSdy_edges,  η₀ * H[2:end-1, 1:end-1]./Δy,  η₀ * H[2:end-1, 2:end]./Δy) 
+    # dSdx_edges .= max.(dSdx_edges, -η₀ * H[1:end-1, 2:end-1]./Δx, -η₀ * H[2:end, 2:end-1]./Δx)
+    # dSdy_edges .= max.(dSdy_edges, -η₀ * H[2:end-1, 1:end-1]./Δy, -η₀ * H[2:end-1, 2:end]./Δy)
 
     Fx = .-avg_y(D) .* dSdx_edges
     Fy = .-avg_x(D) .* dSdy_edges 
@@ -190,11 +179,11 @@ end
 Computes the ice surface velocity for a given glacier state
 """
 function surface_V!(H::Matrix{F}, simulation::SIM) where {F <: AbstractFloat, SIM <: Simulation}
-    params::Parameters = simulation.parameters
+    params::Sleipnir.Parameters = simulation.parameters
     ft = params.simulation.float_type
     it = params.simulation.int_type
     iceflow_model = simulation.model.iceflow
-    glacier::Glacier = simulation.glaciers[iceflow_model.glacier_idx[]]
+    glacier::Sleipnir.Glacier2D = simulation.glaciers[iceflow_model.glacier_idx[]]
     B::Matrix{ft} = glacier.B
     H̄::Matrix{F} = iceflow_model.H̄
     dSdx::Matrix{ft} = iceflow_model.dSdx

src/parameters/SolverParameters.jl

@@ -1,7 +1,9 @@
+export SolverParameters, Parameters
 
-  mutable struct SolverParameters{F <: AbstractFloat, I <: Integer}
+mutable struct SolverParameters{F <: AbstractFloat, I <: Integer} <: AbstractParameters
     solver::OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm
     reltol::F
+    step::F
     tstops::Union{Nothing,Vector{F}} 
     save_everystep::Bool
     progress::Bool
@@ -22,19 +24,50 @@ Keyword arguments
 function SolverParameters(;
             solver::OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm = RDPK3Sp35(),
             reltol::F = 1e-12,
-            save_everystep = false,
+            step::F = 1.0/12.0,
             tstops::Union{Nothing,Vector{F}} = nothing,
+            save_everystep = false,
             progress::Bool = true,
             progress_steps::I = 10
             ) where {F <: AbstractFloat, I <: Integer}
     # Build the solver parameters based on input values
     solver_parameters = SolverParameters(solver, reltol, 
-                                        tstops,
-                                        save_everystep, progress, progress_steps)
+                                         step, tstops,
+                                         save_everystep, progress, progress_steps)
 
     return solver_parameters
 end
 
+"""
+Parameters(;
+        simulation::SimulationParameters = SimulationParameters()
+        physical::PhysicalParameters = PhysicalParameters()
+        OGGM::OGGMparameters = OGGMparameters(),
+        solver::SolverParameters = SolverParameters()
+        )
+Initialize Huginn parameters
+
+Keyword arguments
+=================
+    
+"""
+
+function Parameters(;
+    physical::PhysicalParameters = PhysicalParameters(),
+    simulation::SimulationParameters = SimulationParameters(),
+    OGGM::OGGMparameters = OGGMparameters(),
+    solver::SolverParameters = SolverParameters()
+    ) 
+
+    # Build the parameters based on all the subtypes of parameters
+    parameters = Sleipnir.Parameters(physical, simulation, OGGM,
+                                     nothing, solver, nothing)
+
+    return parameters
+end
+
+# Base.:(==)(a::Parameters, b::Parameters) = a.physical == b.physical && a.simulation == b.simulation && a.OGGM == b.OGGM && a.solver == b.solver
+
 """
     define_callback_steps(tspan::Tuple{Float64, Float64}, step::Float64)
 

src/simulations/predictions/Prediction.jl

@@ -3,17 +3,27 @@ export Prediction
 
 # Subtype composite type for a prediction simulation
 mutable struct Prediction  <: Simulation 
-    model::Model
-    glaciers::Vector{Glacier}
-    parameters::Parameters
+    model::Sleipnir.Model
+    glaciers::Vector{Sleipnir.AbstractGlacier}
+    parameters::Sleipnir.Parameters
     results::Vector{Results}
 end
 
-function Prediction(
-    model::Model,
-    glaciers::Vector{Glacier},
-    parameters::Parameters
+"""
+    function Prediction(
+        model::Sleipnir.Model,
+        glaciers::Vector{Sleipnir.AbstractGlacier},
+        parameters::Sleipnir.Parameters
         )
+Construnctor for Prediction struct with glacier model infomation, glaciers and parameters.
+Keyword arguments
+=================
+"""
+function Prediction(
+    model::Sleipnir.Model,
+    glaciers::Vector{Sleipnir.AbstractGlacier},
+    parameters::Sleipnir.Parameters
+    )
 
     # Build the results struct based on input values
     prediction = Prediction(model,