Skip to content

Commit 3900d24

Browse files
avik-palavik-pal
committed
fix: more options
1 parent db5e40f commit 3900d24

File tree

5 files changed

+66
-26
lines changed

5 files changed

+66
-26
lines changed

‎.gitignore‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,4 +22,5 @@ pdf/
2222
notebook/
2323
markdown/
2424

25-
.vscode
25+
.vscode
26+
*.html

‎benchmarks/NonlinearProblem/Manifest.toml‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2223,9 +2223,9 @@ version = "0.2.0"
22232223

22242224
[[deps.SciMLBase]]
22252225
deps = ["ADTypes", "Accessors", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "Moshi", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"]
2226-
git-tree-sha1 = "c4ce89e19f2a220e34c0f377fb41516b642ec899"
2226+
git-tree-sha1 = "f2db9ab9d33130df3be35be9438da51a3416d528"
22272227
uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
2228-
version = "2.83.1"
2228+
version = "2.84.0"
22292229

22302230
[deps.SciMLBase.extensions]
22312231
SciMLBaseChainRulesCoreExt = "ChainRulesCore"

‎benchmarks/NonlinearProblem/Project.toml‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,3 +59,6 @@ StableRNGs = "1"
5959
StaticArrays = "1"
6060
Sundials = "4.22"
6161
Symbolics = "5, 6"
62+
63+
[sources]
64+
NonlinearSolve = {url = "https://github.com/SciML/NonlinearSolve.jl", rev = "ap/petsc_debug"}

‎benchmarks/NonlinearProblem/bruss.jmd‎

Lines changed: 29 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -360,13 +360,13 @@ solvers_scaling = [
360360
(; pkg = :nonlinearsolve, sparsity = :exact, name = "NR (Exact Sparsity)", alg = NewtonRaphson()),
361361
(; pkg = :wrapper, sparsity = :none, name = "NR [NLsolve.jl]", alg = NLsolveJL(; method = :newton, autodiff = :forward)),
362362
(; pkg = :wrapper, sparsity = :none, name = "Mod. NR [Sundials]", alg = KINSOL()),
363-
(; pkg = :wrapper, sparsity = :none, name = "NR [PETSc] (No Sparsity)", alg = PETScSNES(; snes_type = "newtonls")),
364-
(; pkg = :wrapper, sparsity = :exact, name = "NR [PETSc] (Exact Sparsity)", alg = PETScSNES(; snes_type = "newtonls")),
363+
(; pkg = :wrapper, sparsity = :none, name = "NR [PETSc] (No Sparsity)", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", autodiff = missing)),
364+
(; pkg = :wrapper, sparsity = :exact, name = "NR [PETSc] (Exact Sparsity)", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic")),
365365

366366
(; pkg = :nonlinearsolve, sparsity = :none, name = "TR (No Sparsity)", alg = TrustRegion(; radius_update_scheme = RUS.NLsolve)),
367367
(; pkg = :nonlinearsolve, sparsity = :exact, name = "TR (Exact Sparsity)", alg = TrustRegion(; radius_update_scheme = RUS.NLsolve)),
368368
(; pkg = :wrapper, sparsity = :none, name = "TR [NLsolve.jl]", alg = NLsolveJL(; autodiff = :forward)),
369-
(; pkg = :wrapper, sparsity = :none, name = "TR [PETSc] (No Sparsity)", alg = PETScSNES(; snes_type = "newtontr")),
369+
(; pkg = :wrapper, sparsity = :none, name = "TR [PETSc] (No Sparsity)", alg = PETScSNES(; snes_type = "newtontr", autodiff = missing)),
370370
(; pkg = :wrapper, sparsity = :exact, name = "TR [PETSc] (Exact Sparsity)", alg = PETScSNES(; snes_type = "newtontr")),
371371

372372
(; pkg = :wrapper, sparsity = :none, name = "Mod. Powell [MINPACK]", alg = CMINPACK()),
@@ -405,7 +405,7 @@ for (i, N) in enumerate(Ns)
405405
(alg isa GeneralizedFirstOrderAlgorithm && alg.name == :TrustRegion && N > 64) ||
406406
(alg isa NLsolveJL && N > 128 && alg.method == :newton) ||
407407
(alg isa GeneralizedFirstOrderAlgorithm && alg.name == :NewtonRaphson && N > 128 && ptype == :none) ||
408-
(alg isa PETScSNES && N > 64 && ptype == :exact)
408+
(alg isa PETScSNES && N > 64)
409409
# The last benchmark failed so skip this too
410410
runtimes_scaling[j, i] = NaN
411411
@warn "$(name): Would Have Timed out"
@@ -551,15 +551,23 @@ end
551551
Ns = 2 .^ (2:7)
552552

553553
solvers_scaling_jacobian_free = [
554-
(; pkg = :nonlinearsolve, name = "Newton Krylov", alg = NewtonRaphson(; linsolve = KrylovJL_GMRES())),
555-
(; pkg = :nonlinearsolve, name = "Newton Krylov (ILU)", alg = NewtonRaphson(; linsolve = KrylovJL_GMRES(; precs = incompletelu), concrete_jac = true)),
556-
(; pkg = :nonlinearsolve, name = "Newton Krylov (AMG)", alg = NewtonRaphson(; linsolve = KrylovJL_GMRES(; precs = algebraicmultigrid), concrete_jac = true)),
557-
(; pkg = :nonlinearsolve, name = "Newton Krylov (AMG Jacobi)", alg = NewtonRaphson(; linsolve = KrylovJL_GMRES(; precs = algebraicmultigrid_jacobi), concrete_jac = true)),
558-
(; pkg = :nonlinearsolve, name = "TR Krylov", alg = TrustRegion(; linsolve = KrylovJL_GMRES())),
559-
(; pkg = :nonlinearsolve, name = "TR Krylov (ILU)", alg = TrustRegion(; linsolve = KrylovJL_GMRES(; precs = incompletelu), concrete_jac = true)),
560-
(; pkg = :nonlinearsolve, name = "TR Krylov (AMG)", alg = TrustRegion(; linsolve = KrylovJL_GMRES(; precs = algebraicmultigrid), concrete_jac = true)),
561-
(; pkg = :nonlinearsolve, name = "TR Krylov (AMG Jacobi)", alg = TrustRegion(; linsolve = KrylovJL_GMRES(; precs = algebraicmultigrid_jacobi), concrete_jac = true)),
562-
(; pkg = :wrapper, name = "Newton Krylov [Sundials]", alg = KINSOL(; linear_solver = :GMRES)),
554+
(; pkg = :nonlinearsolve, name = "Newton Krylov", alg = NewtonRaphson(; linsolve = KrylovJL_GMRES())),
555+
(; pkg = :nonlinearsolve, name = "Newton Krylov (ILU)", alg = NewtonRaphson(; linsolve = KrylovJL_GMRES(; precs = incompletelu), concrete_jac = true)),
556+
(; pkg = :nonlinearsolve, name = "Newton Krylov (AMG)", alg = NewtonRaphson(; linsolve = KrylovJL_CG(; precs = algebraicmultigrid), concrete_jac = true)),
557+
(; pkg = :nonlinearsolve, name = "Newton Krylov (AMG Jacobi)", alg = NewtonRaphson(; linsolve = KrylovJL_CG(; precs = algebraicmultigrid_jacobi), concrete_jac = true)),
558+
(; pkg = :wrapper, name = "Newton Krylov [PETSc]", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", ksp_type = "gmres", snes_mf = true)),
559+
(; pkg = :wrapper, name = "Newton Krylov (ILU) [PETSc]", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", ksp_type = "gmres", pc_type = "ilu")),
560+
(; pkg = :wrapper, name = "Newton Krylov (AMG) [PETSc]", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", ksp_type = "cg", pc_type = "gamg")),
561+
(; pkg = :wrapper, name = "Newton Krylov (AMG Jacobi) [PETSc]", alg = PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", ksp_type = "cg", pc_type = "gamg", mg_levels_ksp_type = "richardson", mg_levels_pc_type = "jacobi")),
562+
(; pkg = :nonlinearsolve, name = "TR Krylov", alg = TrustRegion(; linsolve = KrylovJL_GMRES())),
563+
(; pkg = :nonlinearsolve, name = "TR Krylov (ILU)", alg = TrustRegion(; linsolve = KrylovJL_GMRES(; precs = incompletelu), concrete_jac = true)),
564+
(; pkg = :nonlinearsolve, name = "TR Krylov (AMG)", alg = TrustRegion(; linsolve = KrylovJL_CG(; precs = algebraicmultigrid), concrete_jac = true)),
565+
(; pkg = :nonlinearsolve, name = "TR Krylov (AMG Jacobi)", alg = TrustRegion(; linsolve = KrylovJL_CG(; precs = algebraicmultigrid_jacobi), concrete_jac = true)),
566+
(; pkg = :wrapper, name = "TR Krylov [PETSc]", alg = PETScSNES(; snes_type = "newtontr", ksp_type = "gmres", snes_mf = true)),
567+
(; pkg = :wrapper, name = "TR Krylov (ILU) [PETSc]", alg = PETScSNES(; snes_type = "newtontr", ksp_type = "gmres", pc_type = "ilu")),
568+
(; pkg = :wrapper, name = "TR Krylov (AMG) [PETSc]", alg = PETScSNES(; snes_type = "newtontr", ksp_type = "cg", pc_type = "gamg")),
569+
(; pkg = :wrapper, name = "TR Krylov (AMG Jacobi) [PETSc]", alg = PETScSNES(; snes_type = "newtontr", ksp_type = "cg", pc_type = "gamg", mg_levels_ksp_type = "richardson", mg_levels_pc_type = "jacobi")),
570+
(; pkg = :wrapper, name = "Newton Krylov [Sundials]", alg = KINSOL(; linear_solver = :GMRES)),
563571
]
564572

565573
runtimes_scaling = zeros(length(solvers_scaling_jacobian_free), length(Ns)) .- 1
@@ -695,16 +703,24 @@ exploit sparsity will automatically do so.
695703
solvers_all = [
696704
(; pkg = :nonlinearsolve, name = "Default PolyAlg", solver = Dict(:alg => FastShortcutNonlinearPolyalg())),
697705
(; pkg = :nonlinearsolve, name = "RobustMultiNewton (GMRES)", solver = Dict(:alg => RobustMultiNewton(; linsolve = KrylovJL_GMRES()))),
706+
698707
(; pkg = :nonlinearsolve, name = "Newton Raphson", solver = Dict(:alg => NewtonRaphson(; linsolve = nothing))),
699708
(; pkg = :nonlinearsolve, name = "Newton Krylov", solver = Dict(:alg => NewtonRaphson(; linsolve = KrylovJL_GMRES()))),
700709
(; pkg = :nonlinearsolve, name = "Trust Region", solver = Dict(:alg => TrustRegion())),
701710
(; pkg = :nonlinearsolve, name = "TR Krylov", solver = Dict(:alg => TrustRegion(; linsolve = KrylovJL_GMRES()))),
711+
702712
(; pkg = :wrapper, name = "NR [NLsolve.jl]", solver = Dict(:alg => NLsolveJL(; method = :newton, autodiff = :forward))),
703713
(; pkg = :wrapper, name = "TR [NLsolve.jl]", solver = Dict(:alg => NLsolveJL(; autodiff = :forward))),
714+
704715
(; pkg = :wrapper, name = "NR [Sundials]", solver = Dict(:alg => KINSOL())),
705716
(; pkg = :wrapper, name = "Newton Krylov [Sundials]", solver = Dict(:alg => KINSOL(; linear_solver = :GMRES))),
706717

707718
(; pkg = :wrapper, name = "Mod. Powell [MINPACK]", solver = Dict(:alg => CMINPACK())),
719+
720+
(; pkg = :wrapper, name = "NR [PETSc]", solver = Dict(:alg => PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", autodiff = missing))),
721+
(; pkg = :wrapper, name = "TR [PETSc]", solver = Dict(:alg => PETScSNES(; snes_type = "newtontr", autodiff = missing))),
722+
(; pkg = :wrapper, name = "Newton Krylov [PETSc]", solver = Dict(:alg => PETScSNES(; snes_type = "newtonls", snes_linesearch_type = "basic", ksp_type = "gmres", autodiff = missing, snes_mf = true))),
723+
(; pkg = :wrapper, name = "TR Krylov [PETSc]", solver = Dict(:alg => PETScSNES(; snes_type = "newtontr", ksp_type = "gmres", autodiff = missing, snes_mf = true))),
708724
];
709725
```
710726

‎benchmarks/NonlinearProblem/nonlinear_solver_23_tests.jmd‎

Lines changed: 30 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -21,12 +21,12 @@ Furthermore, for NonlinearSolve.jl's Newton Raphson method we try the following
2121
- `BackTracking`
2222

2323
and for NonlinearSolve.jl's Trust Region we try the following Radius Update schemes (in addition to the default):
24-
- `NLsolve`
25-
- `NocedalWright`
26-
- `Hei`
27-
- `Yuan`
28-
- `Bastin`
29-
- `Fan`
24+
- `NLsolve`
25+
- `NocedalWright`
26+
- `Hei`
27+
- `Yuan`
28+
- `Bastin`
29+
- `Fan`
3030
and finally for NonlinearSolve.jl's Levenberg-Marquardt method why try using both the default `α_geodesic` value (`0.75`) and a modified value (`0.5`), and also with and without setting the `CholeskyFactorization` linear solver.
3131

3232
For each benchmarked problem, the second, third, and fourth plots compares the performance of NonlinearSolve's Newton Raphson, Trust Region, and Levenberg-Marquardt methods, respectively. The first plot compares the best methods from each of these categories to the various methods available from other packages. At the end of the benchmarks, we print a summary table of which solvers succeeded for which problems.
@@ -38,7 +38,7 @@ Fetch required packages.
3838
```julia
3939
using NonlinearSolve, LinearSolve, StaticArrays, Sundials, Setfield,
4040
BenchmarkTools, LinearAlgebra, DiffEqDevTools, NonlinearProblemLibrary, CairoMakie,
41-
RecursiveFactorization
41+
RecursiveFactorization, Enzyme
4242
import PolyesterForwardDiff, MINPACK, NLsolve, LineSearches
4343

4444
const RUS = RadiusUpdateSchemes;
@@ -92,13 +92,34 @@ reltols = 1.0 ./ 10.0 .^ (4:12);
9292
Prepares various helper functions for benchmarking a specific problem.
9393

9494
```julia
95+
function set_ad_chunksize(solvers, u0)
96+
ck = NonlinearSolve.pickchunksize(u0)
97+
for i in eachindex(solvers)
98+
@set! solvers[i].solver[:alg] = __set_ad_chunksize(solvers[i].solver[:alg], ck)
99+
end
100+
return solvers
101+
end
102+
103+
function __set_ad_chunksize(solver::GeneralizedFirstOrderAlgorithm, ck)
104+
ad = AutoPolyesterForwardDiff(; chunksize = ck)
105+
return GeneralizedFirstOrderAlgorithm(; solver.descent, solver.linesearch,
106+
solver.trustregion, jvp_autodiff = ad, solver.max_shrink_times, solver.vjp_autodiff,
107+
concrete_jac = solver.concrete_jac, name = solver.name)
108+
end
109+
function __set_ad_chunksize(solver::NonlinearSolvePolyAlgorithm, ck)
110+
algs = [__set_ad_chunksize(alg, ck) for alg in solver.algs]
111+
return NonlinearSolvePolyAlgorithm(algs; solver.start_index)
112+
end
113+
__set_ad_chunksize(solver, ck) = solver
114+
95115
# Benchmarks a specific problem, checks which solvers can solve it and their performance
96116
function benchmark_problem!(prob_name; solver_tracker=solver_tracker)
97117
# Finds the problem and the true solution.
98118
prob = nlprob_23_testcases[prob_name]
99119

100120
# Finds the solvers that can solve the problem
101-
successful_solvers = filter(Base.Fix1(check_solver, prob), solvers_all);
121+
solvers_concrete = set_ad_chunksize(solvers_all, prob.prob.u0);
122+
successful_solvers = filter(solver -> check_solver(prob, solver), solvers_concrete);
102123
push!(solver_tracker, prob_name => successful_solvers);
103124

104125
# Handles the non-general cases.
@@ -161,12 +182,11 @@ function check_solver(prob, solver)
161182
return true
162183
end
163184

164-
# Adds an additional, selected, solver to the general solver set.
165185
# Adds an additional, selected, solver to the general solver set.
166186
function add_solver!(solvers_general, selected_solver_name, additional_solver_set, wp)
167187
if isnothing(selected_solver_name)
168188
isempty(wp.wps) && return
169-
selected_idx = argmin(median.(getfield.(wp.wps, :times)))
189+
selected_idx = argmin(median.(getfield.(wp.wps, :times)))
170190
else
171191
selected_idx = findfirst(s -> s.name==selected_solver_name, additional_solver_set)
172192
isnothing(selected_solver) && error("The $(selected_solver_name) was designated to \

0 commit comments

Comments
 (0)