Add KA example (#6)

Author: vchuravy

examples/Project.toml

@@ -1,10 +1,11 @@
 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 KrylovPreconditioners = "45d422c2-293f-44ce-8315-2cb988662dec"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NewtonKrylov = "0be81120-40bf-4f8b-adf0-26103efb66f1"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
-[sources]
-NewtonKrylov = {path = ".."}
+[sources.NewtonKrylov]
+path = ".."

examples/bratu_ka.jl

@@ -0,0 +1,85 @@
+# # 1D bratu equation from (Kan2022-ko)[@cite]
+
+# ## Necessary packages
+using NewtonKrylov, Krylov
+using KrylovPreconditioners
+using SparseArrays, LinearAlgebra
+using CairoMakie
+using KernelAbstractions
+
+# ## 1D Bratu equations
+# $y′′ + λ * exp(y) = 0$
+
+@kernel function bratu_kernel!(res, y, Δx, λ)
+    i = @index(Global, Linear)
+    N = length(res)
+    y_l = i == 1 ? zero(eltype(y)) : y[i - 1]
+    y_r = i == N ? zero(eltype(y)) : y[i + 1]
+    y′′ = (y_r - 2y[i] + y_l) / Δx^2
+
+    res[i] = y′′ + λ * exp(y[i]) # = 0
+end
+
+function bratu!(res, y, Δx, λ)
+    device = KernelAbstractions.get_backend(res)
+    kernel = bratu_kernel!(device)
+    kernel(res, y, Δx, λ, ndrange = length(res))
+    return nothing
+end
+
+function bratu(y, dx, λ)
+    res = similar(y)
+    bratu!(res, y, dx, λ)
+    return res
+end
+
+# ## Reference solution
+function true_sol_bratu(x)
+    ## for λ = 3.51382, 2nd sol θ = 4.8057
+    θ = 4.79173
+    return -2 * log(cosh(θ * (x - 0.5) / 2) / (cosh(θ / 4)))
+end
+
+# ## Choice of parameters
+const N = 10_000
+const λ = 3.51382
+const dx = 1 / (N + 1) # Grid-spacing
+
+# ### Domain and Inital condition
+x = LinRange(0.0 + dx, 1.0 - dx, N)
+u₀ = sin.(x .* π)
+
+lines(x, u₀, label = "Inital guess")
+
+# ## Reference solution evaluated over domain
+reference = true_sol_bratu.(x)
+
+fig, ax = lines(x, u₀, label = "Inital guess")
+lines!(ax, x, reference, label = "Reference solution")
+axislegend(ax, position = :cb)
+fig
+
+# ## Solving using inplace variant and CG
+uₖ, _ = newton_krylov!(
+    (res, u) -> bratu!(res, u, dx, λ),
+    copy(u₀), similar(u₀);
+    Solver = CgSolver,
+)
+
+ϵ = abs2.(uₖ .- reference)
+
+let
+    fig = Figure(size = (800, 800))
+    ax = Axis(fig[1, 1], title = "", ylabel = "", xlabel = "")
+
+    lines!(ax, x, reference, label = "True solution")
+    lines!(ax, x, u₀, label = "Initial guess")
+    lines!(ax, x, uₖ, label = "Newton-Krylov solution")
+
+    axislegend(ax, position = :cb)
+
+    ax = Axis(fig[1, 2], title = "Error", ylabel = "abs2 err", xlabel = "")
+    lines!(ax, abs2.(uₖ .- reference))
+
+    fig
+end