From 07faebc01796aa31cbff14a2acacc2b51bc2011d Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Fri, 7 Jul 2023 10:39:47 +0000 Subject: [PATCH] Changes in code API and structs to be compatible with Halfar tests --- Project.toml | 8 ++ src/Huginn.jl | 1 - src/models/iceflow/IceflowModel.jl | 23 +++- src/models/iceflow/SIA2D/SIA2D.jl | 24 ++-- src/models/iceflow/SIA2D/SIA2D_utils.jl | 47 +++---- src/parameters/SolverParameters.jl | 41 ++++++- src/simulations/predictions/Prediction.jl | 24 ++-- .../predictions/prediction_utils.jl | 39 +++--- test/halfar.jl | 116 ++++++++++++++++++ test/runtests.jl | 22 ++-- test/utils_test.jl | 46 +++++++ 11 files changed, 312 insertions(+), 79 deletions(-) create mode 100644 test/halfar.jl create mode 100644 test/utils_test.jl diff --git a/Project.toml b/Project.toml index dfe0043..5173c84 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,14 @@ 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" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] \ No newline at end of file diff --git a/src/Huginn.jl b/src/Huginn.jl index a29fd43..130380b 100644 --- a/src/Huginn.jl +++ b/src/Huginn.jl @@ -17,7 +17,6 @@ import Pkg using Distributed using ProgressMeter using PyCall -using TimerOutputs using Reexport ### ODINN.jl dependencies ### diff --git a/src/models/iceflow/IceflowModel.jl b/src/models/iceflow/IceflowModel.jl index 078811d..1ace06f 100644 --- a/src/models/iceflow/IceflowModel.jl +++ b/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") \ No newline at end of file +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 \ No newline at end of file diff --git a/src/models/iceflow/SIA2D/SIA2D.jl b/src/models/iceflow/SIA2D/SIA2D.jl index 1ac6ed2..25cb50a 100644 --- a/src/models/iceflow/SIA2D/SIA2D.jl +++ b/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) diff --git a/src/models/iceflow/SIA2D/SIA2D_utils.jl b/src/models/iceflow/SIA2D/SIA2D_utils.jl index 1fa7ce9..02c69d7 100644 --- a/src/models/iceflow/SIA2D/SIA2D_utils.jl +++ b/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 diff --git a/src/parameters/SolverParameters.jl b/src/parameters/SolverParameters.jl index 2af90c4..dfdb2d3 100644 --- a/src/parameters/SolverParameters.jl +++ b/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) diff --git a/src/simulations/predictions/Prediction.jl b/src/simulations/predictions/Prediction.jl index 0def63a..906a839 100644 --- a/src/simulations/predictions/Prediction.jl +++ b/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, diff --git a/src/simulations/predictions/prediction_utils.jl b/src/simulations/predictions/prediction_utils.jl index 1a7e2ed..85fa2f9 100644 --- a/src/simulations/predictions/prediction_utils.jl +++ b/src/simulations/predictions/prediction_utils.jl @@ -6,7 +6,7 @@ function run!(simulation::Prediction) println("Running forward PDE ice flow model...\n") results_list = @showprogress pmap((glacier_idx) -> batch_iceflow_PDE(glacier_idx, simulation), 1:length(simulation.glaciers)) - save_results_file!(results_list, simulation) + Sleipnir.save_results_file!(results_list, simulation) @everywhere GC.gc() # run garbage collector @@ -22,19 +22,19 @@ function batch_iceflow_PDE(glacier_idx::Int, simulation::Prediction) model = simulation.model params = simulation.parameters glacier = simulation.glaciers[glacier_idx] - println("Processing glacier: ", glacier.gdir.rgi_id) + glacier_id = isnothing(glacier.gdir) ? "unnamed" : glacier.gdir + println("Processing glacier: ", glacier_id) + # Initialize glacier ice flow model initialize_iceflow_model!(model.iceflow, glacier_idx, glacier, params) params.solver.tstops = define_callback_steps(params.simulation.tspan, params.solver.step) - stop_condition(u,t,integrator) = stop_condition_tstops(u,t,integrator, params.solver.tstops) #closure + stop_condition(u,t,integrator) = Sleipnir.stop_condition_tstops(u,t,integrator, params.solver.tstops) #closure function action!(integrator) if params.simulation.use_MB # Compute mass balance - @timeit get_timer("ODINN") "Mass balance" begin MB_timestep!(model, glacier, params.solver, integrator.t) apply_MB_mask!(integrator.u, glacier, model.iceflow) - end end # # Recompute A value # A = context[1] @@ -45,20 +45,30 @@ function batch_iceflow_PDE(glacier_idx::Int, simulation::Prediction) # Run iceflow PDE for this glacier du = params.simulation.use_iceflow ? SIA2D! : noSIA2D! - results = @timeit get_timer("ODINN") "Simulate iceflow PDE" simulate_iceflow_PDE!(simulation, model, params, cb_MB; du = du) + results = simulate_iceflow_PDE!(simulation, model, params, cb_MB; du = du) return results end """ - simulate_iceflow_PDE!(simulation::Simulation, model::Model, params::Parameters, cb::DiscreteCallback; du = SIA2D!) + function simulate_iceflow_PDE!( + simulation::SIM, + model::Sleipnir.Model, + params::Sleipnir.Parameters, + cb::DiscreteCallback; + du = SIA2D!) where {SIM <: Simulation} Make forward simulation of the iceflow PDE determined in `du`. """ -function simulate_iceflow_PDE!(simulation::SIM, model::Model, params::Parameters, cb::DiscreteCallback; du = SIA2D!) where {SIM <: Simulation} +function simulate_iceflow_PDE!( + simulation::SIM, + model::Sleipnir.Model, + params::Sleipnir.Parameters, + cb::DiscreteCallback; + du = SIA2D!) where {SIM <: Simulation} + # Define problem to be solved iceflow_prob = ODEProblem{true,SciMLBase.FullSpecialize}(du, model.iceflow.H, params.simulation.tspan, tstops=params.solver.tstops, simulation) - @timeit get_timer("ODINN") "Solve" begin iceflow_sol = solve(iceflow_prob, params.solver.solver, callback=cb, @@ -67,23 +77,18 @@ function simulate_iceflow_PDE!(simulation::SIM, model::Model, params::Parameters save_everystep=params.solver.save_everystep, progress=params.solver.progress, progress_steps=params.solver.progress_steps) - end # @show iceflow_sol.destats # Compute average ice surface velocities for the simulated period - @timeit get_timer("ODINN") "Postprocessing" begin model.iceflow.H .= iceflow_sol.u[end] map!(x -> ifelse(x>0.0,x,0.0), model.iceflow.H, model.iceflow.H) - @timeit get_timer("ODINN") "Surface v" avg_surface_V!(iceflow_sol[begin], simulation) # Average velocity with average temperature + avg_surface_V!(iceflow_sol[begin], simulation) # Average velocity with average temperature glacier_idx = simulation.model.iceflow.glacier_idx - glacier::Glacier = simulation.glaciers[glacier_idx[]] + glacier::Sleipnir.Glacier2D = simulation.glaciers[glacier_idx[]] model.iceflow.S .= glacier.B .+ model.iceflow.H # Surface topography - end # Update simulation results - @timeit get_timer("ODINN") "Results" begin - results = create_results(simulation, glacier_idx[], iceflow_sol; light=true) - end + results = Sleipnir.create_results(simulation, glacier_idx[], iceflow_sol; light=true) return results end diff --git a/test/halfar.jl b/test/halfar.jl new file mode 100644 index 0000000..56b2dc5 --- /dev/null +++ b/test/halfar.jl @@ -0,0 +1,116 @@ + +using Revise +using Infiltrator +using Test +using Plots +using Reexport + +# @reexport using Sleipnir +using Huginn + +# Include utils for halfar +include("./utils_test.jl") + +""" + unit_halfar_test(; A, t₀, t₁, Δx, Δy, nx, ny, h₀, r₀, rtol=0.02, atol=1.0, distance_to_border=3, save_plot=false) + +Test one single iteration of iceflow with the Halfar solution + +Arguments +================= + - `A`: Glen coefficient + - `t₀`: Initial time in simulation (must be different than zero!) + - `t₁`: Final time in simulation + - `Δx`, `Δy`: Spacial timesteps + - `nx`, `ny`: Number of grid points + - `h₀`, `r₀`: Parameters in the Halfar solutions + - `rtol`: Relative tolerance for the test + - `atol`: Absolute tolerance for the test + - `distance_to_border`: Minimum distance to the border used to evaluate test. Points close to the border are not considered. + - `save_plot`: Save plot with comparision of prediction and true solution. +""" +function unit_halfar_test(; A, t₀, t₁, Δx, Δy, nx, ny, h₀, r₀, rtol=0.02, atol=1.0, distance_to_border=3, save_plot=false) + + # Get parameters for a simulation + parameters = Parameters(simulation=SimulationParameters(tspan=(t₀, t₁), + use_MB=false, + use_iceflow=true, + multiprocessing=true, + workers=1), + physical=PhysicalParameters(A=A), + solver=SolverParameters(reltol=1e-12)) + + # Bed (it has to be flat for the Halfar solution) + B = zeros((nx,ny)) + + model = Model(iceflow = SIA2Dmodel(parameters)) #, + # mass_balance = nothing, # TImodel1(parameters; DDF=8.0/1000.0, acc_factor=1.0/1000.0), + # machine_learning = nothing) # NN(parameters)) + + # Initial condition of the glacier + R₀ = [sqrt((Δx * (i - nx/2))^2 + (Δy * (j - ny/2))^2) for i in 1:nx, j in 1:ny] + H₀ = halfar_solution(R₀, t₀, h₀, r₀, parameters.physical) + S = B + H₀ + # Final expected solution + H₁ = halfar_solution(R₀, t₁, h₀, r₀, parameters.physical) + + # Define glacier object + glacier = Glacier(rgi_id = "toy", H₀ = H₀, S = S, B = B, + Δx=Δx, Δy=Δy, nx=nx, ny=ny) + glaciers = Vector{Sleipnir.AbstractGlacier}([glacier]) + + prediction = Prediction(model, glaciers, parameters) + run!(prediction) + + # Final solution + H₁_pred = prediction.results[1].H[end] + H_diff = H₁_pred .- H₁ + + # Error calculation + absolute_error = maximum(abs.(H_diff[is_border(H₁, distance_to_border)])) + percentage_error = maximum(abs.((H_diff./H₁)[is_border(H₁, distance_to_border)])) + maximum_flow = maximum(abs.(((H₁ .- H₀))[is_border(H₁, distance_to_border)])) + + # Optional plot + if save_plot + fig = Figure(resolution = (800, 800)) + + Axis(fig[1, 1], title = "Initial Condition") + heatmap!(H₀, colormap=:viridis, colorrange=(0, maximum(H₀))) + + Axis(fig[1, 2], title = "Final State") + heatmap!(H₁, colormap=:viridis, colorrange=(0, maximum(H₀))) + + Axis(fig[2,1], title="Prediction") + heatmap!(H₁_pred, colormap=:viridis, colorrange=(0, maximum(H₀))) + + Axis(fig[2,2], title="Difference") + heatmap!(H_diff, colormap=Reverse(:balance), colorrange=(-10, 10)) + + save("test/halfar_test.png", fig) + end + + # @show percentage_error, absolute_error, maximum_flow + @test all([percentage_error < rtol, absolute_error < atol]) +end + +""" + function halfar_test(; rtol, atol) +Multiple tests using Halfar solution. + +Arguments +================= + - `rtol`: Relative tolerance for ice thickness H + - `atol`: Absolute tolerance for ice thickness H +""" +function halfar_test(; rtol, atol) + unit_halfar_test(A=4e-17, t₀=5.0, t₁=10.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=1000, rtol=rtol, atol=atol) + unit_halfar_test(A=8e-17, t₀=5.0, t₁=10.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=1000, rtol=rtol, atol=atol) + unit_halfar_test(A=4e-17, t₀=5.0, t₁=40.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol) + unit_halfar_test(A=8e-17, t₀=5.0, t₁=40.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol) + unit_halfar_test(A=4e-17, t₀=5.0, t₁=100.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol) + unit_halfar_test(A=8e-17, t₀=5.0, t₁=100.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol) + unit_halfar_test(A=4e-17, t₀=5.0, t₁=40.0, Δx=80.0, Δy=80.0, nx=100, ny=100, h₀=300, r₀=1000, rtol=rtol, atol=atol) + unit_halfar_test(A=8e-17, t₀=5.0, t₁=40.0, Δx=80.0, Δy=80.0, nx=100, ny=100, h₀=300, r₀=1000, rtol=rtol, atol=atol) + unit_halfar_test(A=4e-17, t₀=5.0, t₁=10.0, Δx=10.0, Δy=10.0, nx=500, ny=500, h₀=300, r₀=1000, rtol=rtol, atol=atol) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 6591b24..f0e289b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,23 +1,27 @@ import Pkg Pkg.activate(dirname(Base.current_project())) -using PyCall -# Update SSL certificate to avoid issue in GitHub Actions CI -certifi = pyimport("certifi") -ENV["SSL_CERT_FILE"] = certifi.where() -# println("Current SSL certificate: ", ENV["SSL_CERT_FILE"]) - using Revise using Test using JLD2 +using Plots using Infiltrator -include(joinpath("..", "src/Huginn.jl")) +using Huginn -include(joinpath(Sleipnir.root_dir, "test/params_construction.jl")) +# include("PDE_UDE_solve.jl") +include("halfar.jl") +# include("mass_conservation.jl") # Activate to avoid GKS backend Plot issues in the JupyterHub ENV["GKSwstype"]="nul" -@testset "Solver parameters constructors" params_construction() +atol = 0.01 +# @testset "PDE and UDE SIA solvers without MB" pde_solve_test(atol; MB=false, fast=true) + +atol = 2.0 +# @testset "PDE and UDE SIA solvers with MB" pde_solve_test(atol; MB=true, fast=true) + +@testset "Halfar Solution" halfar_test(rtol=0.02, atol=1.0) +# @testset "Conservation of Mass - Flat Bed" unit_mass_flatbed_test(rtol=1.0e-7, atol=1000) \ No newline at end of file diff --git a/test/utils_test.jl b/test/utils_test.jl new file mode 100644 index 0000000..7c882e8 --- /dev/null +++ b/test/utils_test.jl @@ -0,0 +1,46 @@ +# -*- coding: utf-8 -*- +""" + is_border(A::Matrix, distance::Int) + +Return a matrix with booleans indicating if a given pixel is at distance at most `distance` of the end of the +matrix, which is indicated by the first pixel value to reach zero. + +Arguments: + - A: Array + - distance: distance to the border, computed as the number of pixels we need to move to find a pixel with value zero +""" +function is_border(A::Matrix{F}, distance::Int) where {F <: AbstractFloat} + B = copy(A) + for i in 1:distance + B .= min.(B, circshift(B, (1,0)), circshift(B, (-1,0)), circshift(B, (0,1)), circshift(B, (0,-1))) + end + return B .> 0.001 +end + +""" + halfar_solution(t, r, θ) + +Returns the evaluation of the Halfar solutions for the SIA equation. + +Arguments: + - r: radial distance. The solutions have polar symmetry around the center of origin. + - t: time + - ν = (A, H₀, R₀) +""" +function halfar_solution(R, t, h₀, r₀, physical_parameters::PhysicalParameters) + + # parameters of Halfar solutions + ρ = physical_parameters.ρ + g = physical_parameters.g + n = physical_parameters.n + A = physical_parameters.A + + Γ = 2 * A * (ρ * g)^n / (n+2) + # Characteristic time + τ₀ = (7/4)^3 * r₀^4 / ( 18 * Γ * h₀^7 ) + + return [r <= r₀ * (t/τ₀)^(1/18) ? + h₀ * (τ₀/t)^(1/9) * ( 1 - ( (τ₀/t)^(1/18) * (r/r₀) )^(4/3) )^(3/7) : + 0.0 + for r in R] +end