From 91ae280463f4e15110166a346c74170df25a1001 Mon Sep 17 00:00:00 2001
From: Chris Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 4 Dec 2022 02:28:31 +0100
Subject: [PATCH] push new test

---
 Project.toml     |  4 +++-
 test/runtests.jl |  1 +
 test/sparse.jl   | 43 +++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 47 insertions(+), 1 deletion(-)
 create mode 100644 test/sparse.jl

diff --git a/Project.toml b/Project.toml
index d63cf5913..11fb4ef68 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,6 +15,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
@@ -40,7 +41,8 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays"]
+test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics"]
diff --git a/test/runtests.jl b/test/runtests.jl
index 7269ed6bf..37b9957e7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,4 +9,5 @@ const is_APPVEYOR = Sys.iswindows() && haskey(ENV, "APPVEYOR")
 
 if GROUP == "All" || GROUP == "Core"
     @time @safetestset "Basic Tests + Some AD" begin include("basictests.jl") end
+    @time @safetestset "Sparsity Tests" begin include("sparse.jl") end
 end end
diff --git a/test/sparse.jl b/test/sparse.jl
new file mode 100644
index 000000000..24c16408a
--- /dev/null
+++ b/test/sparse.jl
@@ -0,0 +1,43 @@
+using NonlinearSolve, LinearAlgebra, SparseArrays
+
+const N = 32
+const xyd_brusselator = range(0,stop=1,length=N)
+brusselator_f(x, y) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * 5.
+limit(a, N) = a == N+1 ? 1 : a == 0 ? N : a
+function brusselator_2d_loop(du, u, p)
+  A, B, alpha, dx = p
+  alpha = alpha/dx^2
+  @inbounds for I in CartesianIndices((N, N))
+    i, j = Tuple(I)
+    x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
+    ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N)
+    du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
+                B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y)
+    du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
+                A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
+    end
+end
+p = (3.4, 1., 10., step(xyd_brusselator))
+
+function init_brusselator_2d(xyd)
+  N = length(xyd)
+  u = zeros(N, N, 2)
+  for I in CartesianIndices((N, N))
+    x = xyd[I[1]]
+    y = xyd[I[2]]
+    u[I,1] = 22*(y*(1-y))^(3/2)
+    u[I,2] = 27*(x*(1-x))^(3/2)
+  end
+  u
+end
+u0 = init_brusselator_2d(xyd_brusselator)
+prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop,u0,p)
+sol = solve(prob_brusselator_2d, NewtonRaphson())
+
+using Symbolics
+du0 = copy(u0)
+jac_sparsity = Symbolics.jacobian_sparsity((du,u)->brusselator_2d_loop(du,u,p),du0,u0)
+
+f = NonlinearFunction(brusselator_2d_loop;jac_prototype=float.(jac_sparsity))
+prob_brusselator_2d = NonlinearProblem(f,u0,p)
+sol = solve(prob_brusselator_2d, NewtonRaphson())