Add advanced tutorials - Githubissues

ErikQQY commented 1 year ago

What is the correct way of adding a preconditioner for linear solvers in NonlinearSolve.jl? I wrote this tutorial based on solving large stiff ODE in diffeq docs. I think the preconditioners should be the same since we have a common LinearSolve.jl interface, but the preconditioners failed when I passed them to solve.

MWE:

using NonlinearSolve, LinearAlgebra, SparseArrays, LinearSolve

const N = 32
const xyd_brusselator = range(0, stop = 1, length = N)
brusselator_f(x, y) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * 5.0
limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a
function brusselator_2d_loop(du, u, p)
    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)
        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_brusselator_2d = NonlinearProblem(brusselator_2d_loop, u0, p)
using Symbolics
du0 = copy(u0)
jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p),
    du0, u0)
ff = NonlinearFunction(brusselator_2d_loop; jac_prototype = float.(jac_sparsity))
prob_brusselator_2d_sparse = NonlinearProblem(ff, u0, p)
using IncompleteLU
function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
    if newW === nothing || newW
        Pl = ilu(convert(AbstractMatrix, W), τ = 50.0)
    else
        Pl = Plprev
    end
    Pl, nothing
end
solve(prob_brusselator_2d_sparse,
    NewtonRaphson(linsolve = KrylovJL_GMRES(), precs = incompletelu,
        concrete_jac = true));

codecov[bot] commented 1 year ago

Codecov Report

Merging #213 (d1910b4) into master (5c89c59) will increase coverage by 1.08%. Report is 2 commits behind head on master. The diff coverage is n/a.

@@            Coverage Diff             @@
##           master     #213      +/-   ##
==========================================
+ Coverage   19.86%   20.95%   +1.08%     
==========================================
  Files           8        8              
  Lines         735      735              
==========================================
+ Hits          146      154       +8     
+ Misses        589      581       -8

see 2 files with indirect coverage changes

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more