Resolving Issue: Errors in Simple Handwritten PDEs as ODEs

benchmarks/SimpleHandwrittenPDE/Project.toml

@@ -1,30 +1,23 @@
 [deps]
 ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f"
+ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 LSODA = "7f56f5a3-f504-529b-bc02-0b1fe5e64312"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
-ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-ODEInterfaceDiffEq = "09606e27-ecf5-54fc-bb29-004bd9f985bf"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
 SciMLBenchmarks = "31c91b34-3c75-11e9-0341-95557aab0344"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 
 [compat]
 ApproxFun = "0.13"
 DiffEqDevTools = "2.22"
-DomainSets = "0.6, 0.7"
 LSODA = "0.6, 0.7"
-LinearSolve = "2, 3"
-MethodOfLines = "0.9"
-ModelingToolkit = "8, 10"
-ODEInterfaceDiffEq = "3.7"
-OrdinaryDiffEq = "5.41, 6"
 Plots = "1.4"
-RecursiveFactorization = "0.2.23"
 SciMLBenchmarks = "0.1"
 Sundials = "4.2"

benchmarks/SimpleHandwrittenPDE/allen_cahn_fdm_wpd.jmd

@@ -1,233 +1,130 @@
 ---
-title: Allen_Cahn FDM Work-Precision Diagrams
-author: HAO HAO
+title: Allen-Cahn Finite-Difference Method Work-Precision Diagrams
+author: Arjit Seth
 ---
 
-```julia
-using ApproxFun, OrdinaryDiffEq, Sundials, LinearSolve
-using DiffEqDevTools
-using LinearAlgebra, RecursiveFactorization
-using Plots; gr()
-```
-
-Here is the Burgers equation using FDM.
-
-```julia
-# Define the linear and nonlinear terms
-function lin_term(N)
-    dx = 1/(N + 1)
-    du = 1/2 * ones(N - 1) # super diagonal
-    dl = -1/2 * ones(N - 1) # lower diagonal
-    DiffEqArrayOperator(0.01*(1/dx) * diagm(-1 => dl, 1 => du))
-end
-
-function nl_term(N)
-    function (du,u,p,t)
-        du .= u .- u.^3
-    end
-end
-
-# Construct the problem
-function allen_cahn(N)
-    f1 = lin_term(N)
-    f2 = nl_term(N)
-    dx = 1 / (N + 1)
-    xs = (1:N) * dx
-
-    f0 = x -> .53*x + .47*sin(-1.5*pi*x) - x
-    u0 = f0.(xs)
-    prob = SplitODEProblem(f1, f2, u0, (0.0, 1.0))
-    xs, prob
-end;
-```
-
-Reference solution using Vern9 is below:
-
-```julia
-xs, prob = allen_cahn(100)
-sol = solve(prob, RadauIIA5(autodiff=false); abstol=1e-14, reltol=1e-14, dt=1e-4, adaptive=false)
-test_sol = TestSolution(sol);
-
-tslices = [0.0 0.25 0.50 0.75 1.]
-ys = hcat((sol(t) for t in tslices)...)
-labels = ["t = $t" for t in tslices]
-plot(xs, ys, label=labels)
-```
-
-Linear solvers
-
-```julia
-const LS_Dense = LinSolveFactorize(lu)
-```
-
-## High tolerances
+## Equations
 
-## In-family comparisons
-
-1.IMEX methods (dense linear solver)
-
-```julia
-abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
-reltols = 0.1 .^ (1:4)
-multipliers = 0.5 .^ (0:3)
-setups = [Dict(:alg => IMEXEuler(), :dts => 1e-3 * multipliers),
-          Dict(:alg => CNAB2(), :dts => 1e-4 * multipliers),
-          Dict(:alg => CNLF2(), :dts => 1e-4 * multipliers),
-          Dict(:alg => SBDF2(), :dts => 1e-3 * multipliers)]
-labels = ["IMEXEuler" "CNAB2" "CNLF2" "SBDF2"]
-@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
-                            print_names=true, names=labels,
-                            numruns=5, error_estimate=:l2,
-                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
-
-plot(wp, label=labels, markershape=:auto, title="IMEX methods, dense linsolve, low order")
+The Allen-Cahn PDE is posed as follows:
+```math
+\begin{align}
+   \partial_t u = \gamma \nabla^2 u - u^3 + u \\
+\end{align}
 ```
-
-1.IMEX methods (Krylov linear solver)
-
-```julia
-abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
-reltols = 0.1 .^ (1:4)
-multipliers = 0.5 .^ (0:3)
-setups = [Dict(:alg => IMEXEuler(linsolve=KrylovJL_GMRES()), :dts => 1e-3 * multipliers),
-          Dict(:alg => CNAB2(linsolve=KrylovJL_GMRES()), :dts => 1e-4 * multipliers),
-          Dict(:alg => CNLF2(linsolve=KrylovJL_GMRES()), :dts => 1e-4 * multipliers),
-          Dict(:alg => SBDF2(linsolve=KrylovJL_GMRES()), :dts => 1e-3 * multipliers)]
-labels = ["IMEXEuler" "CNAB2" "CNLF2" "SBDF2"]
-@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
-                            print_names=true, names=labels,
-                            numruns=5, error_estimate=:l2,
-                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
-
-plot(wp, label=labels, markershape=:auto, title="IMEX methods, Krylov linsolve, low order")
+This is discretized in space via the finite difference method:
+```math
+\begin{align}
+    \partial_t u = \gamma D_x^2 u - u^3 + u
+\end{align}
 ```
 
