From b1bd458e9e2f30207bc5bdc01c400223aa91a95a Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Tue, 28 May 2024 17:50:29 +0200
Subject: [PATCH 1/3] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 9dfd094..39e8cff 100644
--- a/Project.toml
+++ b/Project.toml
@@ -31,7 +31,7 @@ Infiltrator = "1"
 JLD2 = "0.4"
 Muninn = "0.3"
 Optimization = "3"
-OptimizationOptimJL = "0.3"
+OptimizationOptimJL = "0.1.5"
 OrdinaryDiffEq = "6"
 PlotThemes = "3"
 Plots = "1"

From 2b005eb526120ed59ca93086b044b93aa6805aee Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sat, 1 Jun 2024 16:21:45 +0200
Subject: [PATCH 2/3] Use a function barrier and workaround Tullio

---
 src/models/iceflow/SIA2D/SIA2D_utils.jl | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/src/models/iceflow/SIA2D/SIA2D_utils.jl b/src/models/iceflow/SIA2D/SIA2D_utils.jl
index 3a4841a..80a463e 100644
--- a/src/models/iceflow/SIA2D/SIA2D_utils.jl
+++ b/src/models/iceflow/SIA2D/SIA2D_utils.jl
@@ -106,7 +106,10 @@ function SIA2D(H::Matrix{R}, simulation::SIM, t::R; batch_id::Union{Nothing, I}
     n = SIA2D_model.n
     ρ = params.physical.ρ
     g = params.physical.g
+    _SIA2D(glacier, SIA2D_model, H, B, Δx, Δy, A, n, ρ, g)
+end
 
+function _SIA2D(glacier, SIA2D_model, H, B, Δx, Δy, A, n, ρ, g)
     @views H = ifelse.(H.<0.0, 0.0, H) # prevent values from going negative
 
     # First, enforce values to be positive
@@ -141,9 +144,13 @@ function SIA2D(H::Matrix{R}, simulation::SIM, t::R; batch_id::Union{Nothing, I}
     Fy = .-avg_x(D) .* dSdy_edges 
 
     #  Flux divergence
-    @tullio dH[i,j] := -(diff_x(Fx)[pad(i-1,1,1),pad(j-1,1,1)] / Δx + diff_y(Fy)[pad(i-1,1,1),pad(j-1,1,1)] / Δy) 
-
-    return dH
+    pad(i,lower,upper) = ifelse(i<lower,lower,ifelse(i>upper, upper, i))
+    dH = similar(H)
+    for i in 1:size(H,1), j in 1:size(H,2)
+        dH[i,j] = -(diff_x(Fx)[pad(i-1,1,1),pad(j-1,1,1)] / Δx + diff_y(Fy)[pad(i-1,1,1),pad(j-1,1,1)] / Δy)
+    end
+    dH
+    #@tullio threads=false avx=false tensor=false dH[i,j] := -(diff_x(Fx)[pad(i-1,1,1),pad(j-1,1,1)] / Δx + diff_y(Fy)[pad(i-1,1,1),pad(j-1,1,1)] / Δy) 
 end
 
 """

From 9df258a4bfa6351f3e986d92d7347c5c9ad72111 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 16 Jun 2024 09:58:48 -0400
Subject: [PATCH 3/3] Revert "Use a function barrier and workaround Tullio"

This reverts commit 2b005eb526120ed59ca93086b044b93aa6805aee.
---
 src/models/iceflow/SIA2D/SIA2D_utils.jl | 13 +++----------
 1 file changed, 3 insertions(+), 10 deletions(-)

diff --git a/src/models/iceflow/SIA2D/SIA2D_utils.jl b/src/models/iceflow/SIA2D/SIA2D_utils.jl
index 80a463e..3a4841a 100644
--- a/src/models/iceflow/SIA2D/SIA2D_utils.jl
+++ b/src/models/iceflow/SIA2D/SIA2D_utils.jl
@@ -106,10 +106,7 @@ function SIA2D(H::Matrix{R}, simulation::SIM, t::R; batch_id::Union{Nothing, I}
     n = SIA2D_model.n
     ρ = params.physical.ρ
     g = params.physical.g
-    _SIA2D(glacier, SIA2D_model, H, B, Δx, Δy, A, n, ρ, g)
-end
 
-function _SIA2D(glacier, SIA2D_model, H, B, Δx, Δy, A, n, ρ, g)
     @views H = ifelse.(H.<0.0, 0.0, H) # prevent values from going negative
 
     # First, enforce values to be positive
@@ -144,13 +141,9 @@ function _SIA2D(glacier, SIA2D_model, H, B, Δx, Δy, A, n, ρ, g)
     Fy = .-avg_x(D) .* dSdy_edges 
 
     #  Flux divergence
-    pad(i,lower,upper) = ifelse(i<lower,lower,ifelse(i>upper, upper, i))
-    dH = similar(H)
-    for i in 1:size(H,1), j in 1:size(H,2)
-        dH[i,j] = -(diff_x(Fx)[pad(i-1,1,1),pad(j-1,1,1)] / Δx + diff_y(Fy)[pad(i-1,1,1),pad(j-1,1,1)] / Δy)
-    end
-    dH
-    #@tullio threads=false avx=false tensor=false dH[i,j] := -(diff_x(Fx)[pad(i-1,1,1),pad(j-1,1,1)] / Δx + diff_y(Fy)[pad(i-1,1,1),pad(j-1,1,1)] / Δy) 
+    @tullio dH[i,j] := -(diff_x(Fx)[pad(i-1,1,1),pad(j-1,1,1)] / Δx + diff_y(Fy)[pad(i-1,1,1),pad(j-1,1,1)] / Δy) 
+
+    return dH
 end
 
 """