Skip to content

Commit 457acf6

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #1354 from GodotMisogi/dev-pde-benchmarks
Simple Handwritten PDEs: Benchmarks with preconditioners for Krylov subspace methods + minor typo fixes
2 parents 1c9d66d + bbd9d5e commit 457acf6

File tree

8 files changed

+383
-121
lines changed

8 files changed

+383
-121
lines changed

‎benchmarks/SimpleHandwrittenPDE/Manifest.toml‎

Lines changed: 144 additions & 109 deletions
Large diffs are not rendered by default.

‎benchmarks/SimpleHandwrittenPDE/Project.toml‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,14 @@
11
[deps]
2+
AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
23
ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
34
DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
5+
IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
46
LSODA = "7f56f5a3-f504-529b-bc02-0b1fe5e64312"
57
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
68
LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
79
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
810
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
11+
Preconditioners = "af69fa37-3177-5a40-98ee-561f696e4fcd"
912
SciMLBenchmarks = "31c91b34-3c75-11e9-0341-95557aab0344"
1013
SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
1114
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -15,8 +18,8 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
1518
[compat]
1619
ClassicalOrthogonalPolynomials = "0.15"
1720
DiffEqDevTools = "2.22"
18-
LinearSolve = "3"
1921
LSODA = "0.6, 0.7"
22+
LinearSolve = "3"
2023
OrdinaryDiffEq = "6"
2124
Plots = "1.4"
2225
SciMLBenchmarks = "0.1"

‎benchmarks/SimpleHandwrittenPDE/allen_cahn_fdm_wpd.jmd‎