-2. ExpRK methods
+## Packages
 
 ```julia
-abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
-reltols = 0.1 .^ (1:4)
-multipliers = 0.5 .^ (0:3)
-setups = [Dict(:alg => NorsettEuler(), :dts => 1e-3 * multipliers),
-          Dict(:alg => NorsettEuler(krylov=true, m=5), :dts => 1e-3 * multipliers),
-          Dict(:alg => NorsettEuler(krylov=true, m=20), :dts => 1e-3 * multipliers),
-          Dict(:alg => ETDRK2(), :dts => 1e-3 * multipliers),
-          Dict(:alg => ETDRK2(krylov=true, m=20), :dts => 1e-2 * multipliers),
-          Dict(:alg => ETDRK2(krylov=true, m=20), :dts => 1e-2 * multipliers)]
-labels = hcat("NorsettEuler (caching)", "NorsettEuler (m=5)", "NorsettEuler (m=20)",
-              "ETDRK2 (caching)", "ETDRK2 (m=5)", "ETDRK2 (m=20)")
-@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
-                            print_names=true, names=labels,
-                            numruns=5, error_estimate=:l2,
-                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
-
-plot(wp, label=labels, markershape=:auto, title="ExpRK methods, low order")
+using OrdinaryDiffEq
+using DiffEqDevTools
+using SciMLOperators
+using LinearSolve
+using LinearAlgebra
+using SparseArrays
+using Plots
+gr();
 ```
 
-## Between family comparisons
+## Implementation
 
 ```julia
-abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
-reltols = 0.1 .^ (1:4)
-multipliers = 0.5 .^ (0:3)
-setups = [Dict(:alg => CNAB2(), :dts => 1e-4 * multipliers),
-          Dict(:alg => CNAB2(linsolve=KrylovJL_GMRES()), :dts => 1e-4 * multipliers),
-          Dict(:alg => ETDRK2(), :dts => 1e-3 * multipliers)]
-labels = ["CNAB2 (dense linsolve)" "CNAB2 (Krylov linsolve)" "ETDRK2 (m=5)"]
-@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
-                            print_names=true, names=labels,
-                            numruns=5, error_estimate=:l2,
-                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
-
-plot(wp, label=labels, markershape=:auto, title="Between family, low orders")
-```
+# Linear term: Second derivative
+function linear_term(N, dx, nu)
+    # Second derivative operator (-u_xx)
+    e = ones(N + 1)
+    D2 = spdiagm(1 => e[1:end-1], 0 => -2e, -1 => e[1:end-1])
+    # Periodic boundary conditions
+    # D2[1, end] = 1.0 # Left ghost cell
+    # D2[end, 1] = 1.0 # Right ghost cell
+    D2 /= dx^2 # Central difference
+
+    return nu * D2
+end
 
-## Low tolerances
+# Nonlinear term: u - u^3
+function nonlinear_term(alpha)
+    function (du, u, p, t)
+        du .= alpha * (u .- u .^ 3)
+    end
+end
 
-## In-family comparisons
+function apply_boundary_conditions!(u, u_start, u_end)
+    u[1] = u_start 
+    u[end] = u_end
 
-1.IMEX methods (dense linear solver)
+    return u
+end
 
-```julia
-abstols = 0.1 .^ (7:13)
-reltols = 0.1 .^ (4:10)
-setups = [Dict(:alg => KenCarp3()),
-          Dict(:alg => KenCarp4()),
-          Dict(:alg => KenCarp5()),
-          Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:Dense)),
-          Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:Dense)),
-          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:Dense))]
-labels = hcat("KenCarp3", "KenCarp4", "KenCarp5", "ARKODE3", "ARKODE4", "ARKODE5")
-@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
-                            print_names=true, names=labels,
-                            numruns=5, error_estimate=:l2,
-                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
-
-plot(wp, label=labels, markershape=:auto, title="IMEX methods, dense linsolve, medium order")
+# Construct the problem
+function allen_cahn(N, L)
+    # Discretization
+    x = LinRange(-L, L, N + 1)
+    dx = x[2] - x[1]
+
+    # Linear and nonlinear terms
+    nu = 4.0 # Diffusion coefficient
+    alpha = 4.0 # Nonlinear coefficient
+    f1 = linear_term(N, dx, nu)
+    f2 = nonlinear_term(alpha)
+
+    u0 = @. cos(2π * x) # Initial condition
+
+    # Apply boundary conditions
+    fb1(du, u, p, t) = MatrixOperator(f1)(du, apply_boundary_conditions!(u, u[1], u[end]), p, t)
+    fb2(du, u, p, t) = f2(du, apply_boundary_conditions!(u, u[1], u[end]), p, t)
+
+    prob = SplitODEProblem(fb1, fb2, u0, (0.0, 1.0))
+
+    x, prob
+end;
 ```
 
