@@ -5,48 +5,48 @@ using KrylovPreconditioners
55using SparseArrays, LinearAlgebra
66
77function bratu! (res, y, Δx, λ)
8- N = length (y)
9- for i in 1 : N
10- y_l = i == 1 ? zero (eltype (y)) : y[i - 1 ]
11- y_r = i == N ? zero (eltype (y)) : y[i + 1 ]
12- y′′ = (y_r - 2 y[i] + y_l) / Δx^ 2
13-
14- res[i] = y′′ + λ * exp (y[i]) # = 0
15- end
16- nothing
8+ N = length (y)
9+ for i in 1 : N
10+ y_l = i == 1 ? zero (eltype (y)) : y[i - 1 ]
11+ y_r = i == N ? zero (eltype (y)) : y[i + 1 ]
12+ y′′ = (y_r - 2 y[i] + y_l) / Δx^ 2
13+
14+ res[i] = y′′ + λ * exp (y[i]) # = 0
15+ end
16+ return nothing
1717end
1818
1919function bratu (y, dx, λ)
20- res = similar (y)
21- bratu! (res, y, dx, λ)
22- res
20+ res = similar (y)
21+ bratu! (res, y, dx, λ)
22+ return res
2323end
2424
2525function true_sol_bratu (x)
26- # for λ = 3.51382, 2nd sol θ = 4.8057
27- θ = 4.79173
28- - 2 * log (cosh (θ * (x- 0.5 )/ 2 ) / (cosh (θ/ 4 )))
26+ # for λ = 3.51382, 2nd sol θ = 4.8057
27+ θ = 4.79173
28+ return - 2 * log (cosh (θ * (x - 0.5 ) / 2 ) / (cosh (θ / 4 )))
2929end
3030
3131const N = 10_000
3232const λ = 3.51382
3333const dx = 1 / (N + 1 ) # Grid-spacing
3434
35- x = LinRange (0.0 + dx, 1.0 - dx, N)
36- u₀ = sin .(x.* π)
35+ x = LinRange (0.0 + dx, 1.0 - dx, N)
36+ u₀ = sin .(x .* π)
3737
3838reference = true_sol_bratu .(x)
3939
4040uₖ_1 = newton_krylov! (
41- (res, u) -> bratu! (res, u, dx, λ),
42- copy (u₀), similar (u₀);
43- Solver = CgSolver,
41+ (res, u) -> bratu! (res, u, dx, λ),
42+ copy (u₀), similar (u₀);
43+ Solver = CgSolver,
4444)
4545
4646uₖ_2 = newton_krylov (
47- (u) -> bratu (u, dx, λ),
48- copy (u₀);
49- Solver = CgSolver
47+ (u) -> bratu (u, dx, λ),
48+ copy (u₀);
49+ Solver = CgSolver
5050)
5151
5252ϵ1 = abs2 .(uₖ_1 .- reference)
@@ -55,19 +55,19 @@ uₖ_2 = newton_krylov(
5555# #
5656# Solving with a fixed forcing
5757newton_krylov! (
58- (res, u) -> bratu! (res, u, dx, λ),
59- copy (u₀), similar (u₀);
60- Solver = CgSolver,
61- forcing = NewtonKrylov. Fixed ()
58+ (res, u) -> bratu! (res, u, dx, λ),
59+ copy (u₀), similar (u₀);
60+ Solver = CgSolver,
61+ forcing = NewtonKrylov. Fixed ()
6262)
6363
6464# #
6565# Solving with no forcing
6666newton_krylov! (
67- (res, u) -> bratu! (res, u, dx, λ),
68- copy (u₀), similar (u₀);
69- Solver = CgSolver,
70- forcing = nothing
67+ (res, u) -> bratu! (res, u, dx, λ),
68+ copy (u₀), similar (u₀);
69+ Solver = CgSolver,
70+ forcing = nothing
7171)
7272
7373# #
@@ -81,40 +81,40 @@ newton_krylov!(
8181# #
8282# Solve using GMRES + ILU Preconditoner
8383@time newton_krylov! (
84- (res, u) -> bratu! (res, u, dx, λ),
85- copy (u₀), similar (u₀);
86- Solver = GmresSolver,
87- N = (J) -> ilu (collect (J)), # Assembles the full Jacobian
88- ldiv = true ,
84+ (res, u) -> bratu! (res, u, dx, λ),
85+ copy (u₀), similar (u₀);
86+ Solver = GmresSolver,
87+ N = (J) -> ilu (collect (J)), # Assembles the full Jacobian
88+ ldiv = true ,
8989)
9090
9191# #
9292# Solve using FGMRES + ILU Preconditoner
9393@time newton_krylov! (
94- (res, u) -> bratu! (res, u, dx, λ),
95- copy (u₀), similar (u₀);
96- Solver = FgmresSolver,
97- N = (J) -> ilu (collect (J)), # Assembles the full Jacobian
98- ldiv = true ,
94+ (res, u) -> bratu! (res, u, dx, λ),
95+ copy (u₀), similar (u₀);
96+ Solver = FgmresSolver,
97+ N = (J) -> ilu (collect (J)), # Assembles the full Jacobian
98+ ldiv = true ,
9999)
100100
101101# #
102102# Solve using FGMRES + GMRES Preconditoner
103103struct GmresPreconditioner{JOp}
104- J:: JOp
105- itmax:: Int
104+ J:: JOp
105+ itmax:: Int
106106end
107107
108108function LinearAlgebra. mul! (y, P:: GmresPreconditioner , x)
109- sol, _ = gmres (P. J, x; P. itmax)
110- copyto! (y, sol)
109+ sol, _ = gmres (P. J, x; P. itmax)
110+ return copyto! (y, sol)
111111end
112112
113113@time newton_krylov! (
114- (res, u) -> bratu! (res, u, dx, λ),
115- copy (u₀), similar (u₀);
116- Solver = FgmresSolver,
117- N = (J) -> GmresPreconditioner (J, 30 ),
114+ (res, u) -> bratu! (res, u, dx, λ),
115+ copy (u₀), similar (u₀);
116+ Solver = FgmresSolver,
117+ N = (J) -> GmresPreconditioner (J, 30 ),
118118)
119119
120120# # Explodes..
131131# verbose = 1,
132132# Solver = BicgstabSolver,
133133# η_max = nothing
134- # )
134+ # )
0 commit comments