Trying to solve KU=F, FEA problem, with roughly 5-10million DOFs

General Usage
, ,
6

Hi,
here is some testing code. See https://docs.sciml.ai/SciMLBenchmarksOutput/stable/LinearSolve/SparsePDE/ for the way how the matrix is created.
AMGCLWrap.jl wraps AMGCL implemented in C++ and could take advantage of OpenMP multithreading.

(Disclaimer: I am the maintainer of AMGCLWrap.jl)

module PCGTest
using LinearAlgebra, SparseArrays, LinearSolve
using ILUZero: ilu0
using AlgebraicMultigrid: smoothed_aggregation, aspreconditioner
using AMGCLWrap: CGSolver, AMGPrecon, AMGSolverAlgorithm

A ⊕ B = kron(I(size(B, 1)), A) + kron(B, I(size(A, 1)))

function lattice(n; Tv = Float64)
    d = fill(2 * one(Tv), n)
    d[1] = one(Tv)
    d[end] = one(Tv)
    spdiagm(1 => -ones(Tv, n - 1), 0 => d, -1 => -ones(Tv, n - 1))
end

lattice(L...; Tv = Float64) = lattice(L[1]; Tv) ⊕ lattice(L[2:end]...; Tv)

#
# Create a matrix similar to that of a finite difference discretization in a `dim`-dimensional
# unit cube of  ``-Δu + δu`` with approximately N unknowns. It is strictly diagonally dominant,
# symmetric, and thereforr positive definite
#
function fdmatrix(N; dim = 2, Tv = Float64, δ = 1.0e-2)
    n = N^(1 / dim) |> ceil |> Int
    lattice([n for i in 1:dim]...; Tv) + Tv(δ) * I
end


function test(N)
    @show N
    A=fdmatrix(N)
    u=rand(size(A,1))
    b=A*u
    pr=LinearProblem(A,b)
    @time "                              Plain CG" sol=solve(pr,KrylovJL_CG())
    @show norm(sol.u-u,Inf)
    @time "                           ILUZero.ilu0" sol=solve(pr,KrylovJL_CG(),Pl=ilu0(A))
    @show norm(sol.u-u,Inf)
    @time "AlgebraicMultigrid.smoothed_aggregation" sol=solve(pr,KrylovJL_CG(),Pl=aspreconditioner(smoothed_aggregation(A)))
    @show norm(sol.u-u,Inf)
    @time "                  AMGCLWrap.AMGPrecon()" sol=solve(pr,KrylovJL_CG(),Pl=AMGPrecon(A))
    @show norm(sol.u-u,Inf)
    # below, we use the built-in CG frorm AMGCL.
    @time "         AMGCLWrap.AMGSolverAlgorithm()" sol=solve(pr,AMGSolverAlgorithm(;solver=CGSolver()))
    @show norm(sol.u-u,Inf)
end

function main(;N=10_000_000)
    println("precompile:")
    test(1000)
    println("      run:")
    test(N)
end
end

Note to self: PR to SciMLBenchmarks…

1 Like