Skip to content

Commit 0c6ab16

Browse files
avik-palavik-pal
committed
fix: more fixes
1 parent 991884d commit 0c6ab16

File tree

6 files changed

+322
-189
lines changed

6 files changed

+322
-189
lines changed

‎JuliaFormatter.toml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
style="sciml"

‎benchmarks/NonlinearProblem/Manifest.toml‎

Lines changed: 16 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -515,9 +515,9 @@ version = "1.15.1"
515515

516516
[[deps.DifferentiationInterface]]
517517
deps = ["ADTypes", "LinearAlgebra"]
518-
git-tree-sha1 = "e41b6696c84291c4ad15f5f6eaf071b4dfbfda06"
518+
git-tree-sha1 = "aa87a743e3778d35a950b76fbd2ae64f810a2bb3"
519519
uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
520-
version = "0.6.51"
520+
version = "0.6.52"
521521

522522
[deps.DifferentiationInterface.extensions]
523523
DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore"
@@ -527,6 +527,7 @@ version = "0.6.51"
527527
DifferentiationInterfaceFiniteDiffExt = "FiniteDiff"
528528
DifferentiationInterfaceFiniteDifferencesExt = "FiniteDifferences"
529529
DifferentiationInterfaceForwardDiffExt = ["ForwardDiff", "DiffResults"]
530+
DifferentiationInterfaceGPUArraysCoreExt = "GPUArraysCore"
530531
DifferentiationInterfaceGTPSAExt = "GTPSA"
531532
DifferentiationInterfaceMooncakeExt = "Mooncake"
532533
DifferentiationInterfacePolyesterForwardDiffExt = ["PolyesterForwardDiff", "ForwardDiff", "DiffResults"]
@@ -549,6 +550,7 @@ version = "0.6.51"
549550
FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
550551
FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
551552
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
553+
GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
552554
GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8"
553555
Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
554556
PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b"
@@ -578,9 +580,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
578580

579581
[[deps.Distributions]]
580582
deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
581-
git-tree-sha1 = "0b4190661e8a4e51a842070e7dd4fae440ddb7f4"
583+
git-tree-sha1 = "6d8b535fd38293bc54b88455465a1386f8ac1c3c"
582584
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
583-
version = "0.25.118"
585+
version = "0.25.119"
584586

585587
[deps.Distributions.extensions]
586588
DistributionsChainRulesCoreExt = "ChainRulesCore"
@@ -1674,9 +1676,9 @@ version = "0.1.2"
16741676

16751677
[[deps.NonlinearSolve]]
16761678
deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolveQuasiNewton", "NonlinearSolveSpectralMethods", "PrecompileTools", "Preferences", "Reexport", "SciMLBase", "SimpleNonlinearSolve", "SparseArrays", "SparseMatrixColorings", "StaticArraysCore", "SymbolicIndexingInterface"]
1677-
git-tree-sha1 = "7ae7322d658544bd8f6b24a1a0374a6b4ac1fc7e"
1679+
git-tree-sha1 = "4bb539009cc5be9024fb9cbbe7ea2fd1197aa2aa"
16781680
uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
1679-
version = "4.5.1"
1681+
version = "4.6.0"
16801682

16811683
[deps.NonlinearSolve.extensions]
16821684
NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
@@ -1766,9 +1768,9 @@ uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
17661768
version = "0.5.5"
17671769

17681770
[[deps.OffsetArrays]]
1769-
git-tree-sha1 = "a414039192a155fb38c4599a60110f0018c6ec82"
1771+
git-tree-sha1 = "117432e406b5c023f665fa73dc26e79ec3630151"
17701772
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
1771-
version = "1.16.0"
1773+
version = "1.17.0"
17721774
weakdeps = ["Adapt"]
17731775

17741776
[deps.OffsetArrays.extensions]
@@ -2371,9 +2373,9 @@ version = "1.10.0"
23712373

23722374
[[deps.SparseConnectivityTracer]]
23732375
deps = ["ADTypes", "DocStringExtensions", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"]
2374-
git-tree-sha1 = "15dd194e46a5e74b6f7f361e9eb3ed869617ccd5"
2376+
git-tree-sha1 = "cccc976f8fdd51bb3a6c3dcd9e1e7d110582e083"
23752377
uuid = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
2376-
version = "0.6.16"
2378+
version = "0.6.17"
23772379

23782380
[deps.SparseConnectivityTracer.extensions]
23792381
SparseConnectivityTracerDataInterpolationsExt = "DataInterpolations"
@@ -2567,9 +2569,11 @@ version = "7.2.1+1"
25672569

25682570
[[deps.Sundials]]
25692571
deps = ["CEnum", "DataStructures", "DiffEqBase", "Libdl", "LinearAlgebra", "Logging", "PrecompileTools", "Reexport", "SciMLBase", "SparseArrays", "Sundials_jll"]
2570-
git-tree-sha1 = "c135b599cec3558be36eaf86ab1ce7e259ef9534"
2572+
git-tree-sha1 = "87b22845424e6c79fa0dadea81880a17759e701a"
2573+
repo-rev = "ap/kinsol_newton"
2574+
repo-url = "https://github.com/SciML/Sundials.jl.git"
25712575
uuid = "c3572dad-4567-51f8-b174-8c6c989267f4"
2572-
version = "4.27.0"
2576+
version = "4.28.0"
25732577

25742578
[[deps.Sundials_jll]]
25752579
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "SuiteSparse_jll", "libblastrampoline_jll"]

