-
-
Notifications
You must be signed in to change notification settings - Fork 78
/
precs.jl
122 lines (105 loc) · 4.81 KB
/
precs.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
using Sundials, Test, ModelingToolkit, LinearAlgebra, IncompleteLU
import AlgebraicMultigrid
println("precs.jl start tests...")
const N = 32
const xyd_brusselator = range(0; stop = 1, length = N)
brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0
limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a
function brusselator_2d_loop(du, u, p, t)
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, t)
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.0, 10.0, 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_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
u0, (0.0, 11.5), p)
println("rob_ode_brusselator_2d_mtk")
@time prob_ode_brusselator_2d_mtk = ODEProblem(modelingtoolkitize(prob_ode_brusselator_2d), [],
(0.0, 11.5); jac = true, sparse = true);
u0 = prob_ode_brusselator_2d_mtk.u0
p = prob_ode_brusselator_2d_mtk.p
const jaccache = prob_ode_brusselator_2d_mtk.f.jac(u0, p, 0.0)
const W = I - 1.0 * jaccache
prectmp = ilu(W; τ = 50.0)
const preccache = Ref(prectmp)
function psetupilu(p, t, u, du, jok, jcurPtr, gamma)
if jok
prob_ode_brusselator_2d_mtk.f.jac(jaccache, u, p, t)
jcurPtr[] = true
# W = I - gamma*J
@. W = -gamma * jaccache
idxs = diagind(W)
@. @view(W[idxs]) = @view(W[idxs]) + 1
# Build preconditioner on W
preccache[] = ilu(W; τ = 5.0)
end
end
function precilu(z, r, p, t, y, fy, gamma, delta, lr)
ldiv!(z, preccache[], r)
end
println("aspreconditioner")
@time prectmp2 = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(W;
presmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
1))),
postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
1)))))
const preccache2 = Ref(prectmp2)
function psetupamg(p, t, u, du, jok, jcurPtr, gamma)
if jok
prob_ode_brusselator_2d_mtk.f.jac(jaccache, u, p, t)
jcurPtr[] = true
# W = I - gamma*J
@. W = -gamma * jaccache
idxs = diagind(W)
@. @view(W[idxs]) = @view(W[idxs]) + 1
# Build preconditioner on W
preccache2[] = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(W;
presmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
1))),
postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
1)))))
end
end
function precamg(z, r, p, t, y, fy, gamma, delta, lr)
ldiv!(z, preccache2[], r)
end
#=
# Remove for testing
println("sol1:")
@time sol1 = solve(prob_ode_brusselator_2d, CVODE_BDF(; linear_solver = :GMRES);
save_everystep = false);
println("sol2:")
@time sol2 = solve(prob_ode_brusselator_2d,
CVODE_BDF(; linear_solver = :GMRES, prec = precilu, psetup = psetupilu,
prec_side = 1); save_everystep = false);
println("sol3:")
@time sol3 = solve(prob_ode_brusselator_2d,
CVODE_BDF(; linear_solver = :GMRES, prec = precamg, psetup = psetupamg,
prec_side = 1); save_everystep = false);
@test sol1.destats.nf > sol2.destats.nf
@test sol1.destats.nf > sol3.destats.nf
=#