Skip to content

Commit

Permalink
Made out place SIA2D compatibile with AD (#48)
Browse files Browse the repository at this point in the history 
* Plotting utility fix

* Changed type compatibility of SIA2D constructor

* Compatibility with optimizing

* Updated effiency flow param plotting function

* Removed type specifications in SIA2D utilities

* Added enforce positive in out of place SIA2D

* Added C to iceflow model + reduced Real type

* Added utility for calculating H from V

* Added basal sliding as coefficient

* Added sliding as a function of basal shear stress

* Removed complex C

* small test fixes

* update Sleipnir to v0.5

* update Muninn to 0.2.7

* fixed issue with parameters

* Added surface to outplace iceflow model
  • Loading branch information
@vivekag7
vivekag7 committed Mar 28, 2024
1 parent 95f0ce0 commit 6ead632
Show file tree
Hide file tree
Showing 9 changed files with 175 additions and 141 deletions.
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,15 +27,15 @@ BenchmarkTools = "1"
CairoMakie = "0.11"
Infiltrator = "1"
JLD2 = "0.4"
Muninn = "0.2"
Muninn = "0.2.7"
OrdinaryDiffEq = "6"
PlotThemes = "3"
Plots = "1"
ProgressMeter = "1"
PyCall = "1"
Reexport = "1"
Revise = "3"
Sleipnir = "0.4"
Sleipnir = "0.5"
Tullio = "0.3"
julia = "1.7"

Expand Down
48 changes: 31 additions & 17 deletions src/models/iceflow/SIA2D/SIA2D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,9 @@ include("SIA2D_utils.jl")
###############################################

mutable struct SIA2Dmodel{F <: AbstractFloat, I <: Integer} <: SIAmodel
A::Union{Ref{F}, Nothing}
n::Union{Ref{F}, Nothing}
A::Union{Ref{Real}, Nothing}
n::Union{Ref{Real}, Nothing}
C::Union{Ref{Real}, Nothing}
H₀::Union{Matrix{F}, Nothing}
H::Union{Matrix{F}, Nothing}
H̄::Union{Matrix{F}, Nothing}
Expand All @@ -31,16 +32,17 @@ mutable struct SIA2Dmodel{F <: AbstractFloat, I <: Integer} <: SIAmodel
V::Union{Matrix{F}, Nothing}
Vx::Union{Matrix{F}, Nothing}
Vy::Union{Matrix{F}, Nothing}
Γ::Union{Ref{F}, Nothing}
Γ::Union{Ref{Real}, Nothing}
MB::Union{Matrix{F}, Nothing}
MB_mask::Union{BitMatrix, Nothing}
MB_total::Union{Matrix{F}, Nothing}
glacier_idx::Union{Ref{I}, Nothing}
end

function SIA2Dmodel(params::Sleipnir.Parameters;
A::Union{F, Nothing} = nothing,
n::Union{F, Nothing} = nothing,
A::Union{Real, Nothing} = nothing,
n::Union{Real, Nothing} = nothing,
C::Union{Real, Nothing} = nothing,
H₀::Union{Matrix{F}, Nothing} = nothing,
H::Union{Matrix{F}, Nothing} = nothing,
H̄::Union{Matrix{F}, Nothing} = nothing,
Expand Down Expand Up @@ -76,14 +78,17 @@ function SIA2Dmodel(params::Sleipnir.Parameters;
if !isnothing(n)
n = Ref{F}(n)
end
if !isnothing(C)
A = Ref{F}(C)
end
if !isnothing(Γ)
Γ = Ref{F}(Γ)
end
if !isnothing(glacier_idx)
glacier_idx = Ref{I}(glacier_idx)
end

SIA2D_model = SIA2Dmodel{ft,it}(A, n, 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
Expand Down Expand Up @@ -112,12 +117,13 @@ function initialize_iceflow_model!(iceflow_model::IF,
) where {IF <: IceflowModel, I <: Int, G <: Sleipnir.AbstractGlacier}
nx, ny = glacier.nx, glacier.ny
F = params.simulation.float_type
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.H₀ = deepcopy(glacier.H₀)::Matrix{F}
iceflow_model.H = deepcopy(glacier.H₀)::Matrix{F}
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.H₀ = deepcopy(glacier.H₀)
iceflow_model.H = deepcopy(glacier.H₀)
iceflow_model.H̄ = zeros(F,nx-1,ny-1)
iceflow_model.S = deepcopy(glacier.S)::Matrix{F}
iceflow_model.S = deepcopy(glacier.S)
iceflow_model.dSdx = zeros(F,nx-1,ny)
iceflow_model.dSdy= zeros(F,nx,ny-1)
iceflow_model.D = zeros(F,nx-1,ny-1)
Expand All @@ -135,7 +141,7 @@ function initialize_iceflow_model!(iceflow_model::IF,
iceflow_model.V = zeros(F,nx-1,ny-1)
iceflow_model.Vx = zeros(F,nx-1,ny-1)
iceflow_model.Vy = zeros(F,nx-1,ny-1)
iceflow_model.Γ = isnothing(iceflow_model.Γ) ? Ref{F}(0.0) : iceflow_model.Γ
iceflow_model.Γ = isnothing(iceflow_model.Γ) ? Ref{Real}(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)
Expand Down Expand Up @@ -165,11 +171,19 @@ function initialize_iceflow_model(iceflow_model::IF,
) where {IF <: IceflowModel, I <: Int}
nx, ny = glacier.nx, glacier.ny
F = params.simulation.float_type
iceflow_model.A = Ref{F}(glacier.A)
iceflow_model.n = Ref{F}(glacier.n)
iceflow_model.glacier_idx = Ref{I}(glacier_idx)
iceflow_model.A = glacier.A
iceflow_model.n = glacier.n
iceflow_model.C = glacier.C
iceflow_model.S = glacier.S
iceflow_model.glacier_idx = glacier_idx
# We just need initial condition to run out-of-place forward model
iceflow_model.H₀ = deepcopy(glacier.H₀)::Matrix{F}
iceflow_model.H = deepcopy(glacier.H₀)::Matrix{F}
iceflow_model.H₀ = deepcopy(glacier.H₀)
iceflow_model.H = deepcopy(glacier.H₀)

iceflow_model.MB = zeros(F,nx,ny)
iceflow_model.MB_mask = zeros(F,nx,ny)
iceflow_model.MB_total = zeros(F,nx,ny)


end

178 changes: 107 additions & 71 deletions src/models/iceflow/SIA2D/SIA2D_utils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,37 +5,37 @@ SIA2D!(dH, H, SIA2Dmodel)
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}
function SIA2D!(dH::Matrix{<:Real}, H::Matrix{<:Real}, simulation::SIM, t::F) where {F <: AbstractFloat, SIM <: Simulation}

# Retrieve parameters
SIA2D_model::SIA2Dmodel = simulation.model.iceflow
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
n::Ref{F} = SIA2D_model.n
B::Matrix{F} = glacier.B
S::Matrix{F} = SIA2D_model.S
dSdx::Matrix{F} = SIA2D_model.dSdx
dSdy::Matrix{F} = SIA2D_model.dSdy
D::Matrix{F} = SIA2D_model.D
Dx::Matrix{F} = SIA2D_model.Dx
Dy::Matrix{F} = SIA2D_model.Dy
dSdx_edges::Matrix{F} = SIA2D_model.dSdx_edges
dSdy_edges::Matrix{F} = SIA2D_model.dSdy_edges
∇S::Matrix{F} = SIA2D_model.∇S
∇Sx::Matrix{F} = SIA2D_model.∇Sx
∇Sy::Matrix{F} = SIA2D_model.∇Sy
Fx::Matrix{F} = SIA2D_model.Fx
Fy::Matrix{F} = SIA2D_model.Fy
Fxx::Matrix{F} = SIA2D_model.Fxx
Fyy::Matrix{F} = SIA2D_model.Fyy
Δx::F = glacier.Δx
Δy::F = glacier.Δy
Γ::Ref{F} = SIA2D_model.Γ
ρ::F = simulation.parameters.physical.ρ
g::F = simulation.parameters.physical.g
H̄ = SIA2D_model.H̄
A = SIA2D_model.A
n = SIA2D_model.n
B = glacier.B
S = SIA2D_model.S
dSdx = SIA2D_model.dSdx
dSdy = SIA2D_model.dSdy
D = SIA2D_model.D
Dx = SIA2D_model.Dx
Dy = SIA2D_model.Dy
dSdx_edges = SIA2D_model.dSdx_edges
dSdy_edges = SIA2D_model.dSdy_edges
∇S = SIA2D_model.∇S
∇Sx = SIA2D_model.∇Sx
∇Sy = SIA2D_model.∇Sy
Fx = SIA2D_model.Fx
Fy = SIA2D_model.Fy
Fxx = SIA2D_model.Fxx
Fyy = SIA2D_model.Fyy
Δx = glacier.Δx
Δy = glacier.Δy
Γ = SIA2D_model.Γ
ρ = simulation.parameters.physical.ρ
g = simulation.parameters.physical.g

# First, enforce values to be positive
map!(x -> ifelse(x>0.0,x,0.0), H, H)
Expand Down Expand Up @@ -88,22 +88,26 @@ end
Compute a step of the Shallow Ice Approximation UDE in a forward model. Allocates memory.
"""
function SIA2D(H::Matrix{F}, simulation::SIM, t::F) where {F <: AbstractFloat, SIM <: Simulation}
function SIA2D(H::Matrix{<:Real}, simulation::SIM, t) where { SIM <: Simulation}
# function SIA2D(H, SIA2Dmodel)
# Retrieve parameters
SIA2D_model::SIA2Dmodel = simulation.model.iceflow
glacier::Sleipnir.Glacier2D = simulation.glaciers[simulation.model.iceflow.glacier_idx[]]
params::Sleipnir.Parameters = simulation.parameters
int_type = simulation.parameters.simulation.int_type
int_type = simulation.parameters.simulation.int_type
# Retrieve parameters
B::Matrix{F} = glacier.B
Δx::F = glacier.Δx
Δy::F = glacier.Δy
A::Ref{F} = SIA2D_model.A
n::Ref{F} = SIA2D_model.n
ρ::F = params.physical.ρ
g::F = params.physical.g
B = glacier.B
Δx = glacier.Δx
Δy = glacier.Δy
A = SIA2D_model.A
n = SIA2D_model.n
C = SIA2D_model.C
ρ = params.physical.ρ
g = params.physical.g

# First, enforce values to be positive
map!(x -> ifelse(x>0.0,x,0.0), H, H)

# Update glacier surface altimetry
S = B .+ H

Expand Down Expand Up @@ -146,11 +150,10 @@ based on the initial and final ice thickness states.
"""
function avg_surface_V!(simulation::SIM) where {SIM <: Simulation}
# TODO: Add more datapoints to better interpolate this
ft = simulation.parameters.simulation.float_type
iceflow_model = simulation.model.iceflow

Vx₀::Matrix{ft}, Vy₀::Matrix{ft} = surface_V!(iceflow_model.H₀, simulation)
Vx::Matrix{ft}, Vy::Matrix{ft} = surface_V!(iceflow_model.H, 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
Expand All @@ -163,14 +166,12 @@ end
Computes the average ice surface velocity for a given glacier evolution period
based on the initial and final ice thickness states.
"""
function avg_surface_V(H::Matrix{F}, simulation::SIM) where {F <: AbstractFloat, SIM <: Simulation}

ft = simulation.parameters.simulation.float_type
function avg_surface_V(H::Matrix{<:Real}, simulation::SIM) where {SIM <: Simulation}
iceflow_model = simulation.model.iceflow

# We compute the initial and final surface velocity and average them
Vx₀::Matrix{ft}, Vy₀::Matrix{ft} = surface_V(iceflow_model.H₀, simulation)
Vx::Matrix{ft}, Vy::Matrix{ft} = surface_V(H, simulation)
Vx₀, Vy₀ = surface_V(iceflow_model.H₀, simulation)
Vx, Vy = surface_V(H, simulation)

V̄x = (Vx₀ .+ Vx)./2.0
V̄y = (Vy₀ .+ Vy)./2.0
Expand All @@ -185,32 +186,31 @@ end
Computes the ice surface velocity for a given glacier state
"""
function surface_V!(H::Matrix{F}, simulation::SIM) where {F <: AbstractFloat, SIM <: Simulation}
function surface_V!(H::Matrix{<:Real}, simulation::SIM) where {SIM <: Simulation}
params::Sleipnir.Parameters = simulation.parameters
ft = params.simulation.float_type
it = params.simulation.int_type
iceflow_model = simulation.model.iceflow
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
dSdy::Matrix{ft} = iceflow_model.dSdy
∇S::Matrix{ft} = iceflow_model.∇S
∇Sx::Matrix{ft} = iceflow_model.∇Sx
∇Sy::Matrix{ft} = iceflow_model.∇Sy
Γꜛ::Ref{ft} = iceflow_model.Γ
D::Matrix{ft} = iceflow_model.D
Dx::Matrix{ft} = iceflow_model.Dx
Dy::Matrix{ft} = iceflow_model.Dy
A::Ref{ft} = iceflow_model.A
n::Ref{ft} = iceflow_model.n
Δx::ft = glacier.Δx
Δy::ft = glacier.Δy
ρ::ft = params.physical.ρ
g::ft = params.physical.g

B = glacier.B
H̄ = iceflow_model.H̄
dSdx = iceflow_model.dSdx
dSdy = iceflow_model.dSdy
∇S = iceflow_model.∇S
∇Sx = iceflow_model.∇Sx
∇Sy = iceflow_model.∇Sy
Γꜛ = iceflow_model.Γ
D = iceflow_model.D
Dx = iceflow_model.Dx
D = iceflow_model.Dy
A = iceflow_model.A
n = iceflow_model.n
Δx = glacier.Δx
Δy = glacier.Δy
ρ = params.physical.ρ
g = params.physical.g

# Update glacier surface altimetry
S::Matrix{ft} = B .+ H
S = B .+ H

# All grid variables computed in a staggered grid
# Compute surface gradients on edges
Expand All @@ -222,7 +222,7 @@ function surface_V!(H::Matrix{F}, simulation::SIM) where {F <: AbstractFloat, SI

avg!(H̄, H)
Γꜛ[] = 2.0 * A[] * (ρ * g)^n[] / (n[]+1) # surface stress (not average) # 1 / m^3 s
D .= Γꜛ[] .* H̄.^(n[] + 1) .* ∇S
D = Γꜛ[] .* H̄.^(n[] + 1) .* ∇S

# Compute averaged surface velocities
Vx = .-D .* ∇Sx
Expand All @@ -236,18 +236,18 @@ end
Computes the ice surface velocity for a given glacier state
"""
function surface_V(H::Matrix{F}, simulation::SIM) where {F <: AbstractFloat, SIM <: Simulation}
function surface_V(H::Matrix{<:Real}, simulation::SIM) where {SIM <: Simulation}
params::Sleipnir.Parameters = simulation.parameters
ft = params.simulation.float_type

iceflow_model = simulation.model.iceflow
glacier::Sleipnir.Glacier2D = simulation.glaciers[iceflow_model.glacier_idx[]]
B = glacier.B
Δx = glacier.Δx
Δy = glacier.Δy
A::Ref{ft} = iceflow_model.A
n::Ref{ft} = iceflow_model.n
ρ::ft = params.physical.ρ
g::ft = params.physical.g
A = iceflow_model.A
n = iceflow_model.n
ρ = params.physical.ρ
g = params.physical.g

# Update glacier surface altimetry
S = B .+ H
Expand All @@ -266,4 +266,40 @@ function surface_V(H::Matrix{F}, simulation::SIM) where {F <: AbstractFloat, SIM
Vy = - D .* avg_x(dSdy)

return Vx, Vy
end
end

function H_from_V(V::Matrix{<:Real}, simulation::SIM) where {SIM <: Simulation}
params::Sleipnir.Parameters = simulation.parameters

iceflow_model = simulation.model.iceflow
glacier::Sleipnir.Glacier2D = simulation.glaciers[iceflow_model.glacier_idx[]]
B = glacier.B
Δx = glacier.Δx
Δy = glacier.Δy
A = iceflow_model.A
n = iceflow_model.n
C = iceflow_model.C
ρ = params.physical.ρ
g = params.physical.g
H₀ = glacier.H₀

# Update glacier surface altimetry
S = iceflow_model.S
V = Huginn.avg(V)

# All grid variables computed in a staggered grid
# Compute surface gradients on edges
dSdx = Huginn.diff_x(S) / Δx
dSdy = Huginn.diff_y(S) / Δy
∇S = (Huginn.avg_y(dSdx).^2 .+ Huginn.avg_x(dSdy).^2).^(1/2)


Γꜛ = (2.0 * A[] * (ρ * g)^n[]) / (n[]+1) # surface stress (not average) # 1 / m^3 s

H = ( V + C[] ./ (Γꜛ .*(∇S .^ n[]))) .^ (1 / (n[] + 1))

replace!(H, NaN=>0.0)
replace!(H, Inf=>0.0)
return H
end

2 changes: 1 addition & 1 deletion src/parameters/SolverParameters.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ function Parameters(;

# Build the parameters based on all the subtypes of parameters
parameters = Sleipnir.Parameters(physical, simulation, OGGM,
nothing, solver, nothing)
nothing, solver, nothing, nothing)

if parameters.simulation.multiprocessing
enable_multiprocessing(parameters)
Expand Down
Loading

0 comments on commit 6ead632

Please sign in to comment.