‎benchmarks/NonlinearProblem/Project.toml‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,5 @@ StaticArrays = "1"
6161
Sundials = "4.22"
6262
Symbolics = "5, 6"
6363

64-
[sources.NonlinearSolve]
65-
rev = "ap/petsc_debug"
66-
url = "https://github.com/SciML/NonlinearSolve.jl"
64+
[sources]
65+
Sundials = {url = "https://github.com/SciML/Sundials.jl", rev = "ap/kinsol_newton"}

‎benchmarks/NonlinearProblem/bruss.jmd‎

Lines changed: 2 additions & 173 deletions
Original file line numberDiff line numberDiff line change
@@ -370,7 +370,7 @@ solvers_scaling = [
370370
(; pkg = :nonlinearsolve, sparsity = :none, name = "NR (No Sparsity)", alg = NewtonRaphson()),
371371
(; pkg = :nonlinearsolve, sparsity = :exact, name = "NR (Exact Sparsity)", alg = NewtonRaphson()),
372372
(; pkg = :wrapper, sparsity = :none, name = "NR [NLsolve.jl]", alg = NLsolveJL(; method = :newton, autodiff = :forward)),
373-
(; pkg = :wrapper, sparsity = :none, name = "Mod. NR [Sundials]", alg = KINSOL()),
373+
(; pkg = :wrapper, sparsity = :none, name = "NR [Sundials]", alg = KINSOL(; linear_solver = :LapackDense, maxsetupcalls=1)),
374374
(; pkg = :wrapper, sparsity = :none, name = "NR [PETSc] (No Sparsity)", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", autodiff = missing)),
375375
(; pkg = :wrapper, sparsity = :exact, name = "NR [PETSc] (Exact Sparsity)", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic")),
376376

@@ -534,177 +534,6 @@ end
534534
save("brusselator_scaling.svg", fig)
535535
```
536536

537-
# Jacobian-Free Newton / TR Krylov Methods
538-
539-
In this section, we will benchmark jacobian-free nonlinear solvers with Krylov methods. We
540-
will use preconditioning from `AlgebraicMultigrid.jl` and `IncompleteLU.jl`. Unfortunately,
541-
our ability to use 3rd party software is limited here, since only `Sundials.jl` supports
542-
jacobian-free methods via `:GMRES`.
543-
544-
```julia
545-
using AlgebraicMultigrid, IncompleteLU
546-
547-
incompletelu(W, p = nothing) = ilu(W, τ = 50.0), LinearAlgebra.I
548-
549-
function algebraicmultigrid(W, p = nothing)
550-
return aspreconditioner(ruge_stuben(convert(AbstractMatrix, W))), LinearAlgebra.I
551-
end
552-
553-
function algebraicmultigrid_jacobi(W, p = nothing)
554-
A = convert(AbstractMatrix, W)
555-
Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(
556-
A, presmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1))),
557-
postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1)))
558-
))
559-
return Pl, LinearAlgebra.I
560-
end
561-
562-
Ns = 2 .^ (2:7)
563-
564-
solvers_scaling_jacobian_free = [
565-
(; pkg = :nonlinearsolve, name = "Newton Krylov", alg = NewtonRaphson(; linsolve = KrylovJL_GMRES())),
566-
(; pkg = :nonlinearsolve, name = "Newton Krylov (ILU)", alg = NewtonRaphson(; linsolve = KrylovJL_GMRES(; precs = incompletelu), concrete_jac = true)),
567-
(; pkg = :nonlinearsolve, name = "Newton Krylov (AMG)", alg = NewtonRaphson(; linsolve = KrylovJL_GMRES(; precs = algebraicmultigrid), concrete_jac = true)),
568-
# (; pkg = :nonlinearsolve, name = "Newton Krylov (AMG Jacobi)", alg = NewtonRaphson(; linsolve = KrylovJL_GMRES(; precs = algebraicmultigrid_jacobi), concrete_jac = true)),
569-
(; pkg = :wrapper, name = "Newton Krylov [PETSc]", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", ksp_type = "gmres", snes_mf = true)),
570-
(; pkg = :wrapper, name = "Newton Krylov (ILU) [PETSc]", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", ksp_type = "gmres", pc_type = "ilu")),
571-
(; pkg = :wrapper, name = "Newton Krylov (AMG) [PETSc]", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", ksp_type = "gmres", pc_type = "gamg")),
572-
# (; pkg = :wrapper, name = "Newton Krylov (AMG Jacobi) [PETSc]", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", ksp_type = "gmres", pc_type = "gamg", mg_levels_ksp_type = "richardson", mg_levels_pc_type = "jacobi")),
573-
(; pkg = :nonlinearsolve, name = "TR Krylov", alg = TrustRegion(; linsolve = KrylovJL_GMRES())),
574-
(; pkg = :nonlinearsolve, name = "TR Krylov (ILU)", alg = TrustRegion(; linsolve = KrylovJL_GMRES(; precs = incompletelu), concrete_jac = true)),
575-
(; pkg = :nonlinearsolve, name = "TR Krylov (AMG)", alg = TrustRegion(; linsolve = KrylovJL_GMRES(; precs = algebraicmultigrid), concrete_jac = true)),
576-
# (; pkg = :nonlinearsolve, name = "TR Krylov (AMG Jacobi)", alg = TrustRegion(; linsolve = KrylovJL_GMRES(; precs = algebraicmultigrid_jacobi), concrete_jac = true)),
577-
# (; pkg = :wrapper, name = "TR Krylov [PETSc]", alg = PETScSNES(; snes_type = "newtontr", ksp_type = "gmres", snes_mf = true)),
578-
# (; pkg = :wrapper, name = "TR Krylov (ILU) [PETSc]", alg = PETScSNES(; snes_type = "newtontr", ksp_type = "gmres", pc_type = "ilu")),
579-
# (; pkg = :wrapper, name = "TR Krylov (AMG) [PETSc]", alg = PETScSNES(; snes_type = "newtontr", ksp_type = "gmres", pc_type = "gamg")),
580-
# (; pkg = :wrapper, name = "TR Krylov (AMG Jacobi) [PETSc]", alg = PETScSNES(; snes_type = "newtontr", ksp_type = "gmres", pc_type = "gamg", mg_levels_ksp_type = "richardson", mg_levels_pc_type = "jacobi")),
581-
(; pkg = :wrapper, name = "Newton Krylov [Sundials]", alg = KINSOL(; linear_solver = :GMRES)),
582-
]
583-
584-
runtimes_scaling = zeros(length(solvers_scaling_jacobian_free), length(Ns)) .- 1
585-
586-
for (i, N) in enumerate(Ns)
587-
prob = generate_brusselator_problem(
588-
N; sparsity = TracerSparsityDetector()
589-
)
590-
591-
@info "Benchmarking N = $N"
592-
593-
for (j, solver) in enumerate(solvers_scaling_jacobian_free)
594-
alg = solver.alg
595-
name = solver.name
596-
597-
if (j > 1 && runtimes_scaling[j - 1, i] == -1)
598-
# The last benchmark failed so skip this too
599-
runtimes_scaling[j, i] = NaN
600-
@warn "$(name): Would Have Timed out"
601-
else
602-
function benchmark_function()
603-
termination_condition = (alg isa PETScSNES || alg isa KINSOL) ?
604-
nothing :
605-
NonlinearSolveBase.AbsNormTerminationMode(Base.Fix1(maximum, abs))
606-
sol = solve(prob, alg; abstol=1e-5, reltol=1e-5,
607-
linsolve_kwargs = (; abstol = 1e-6, reltol = 1e-6),
608-
termination_condition)
609-
if SciMLBase.successful_retcode(sol) || norm(sol.resid, Inf) ≤ 1e-4
610-
runtimes_scaling[j, i] = @belapsed solve($prob, $alg; abstol=1e-5,
611-
reltol=1e-5,
612-
linsolve_kwargs = (; abstol = 1e-6, reltol = 1e-6),
613-
termination_condition=$termination_condition)
614-
else
615-
runtimes_scaling[j, i] = NaN
616-
end
617-
@info "$(name): $(runtimes_scaling[j, i]) | $(norm(sol.resid, Inf)) | $(sol.retcode)"
618-
end
619-
620-
timeout(benchmark_function, 600)
621-
622-
if runtimes_scaling[j, i] == -1
623-
@warn "$(name): Timed out"
624-
runtimes_scaling[j, i] = NaN
625-
end
626-
end
627-
end
628-
629-
println()
630-
end
631-
```
632-
633-
Plot the results.
634-
635-
```julia
636-
fig = begin
637-
ASPECT_RATIO = 0.7
638-
WIDTH = 1200
639-
HEIGHT = round(Int, WIDTH * ASPECT_RATIO)
640-
STROKEWIDTH = 2.5
641-
642-
cycle = Cycle([:marker], covary = true)
643-
colors = cgrad(:tableau_20, length(solvers_scaling_jacobian_free); categorical = true)
644-
theme = Theme(Lines = (cycle = cycle,), Scatter = (cycle = cycle,))
645-
LINESTYLES = Dict(
646-
(:nonlinearsolve, :none) => :solid,
647-
(:nonlinearsolve, :amg) => :dot,
648-
(:nonlinearsolve, :amg_jacobi) => :dash,
649-
(:nonlinearsolve, :ilu) => :dashdot,
650-
)
651-
652-
Ns_ = Ns .^ 2 .* 2
653-
654-
with_theme(theme) do
655-
fig = Figure(; size = (WIDTH, HEIGHT))
656-
657-
ax = Axis(fig[1, 1:2], ylabel = L"Time ($s$)", xlabel = L"Problem Size ($N$)",
658-
xscale = log10, yscale = log10, xlabelsize = 22, ylabelsize = 22,
659-
xticklabelsize = 20, yticklabelsize = 20, xtickwidth = STROKEWIDTH,
660-
ytickwidth = STROKEWIDTH, spinewidth = STROKEWIDTH)
661-
662-
idxs = get_ordering(runtimes_scaling)
663-
664-
ls, scs, labels = [], [], []
665-
for (i, solver) in zip(idxs, solvers_scaling_jacobian_free[idxs])
666-
all(isnan, runtimes_scaling[i, :]) && continue
667-
precon = occursin("AMG Jacobi", solver.name) ? :amg_jacobi : occursin("AMG", solver.name) ? :amg : occursin("ILU", solver.name) ? :ilu : :none
668-
linestyle = LINESTYLES[(solver.pkg, precon)]
669-
l = lines!(Ns_, runtimes_scaling[i, :]; linewidth = 5, color = colors[i],
670-
linestyle)
671-
sc = scatter!(Ns_, runtimes_scaling[i, :]; markersize = 16, strokewidth = 2,
672-
color = colors[i])
673-
push!(ls, l)
674-
push!(scs, sc)
675-
push!(labels, solver.name)
676-
end
677-
678-
axislegend(ax, [[l, sc] for (l, sc) in zip(ls, scs)], labels,
679-
"Successful Solvers";
680-
framevisible=true, framewidth = STROKEWIDTH, orientation = :vertical,
681-
titlesize = 20, labelsize = 16, position = :rb,
682-
tellheight = true, tellwidth = false, patchsize = (40.0f0, 20.0f0))
683-
684-
axislegend(ax, [
685-
LineElement(; linestyle = :solid, linewidth = 5),
686-
LineElement(; linestyle = :dot, linewidth = 5),
687-
LineElement(; linestyle = :dash, linewidth = 5),
688-
LineElement(; linestyle = :dashdot, linewidth = 5),
689-
], ["No Preconditioning", "AMG", "AMG Jacobi", "Incomplete LU"],
690-
"Preconditioning"; framevisible=true, framewidth = STROKEWIDTH,
691-
orientation = :vertical, titlesize = 20, labelsize = 16,
692-
tellheight = true, tellwidth = true, patchsize = (40.0f0, 20.0f0),
693-
position = :lt)
694-
695-
fig[0, :] = Label(fig,
696-
"Brusselator 2D: Scaling of Jacobian-Free Nonlinear Solvers with Problem Size",
697-
fontsize = 24, tellwidth = false, font = :bold)
698-
699-
return fig
700-
end
701-
end
702-
```
703-
704-
```julia
705-
save("brusselator_krylov_methods_scaling.svg", fig)
706-
```
707-
708537
# Work-Precision Diagram
709538

710539
In this section, we will generate the work-precision of the solvers. All solvers that can
@@ -723,7 +552,7 @@ solvers_all = [
723552
(; pkg = :wrapper, name = "NR [NLsolve.jl]", solver = Dict(:alg => NLsolveJL(; method = :newton, autodiff = :forward))),
724553
(; pkg = :wrapper, name = "TR [NLsolve.jl]", solver = Dict(:alg => NLsolveJL(; autodiff = :forward))),
725554

726-
(; pkg = :wrapper, name = "NR [Sundials]", solver = Dict(:alg => KINSOL())),
555+
(; pkg = :wrapper, name = "NR [Sundials]", solver = Dict(:alg => KINSOL(; linear_solver = :LapackDense, maxsetupcalls=1))),
727556
(; pkg = :wrapper, name = "Newton Krylov [Sundials]", solver = Dict(:alg => KINSOL(; linear_solver = :GMRES))),
728557

729558
(; pkg = :wrapper, name = "Mod. Powell [MINPACK]", solver = Dict(:alg => CMINPACK())),

0 commit comments

Comments
 (0)