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

General Usage
, ,
1

Hi, Im new to Julia, but i have experience with MATLAB, im trying to solve a KU=F, FEA type problem , with over 1 million DOFs. Note: the K matrix is sparse, symmetric, positive definite and banded. In matlab, the command U = pcg(K, F ,tol ,maxit=1000 , M1, M2, x0), converges very quickly. However, in Julia im finding it very difficult to pick a PCG solver of similar performance. Im using the following code.

using MAT, LinearAlgebra, LinearSolve , SparseArrays
prob = LinearProblem(K_stiffness, Force);
sol = solve(prob,LinearSolve.SimpleGMRES());

I have also tried the below code: (using cg) but it isalso horribly slow:

sol = IterativeSolvers.cg(K_stiffness, Force)

Please suggest packages and corrections to solve this linar problem quickly. Also suggest packages that may use PCG() method specifically, because it has proven its competence (in matlab) for this specific problem.

1 Like
2

It seems, you are running without a preconditioner here. What is M1, M2 in the matlab example?

Can you say more about your problem ? Is it elasticity ? Diffusion ? 2D or 3D ?

You can run LinearSolve with KrylovJL_CG instead of GMRES. You can find preconditioners in ILUZero.jl and AlgebraicMultigrid.jl .

3 Likes
3

The LinearSolve documentation explicitly recommends against the IterativeSolvers.jl implementations and recommends the Krylov.jl ones for this reason. So did you try sol = solve(prob,LinearSolve.KrylovJL_CG())?

3 Likes
4

M1 and M2 are not fed in matlab. In the place of M1,M2 in the matlab pcg syntax I just give [ ], [ ] . Also my problem is an electrostatic problem. So it’s more like [K]*[Phi] =[ F ]. Where phi is nodal electric potential. . It’s a 3D problem with uniform, linear cuboidal mesh elements, hence a lot of DOFs. Thank you for the recommendation for gmres and krylov. Can you please give me the code for solving this problem. I’ll try out and let you know. Thanks in advance.

5

Oh chris, I love your videos. I’m a big fan. Thanks for your suggestion. I will try it out and let you know. How do I feed things like max iteration and tolerance?

Ps: Also your video on solving differential equations with Julia saved me a lot of time in my research. Thanks a lot. Ho

1 Like
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
7

Thanks Chris, I tried your method and the bunch of other methods given by j-Fu. It turns out, there was a scaling issue. Once I scaled my force vector with an appropriate constant, I tried all these methods and got good results.

8

Much appreciated!!!. Thanks. There was a scaling issue. In addition I tried all of your specified algorithms, Plain CG and AMGCLWrap.AMGPrecon were the fastest.
Also @ChrisRackauckas and @j-fu , I work in the field of electromagnetics, high voltage physics, numerical methods, wireless power transfer etc, if you need to know anything in these specific areas I would love to return the favor. Dropping my email ID for the same:
Balajis@iitk.ac.in, Balajisriram15@gmail.com.