Lines changed: 113 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -223,6 +223,119 @@ labels = hcat(
223223
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Krylov Linsolve, Low Tolerances")
224224
```
225225

226+
Krylov solvers with preconditioners.
227+
```julia
228+
# %%
229+
# Weighted diagonal preconditioner
230+
import LinearAlgebra as LA
231+
Base.@kwdef struct WeightedDiagonalPreconBuilder
232+
w::Float64
233+
end
234+
235+
(builder::WeightedDiagonalPreconBuilder)(A, du, u, p, t, newW, Plprev, Prprev, solverdata) = (builder.w * LA.Diagonal(convert(AbstractMatrix, A)), LA.I)
236+
237+
# Incomplete LU factorization
238+
using IncompleteLU
239+
function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
240+
if newW === nothing || newW
241+
Pl = ilu(convert(AbstractMatrix, W), τ = 50.0)
242+
else
243+
Pl = Plprev
244+
end
245+
Pl, nothing
246+
end
247+
248+
# Algebraic multigrid
249+
using AlgebraicMultigrid
250+
function algebraicmultigrid(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
251+
if newW === nothing || newW
252+
Pl = aspreconditioner(ruge_stuben(convert(AbstractMatrix, W)))
253+
else
254+
Pl = Plprev
255+
end
256+
Pl, nothing
257+
end
258+
function algebraicmultigrid2(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
259+
if newW === nothing || newW
260+
A = convert(AbstractMatrix, W)
261+
Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A,
262+
presmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
263+
1))),
264+
postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
265+
1)))))
266+
else
267+
Pl = Plprev
268+
end
269+
Pl, nothing
270+
end
271+
272+
# Aliases
273+
using Preconditioners
274+
function cholesky(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
275+
W = convert(AbstractMatrix, W)
276+
if newW === nothing || newW
277+
Pl = CholeskyPreconditioner(W)
278+
else
279+
Pl = Plprev
280+
end
281+
Pl, nothing
282+
end
283+
chol = cholesky
284+
inc_lu = incompletelu
285+
w_diag = WeightedDiagonalPreconBuilder(w = 0.9)
286+
amg = algebraicmultigrid
287+
amg2 = algebraicmultigrid2
288+
289+
abstols = 0.1 .^ (8:13)
290+
reltols = 0.1 .^ (5:10)
291+
N = length(abstols)
292+
setups = [
293+
# Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES(), :dts => 1e-2 * ones(N), :adaptive => true),
294+
# Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES(), :dts => 1e-2 * ones(N), :adaptive => true),
295+
# Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES(), :dts => 1e-2 * ones(N), :adaptive => true),
296+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES(), precs = inc_lu, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
297+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES(), precs = inc_lu, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
298+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES(), precs = inc_lu, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
299+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES(), precs = w_diag, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
300+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES(), precs = w_diag, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
301+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES(), precs = w_diag, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
302+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES(), precs = amg, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
303+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES(), precs = amg, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
304+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES(), precs = amg, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
305+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES(), precs = amg2, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
306+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES(), precs = amg2, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
307+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES(), precs = amg2, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
308+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES()), :dts => 1e-2 * ones(N), :adaptive => true),
309+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES()), :dts => 1e-2 * ones(N), :adaptive => true),
310+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES()), :dts => 1e-2 * ones(N), :adaptive => true),
311+
]
312+
labels = hcat(
313+
# "KenCarp3 (Cholesky)",
314+
# "KenCarp4 (Cholesky)",
315+
# "KenCarp5 (Cholesky)",
316+
"KenCarp3 (ILU)",
317+
"KenCarp4 (ILU)",
318+
"KenCarp5 (ILU)",
319+
"KenCarp3 (Diagonal, w = $(w_diag.w))",
320+
"KenCarp4 (Diagonal, w = $(w_diag.w))",
321+
"KenCarp5 (Diagonal, w = $(w_diag.w))",
322+
"KenCarp3 (AMG)",
323+
"KenCarp4 (AMG)",
324+
"KenCarp5 (AMG)",
325+
"KenCarp3 (AMG-Jacobi)",
326+
"KenCarp4 (AMG-Jacobi)",
327+
"KenCarp5 (AMG-Jacobi)",
328+
"KenCarp3 (Identity)",
329+
"KenCarp4 (Identity)",
330+
"KenCarp5 (Identity)",
331+
)
332+
@time wp = WorkPrecisionSet(prob, abstols, reltols, setups;
333+
print_names=true, names=labels, numruns=5, error_estimate=:l2,
334+
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
335+
336+
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Krylov Linsolve, Low Tolerances")
337+
```
338+
226339
#### Exponential Integrators
227340

228341
```julia

‎benchmarks/SimpleHandwrittenPDE/burgers_fdm_wpd.jmd‎

Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -223,6 +223,123 @@ labels = hcat(
223223
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Krylov Linsolve, Low Tolerances")
224224
```
225225

226+
Krylov solvers with preconditioners.
227+
```julia
228+
# Weighted diagonal preconditioner
229+
import LinearAlgebra as LA
230+
Base.@kwdef struct WeightedDiagonalPreconBuilder
231+
w::Float64
232+
end
233+
234+
(builder::WeightedDiagonalPreconBuilder)(A, du, u, p, t, newW, Plprev, Prprev, solverdata) = (builder.w * LA.Diagonal(convert(AbstractMatrix, A)), LA.I)
235+
236+
# Incomplete LU factorization
237+
using IncompleteLU
238+
Base.@kwdef struct IncompleteLUPreconBuilder
239+
τ::Float64
240+
end
241+
function (builder::IncompleteLUPreconBuilder)(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
242+
if newW === nothing || newW
243+
Pl = ilu(convert(AbstractMatrix, W), τ = builder.τ)
244+
else
245+
Pl = Plprev
246+
end
247+
Pl, nothing
248+
end
249+
250+
# Algebraic multigrid
251+
using AlgebraicMultigrid
252+
function algebraicmultigrid(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
253+
if newW === nothing || newW
254+
Pl = aspreconditioner(ruge_stuben(convert(AbstractMatrix, W)))
255+
else
256+
Pl = Plprev
257+
end
258+
Pl, nothing
259+
end
260+
function algebraicmultigrid2(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
261+
if newW === nothing || newW
262+
A = convert(AbstractMatrix, W)
263+
Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A,
264+
presmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
265+
1))),
266+
postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
267+
1)))))
268+
else
269+
Pl = Plprev
270+
end
271+
Pl, nothing
272+
end
273+
274+
# Cholesky
275+
using Preconditioners
276+
function cholesky(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
277+
W = convert(AbstractMatrix, W)
278+
if newW === nothing || newW
279+
Pl = CholeskyPreconditioner(W)
280+
else
281+
Pl = Plprev
282+
end
283+
Pl, nothing
284+
end
285+
286+
# Aliases
287+
chol = cholesky
288+
inc_lu = IncompleteLUPreconBuilder(τ = 50.0)
289+
w_diag = WeightedDiagonalPreconBuilder(w = 0.9)
290+
amg = algebraicmultigrid
291+
amg2 = algebraicmultigrid2
292+
293+
abstols = 0.1 .^ (8:13)
294+
reltols = 0.1 .^ (5:10)
295+
N = length(abstols)
296+
setups = [
297+
# Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES(), :dts => 1e-2 * ones(N), :adaptive => true),
298+
# Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES(), :dts => 1e-2 * ones(N), :adaptive => true),
299+
# Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES(), :dts => 1e-2 * ones(N), :adaptive => true),
300+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES(), precs = inc_lu, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
301+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES(), precs = inc_lu, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
302+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES(), precs = inc_lu, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
303+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES(), precs = w_diag, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
304+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES(), precs = w_diag, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
305+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES(), precs = w_diag, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
306+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES(), precs = amg, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
307+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES(), precs = amg, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
308+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES(), precs = amg, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
309+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES(), precs = amg2, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
310+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES(), precs = amg2, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
311+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES(), precs = amg2, concrete_jac = true), :dts => 1e-2 * ones(N), :adaptive => true),
312+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES()), :dts => 1e-2 * ones(N), :adaptive => true),
313+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES()), :dts => 1e-2 * ones(N), :adaptive => true),
314+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES()), :dts => 1e-2 * ones(N), :adaptive => true),
315+
]
316+
labels = hcat(
317+
# "KenCarp3 (Cholesky)",
318+
# "KenCarp4 (Cholesky)",
319+
# "KenCarp5 (Cholesky)",
320+
"KenCarp3 (ILU, τ = $(inc_lu.τ))",
321+
"KenCarp4 (ILU, τ = $(inc_lu.τ))",
322+
"KenCarp5 (ILU, τ = $(inc_lu.τ))",
323+
"KenCarp3 (Diagonal, w = $(w_diag.w))",
324+
"KenCarp4 (Diagonal, w = $(w_diag.w))",
325+
"KenCarp5 (Diagonal, w = $(w_diag.w))",
326+
"KenCarp3 (AMG)",
327+
"KenCarp4 (AMG)",
328+
"KenCarp5 (AMG)",
329+
"KenCarp3 (AMG-Jacobi)",
330+
"KenCarp4 (AMG-Jacobi)",
331+
"KenCarp5 (AMG-Jacobi)",
332+
"KenCarp3 (Identity)",
333+
"KenCarp4 (Identity)",
334+
"KenCarp5 (Identity)",
335+
)
336+
@time wp = WorkPrecisionSet(prob, abstols, reltols, setups;
337+
print_names=true, names=labels, numruns=5, error_estimate=:l2,
338+
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
339+
340+
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Krylov Linsolve, Low Tolerances")
341+
```
342+
226343
#### Exponential Integrators
227344

228345
```julia

‎benchmarks/SimpleHandwrittenPDE/kdv_fdm_wpd.jmd‎

Lines changed: 2 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ author: Arjit Seth
55

66
### Problem Description
77

8-
The Kortwrieg de Vries (KdV) partial differential equation is solved on the domain $[-L, L] \times [0, T] \in \mathbb R \times R,~L = 16,~T = 20$, with the following initial and periodic boundary conditions:
8+
The Kortwrieg de Vries (KdV) partial differential equation is solved on the domain $[-L, L] \times [0, T] \in \mathbb R \times R,~L = 16,~T = 5$, with the following initial and periodic boundary conditions:
99
```math
1010
\begin{align}
1111
\partial_t u(t,x) & = -6u(t,x)\partial_x u(t,x) -\partial_x^3 u(t,x), \\
@@ -177,9 +177,6 @@ setups = [
177177
Dict(:alg => KenCarp3()),
178178
Dict(:alg => KenCarp4()),
179179
Dict(:alg => KenCarp5()),
180-
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES())),
181-
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES())),
182-
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES())),
183180
Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:GMRES)),
184181
Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:GMRES)),
185182
Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES)),
@@ -188,9 +185,6 @@ labels = hcat(
188185
"KenCarp3 (dense)",
189186
"KenCarp4 (dense)",
190187
"KenCarp5 (dense)",
191-
"KenCarp3 (Krylov)",
192-
"KenCarp4 (Krylov)",
193-
"KenCarp5 (Krylov)",
194188
"ARKODE3 (Krylov)",
195189
"ARKODE4 (Krylov)",
196190
"ARKODE5 (Krylov)",
@@ -199,7 +193,7 @@ labels = hcat(
199193
print_names=true, names=labels, numruns=5, error_estimate=:l2,
200194
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
201195

202-
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Krylov Linsolve, Low Tolerances")
196+
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Low Tolerances")
203197
```
204198

205199
#### Exponential Integrators

‎benchmarks/SimpleHandwrittenPDE/kdv_spectral_wpd.jmd‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -194,7 +194,7 @@ labels = hcat(
194194
print_names=true, names=labels, numruns=5, error_estimate=:l2,
195195
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
196196

197-
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Krylov Linsolve, Low Tolerances")
197+
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Low Tolerances")
198198
```
199199

200200
#### Exponential Integrators

‎benchmarks/SimpleHandwrittenPDE/ks_fdm_wpd.jmd‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -197,7 +197,7 @@ labels = hcat(
197197
print_names=true, names=labels, numruns=5, error_estimate=:l2,
198198
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
199199

200-
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Krylov Linsolve, Low Tolerances")
200+
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Low Tolerances")
201201
```
202202

203203
#### Exponential Integrators

‎benchmarks/SimpleHandwrittenPDE/ks_spectral_wpd.jmd‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -194,7 +194,7 @@ labels = hcat(
194194
print_names=true, names=labels, numruns=5, error_estimate=:l2,
195195
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
196196

197-
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Krylov Linsolve, Low Tolerances")
197+
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Low Tolerances")
198198
```
199199

200200
#### Exponential Integrators

0 commit comments

Comments
 (0)