-1.IMEX methods (Krylov linear solver)
+## Reference Solution
 
 ```julia
-abstols = 0.1 .^ (7:13)
-reltols = 0.1 .^ (4:10)
-setups = [Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES())),
-          Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES())),
-          Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES())),
-          Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:GMRES)),
-          Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:GMRES)),
-          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES))]
-labels = ["KenCarp3" "KenCarp4" "KenCarp5" "ARKODE3" "ARKODE4" "ARKODE5"]
-@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
-                            print_names=true, names=labels,
-                            numruns=5, error_estimate=:l2,
-                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
-
-plot(wp, label=labels, markershape=:auto, title="IMEX methods, medium order")
+N = 128 # Number of grid points
+L = 2.0  # Domain length
+xs, prob = allen_cahn(N, L)
+# @time sol = solve(prob, AutoVern7(RadauIIA5(autodiff=false)); dt=1e-2, abstol=1e-14, reltol=1e-14, adaptive=true)
+@time sol = solve(prob, Tsit5(); abstol=1e-14, reltol=1e-14, adaptive=true)
+
+tslices = LinRange(prob.tspan..., 50)
+ys = mapreduce(sol, hcat, tslices)
+plt = heatmap(xs, tslices, ys', xlabel="x", ylabel="t")
 ```
 
-2.ExpRK methods
+## Work-Precision Diagrams
 
 ```julia
-abstols = 0.1 .^ (7:11) # all fixed dt methods so these don't matter much
-reltols = 0.1 .^ (4:8)
-multipliers = 0.5 .^ (0:4)
-setups = [Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
-          Dict(:alg => ETDRK3(krylov=true, m=5), :dts => 1e-2 * multipliers),
-          Dict(:alg => ETDRK4(), :dts => 1e-2 * multipliers),
-          Dict(:alg => ETDRK4(krylov=true, m=5), :dts => 1e-2 * multipliers),
-          Dict(:alg => HochOst4(), :dts => 1e-2 * multipliers),
-          Dict(:alg => HochOst4(krylov=true, m=5), :dts => 1e-2 * multipliers)]
-labels = hcat("ETDRK3 (caching)", "ETDRK3 (m=5)", "ETDRK4 (caching)",
-              "ETDRK4 (m=5)", "HochOst4 (caching)", "HochOst4 (m=5)")
-@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
-                            print_names=true, names=labels,
-                            numruns=5, error_estimate=:l2,
-                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
-
-plot(wp, label=labels, markershape=:auto, title="ExpRK methods, medium order")
-```
-
-## Between family comparisons
+test_sol = TestSolution(sol); 
 
-```julia
-abstols = 0.1 .^ (7:11)
-reltols = 0.1 .^ (4:8)
-multipliers = 0.5 .^ (0:4)
-setups = [Dict(:alg => KenCarp5()),
-          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:Dense)),
-          Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES())),
-          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES)),
-          Dict(:alg => ETDRK3(krylov=true, m=5), :dts => 1e-2 * multipliers),
-          Dict(:alg => ETDRK4(krylov=true, m=5), :dts => 1e-2 * multipliers)]
-labels = hcat("KenCarp5 (dense linsolve)", "ARKODE (dense linsolve)", "KenCarp5 (Krylov linsolve)",
-              "ARKODE (Krylov linsolve)", "ETDRK3 (m=5)", "ETDRK4 (m=5)")
-@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
-                            print_names=true, names=labels,
-                            numruns=5, error_estimate=:l2,
-                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
-
-plot(wp, label=labels, markershape=:auto, title="Between family, medium order")
+abstols = 0.1 .^ (5:7) # all fixed dt methods so these don't matter much
+reltols = 0.1 .^ (1:3)
+multipliers = 0.5 .^ (0:3)
+setups = [
+    Dict(:alg => IMEXEuler(), :dts => 1e-4 * multipliers),
+    Dict(:alg => CNAB2(), :dts => 1e-4 * multipliers),
+    # Dict(:alg => CNLF2(), :dts => 1e-2 * multipliers), # Convergence fails
+    Dict(:alg => SBDF2(), :dts => 1e-4 * multipliers)
+]
+labels = hcat(
+    "IMEXEuler",
+    "CNAB2",
+    # "CNLF2",
+    "SBDF2",
+
+)
+@time wp = WorkPrecisionSet(prob, abstols, reltols, setups;
+    print_names=true, names=labels, numruns=5, error_estimate=:l2,
+    save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
+
+plot(wp, label=labels, markershape=:auto, title="Work-Precision Diagram, High Tolerance")
 ```
 
 ```julia, echo = false