@@ -150,8 +150,7 @@ reltols = 1.0 ./ 10.0 .^ (5:8);
150150
151151## Failures
152152
153- Before proceeding to the results, we note the notable omissions. CVODE with KLU diverges in the solution, and
154- thus it is omitted from the results:
153+ CVODE with KLU diverges on this problem and is omitted from all suites:
155154
156155```julia
157156try
@@ -161,78 +160,37 @@ catch e
161160end
162161```
163162
164- ## Work-Precision Diagrams (CVODE and lsoda solvers )
163+ ## Work-Precision Diagrams (competitive CVODE / lsoda)
165164
166- #### Declare solvers.
165+ Main suite: methods that remain competitive on this large sparse stiff system.
166+ Dense Lapack and bare GMRES (no preconditioner) are timed only in the loser
167+ section below.
167168
168169```julia
169170setups = [
170171 Dict(:alg=>lsoda(), :prob_choice => 1),
171172 Dict(:alg=>CVODE_BDF(), :prob_choice => 1),
172- Dict(:alg=>CVODE_BDF(linear_solver = :LapackDense), :prob_choice => 1),
173- Dict(:alg=>CVODE_BDF(linear_solver = :GMRES), :prob_choice => 1),
174173 Dict(
175174 :alg=>CVODE_BDF(linear_solver = :GMRES, prec = precilu, psetup = psetupilu, prec_side = 1),
176175 :prob_choice => 2)
177176];
178177```
179178
180- #### Plot Work-Precision Diagram.
181-
182179```julia
183180wp = WorkPrecisionSet(
184181 [oprob, oprob_sparse, sparsejacprob], abstols, reltols, setups; error_estimate = :l2,
185182 saveat = tf/10000.0, appxsol = [test_sol, test_sol, test_sol], maxiters = Int(1e6), numruns = 1)
186183
187- names = ["lsoda" "CVODE_BDF" "CVODE_BDF (LapackDense)" "CVODE_BDF (GMRES)" "CVODE_BDF (GMRES, iLU)" "CVODE_BDF (KLU, sparse jac)"]
188- plot(wp; label = names)
189- ```
190-
191- ## Work-Precision Diagrams (various Julia solvers)
192-
193- #### Declare solvers (using default linear solver).
194-
195- ```julia
196- setups = [
197- Dict(:alg=>TRBDF2(autodiff = AutoFiniteDiff())),
198- Dict(:alg=>QNDF(autodiff = AutoFiniteDiff())),
199- Dict(:alg=>FBDF(autodiff = AutoFiniteDiff())),
200- Dict(:alg=>KenCarp4(autodiff = AutoFiniteDiff()))
201- ];
202- ```
203-
204- #### Plot Work-Precision Diagram (using default linear solver).
205-
206- ```julia
207- wp = WorkPrecisionSet(oprob, abstols, reltols, setups; error_estimate = :l2,
208- saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e6), numruns = 1)
209-
210- names = ["TRBDF2" "QNDF" "FBDF" "KenCarp4"]
184+ names = ["lsoda" "CVODE_BDF" "CVODE_BDF (GMRES, iLU)"]
211185plot(wp; label = names)
212186```
213187
214- #### Declare solvers (using GMRES linear solver).
215-
216- ```julia
217- setups = [
218- Dict(:alg=>TRBDF2(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff())),
219- Dict(:alg=>QNDF(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff())),
220- Dict(:alg=>FBDF(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff())),
221- Dict(:alg=>KenCarp4(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()))
222- ];
223- ```
224-
225- #### Plot Work-Precision Diagram (using GMRES linear solver).
188+ ## Work-Precision Diagrams (competitive Julia solvers)
226189
227- ```julia
228- wp = WorkPrecisionSet(oprob, abstols, reltols, setups; error_estimate = :l2,
229- saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e6), numruns = 1)
230-
231- names = ["TRBDF2 (GMRES)" "QNDF (GMRES)" "FBDF (GMRES)" "KenCarp4 (GMRES)"]
232- plot(wp; label = names)
233- ```
190+ Only preconditioned GMRES and sparse KLU variants; default dense factorizations
191+ and unpreconditioned GMRES are in the loser section.
234192
235- #### Declare solvers (using GMRES linear solver, with pre-conditioner).
193+ #### GMRES + incomplete LU
236194
237195```julia
238196setups = [
@@ -245,8 +203,6 @@ setups = [
245203];
246204```
247205
248- #### Plot Work-Precision Diagram (using GMRES linear solver, with pre-conditioner).
249-
250206```julia
251207wp = WorkPrecisionSet(sparsejacprob, abstols, reltols, setups; error_estimate = :l2,
252208 saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e6), numruns = 1)
@@ -255,9 +211,7 @@ names = ["TRBDF2 (GMRES, iLU)" "QNDF (GMRES, iLU)" "FBDF (GMRES, iLU)" "KenCarp4
255211plot(wp; label = names)
256212```
257213
258- #### Declare solvers (using sparse jacobian)
259-
260- We designate the solvers we wish to use.
214+ #### Sparse Jacobian + KLU
261215
262216```julia
263217setups = [
@@ -268,10 +222,6 @@ setups = [
268222];
269223```
270224
271- #### Plot Work-Precision Diagram (using sparse jacobian)
272-
273- Finally, we generate a work-precision diagram for the selection of solvers.
274-
275225```julia
276226wp = WorkPrecisionSet(sparsejacprob, abstols, reltols, setups; error_estimate = :l2,
277227 saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e6), numruns = 1)
@@ -280,6 +230,69 @@ names = ["TRBDF2 (KLU, sparse jac)" "QNDF (KLU, sparse jac)" "FBDF (KLU, sparse
280230plot(wp; label = names)
281231```
282232
233+ ## Loser methods (large cost in isolation)
234+
235+ On this ~1122-ODE sparse chemistry system, the following are not competitive:
236+ dense Lapack factorization, default dense Julia Newton/linear solves, and GMRES
237+ without a preconditioner. We do **not** fold them into the main work-precision
238+ suites. Instead each is timed **once**, in isolation, at a fixed tolerance, next
239+ to a competitive sparse reference so the wall-time gap is obvious.
240+
241+ ```julia
242+ const _loser_tol = 1e-6
243+ const _loser_maxiters = Int(1e6)
244+ _solve_kwargs = (; abstol = _loser_tol, reltol = _loser_tol, maxiters = _loser_maxiters,
245+ save_everystep = false)
246+
247+ loser_labels = String[]
248+ loser_elapsed = Float64[]
249+
250+ function _time_loser!(label, prob, alg)
251+ println("--- $label ---")
252+ t = @elapsed sol = solve(prob, alg; _solve_kwargs...)
253+ @show sol.retcode
254+ println("elapsed = ", t, " s")
255+ push!(loser_labels, label)
256+ push!(loser_elapsed, t)
257+ return sol
258+ end
259+
260+ # Competitive reference (sparse KLU)
261+ _time_loser!("FBDF + KLU (reference)", sparsejacprob,
262+ FBDF(linsolve = KLUFactorization(), autodiff = AutoFiniteDiff()))
263+
264+ # Dense CVODE Lapack
265+ _time_loser!("CVODE_BDF LapackDense", oprob, CVODE_BDF(linear_solver = :LapackDense))
266+
267+ # Bare CVODE GMRES (no preconditioner)
268+ _time_loser!("CVODE_BDF GMRES (no prec)", oprob, CVODE_BDF(linear_solver = :GMRES))
269+
270+ # Default dense Julia factorizations on the non-sparse problem
271+ _time_loser!("TRBDF2 (default dense)", oprob, TRBDF2(autodiff = AutoFiniteDiff()))
272+ _time_loser!("QNDF (default dense)", oprob, QNDF(autodiff = AutoFiniteDiff()))
273+ _time_loser!("FBDF (default dense)", oprob, FBDF(autodiff = AutoFiniteDiff()))
274+ _time_loser!("KenCarp4 (default dense)", oprob, KenCarp4(autodiff = AutoFiniteDiff()))
275+
276+ # Unpreconditioned GMRES on the dense residual problem
277+ _time_loser!("TRBDF2 GMRES (no prec)", oprob,
278+ TRBDF2(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()))
279+ _time_loser!("QNDF GMRES (no prec)", oprob,
280+ QNDF(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()))
281+ _time_loser!("FBDF GMRES (no prec)", oprob,
282+ FBDF(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()))
283+ _time_loser!("KenCarp4 GMRES (no prec)", oprob,
284+ KenCarp4(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()))
285+ ```
286+
287+ ```julia
288+ # Relative cost vs the sparse KLU reference (first entry)
289+ ref_t = loser_elapsed[1]
290+ bar(loser_labels, loser_elapsed ./ ref_t; xrotation = 45, legend = false,
291+ ylabel = "wall time / (FBDF+KLU reference)",
292+ title = "BCR loser isolation (tol=$_loser_tol, one solve each)",
293+ size = (900, 500), left_margin = 5Plots.mm, bottom_margin = 15Plots.mm)
294+ ```
295+
283296## Summary of results
284297
285298Finally, we compute a single diagram comparing the various solvers used.
0 commit comments