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

**URL:** <https://discourse.julialang.org/t/trying-to-solve-ku-f-fea-problem-with-roughly-5-10million-dofs/112480>\
**Category:** General Usage\
**Tags:** package, linearalgebra, linearsolve\
**Created:** [April 3, 2024, 5:47pm UTC](https://discourse.julialang.org/t/trying-to-solve-ku-f-fea-problem-with-roughly-5-10million-dofs/112480 "2024-04-03T17:47:08Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![balaji\_sriram](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balaji_sriram/32/208000_2.png) [@balaji\_sriram](https://discourse.julialang.org/u/balaji_sriram)\
**Post date:** [April 3, 2024, 5:47pm UTC](https://discourse.julialang.org/t/trying-to-solve-ku-f-fea-problem-with-roughly-5-10million-dofs/112480/1 "2024-04-03T17:47:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![j-fu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j-fu/32/11373_2.png) [@j-fu](https://discourse.julialang.org/u/j-fu)\
**Post date:** [April 3, 2024, 7:22pm UTC](https://discourse.julialang.org/t/trying-to-solve-ku-f-fea-problem-with-roughly-5-10million-dofs/112480/2 "2024-04-03T19:22:32Z")

</div>

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 .

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 4, 2024, 8:11am UTC](https://discourse.julialang.org/t/trying-to-solve-ku-f-fea-problem-with-roughly-5-10million-dofs/112480/3 "2024-04-04T08:11:37Z")

</div>

> [@balaji\_sriram](#):
>
> I have also tried the below code: (using cg) but it isalso horribly slow:
> 
> sol = IterativeSolvers.cg(K\_stiffness, Force)

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())`?

---

<div class="post-metadata">

**Author:** ![Balaji\_Lakshminaraya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balaji_lakshminaraya/32/205881_2.png) [@Balaji\_Lakshminaraya](https://discourse.julialang.org/u/Balaji_Lakshminaraya)\
**Post date:** [April 4, 2024, 8:14am UTC](https://discourse.julialang.org/t/trying-to-solve-ku-f-fea-problem-with-roughly-5-10million-dofs/112480/4 "2024-04-04T08:14:55Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Balaji\_Lakshminaraya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balaji_lakshminaraya/32/205881_2.png) [@Balaji\_Lakshminaraya](https://discourse.julialang.org/u/Balaji_Lakshminaraya)\
**Post date:** [April 4, 2024, 8:16am UTC](https://discourse.julialang.org/t/trying-to-solve-ku-f-fea-problem-with-roughly-5-10million-dofs/112480/5 "2024-04-04T08:16:38Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![j-fu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j-fu/32/11373_2.png) [@j-fu](https://discourse.julialang.org/u/j-fu)\
**Post date:** [April 4, 2024, 2:11pm UTC](https://discourse.julialang.org/t/trying-to-solve-ku-f-fea-problem-with-roughly-5-10million-dofs/112480/6 "2024-04-04T14:11:44Z")

</div>

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

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

```julia
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…

---

<div class="post-metadata">

**Author:** ![balaji\_sriram](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balaji_sriram/32/208000_2.png) [@balaji\_sriram](https://discourse.julialang.org/u/balaji_sriram)\
**Post date:** [April 5, 2024, 2:15pm UTC](https://discourse.julialang.org/t/trying-to-solve-ku-f-fea-problem-with-roughly-5-10million-dofs/112480/7 "2024-04-05T14:15:33Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![balaji\_sriram](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balaji_sriram/32/208000_2.png) [@balaji\_sriram](https://discourse.julialang.org/u/balaji_sriram)\
**Post date:** [April 5, 2024, 2:20pm UTC](https://discourse.julialang.org/t/trying-to-solve-ku-f-fea-problem-with-roughly-5-10million-dofs/112480/8 "2024-04-05T14:20:48Z")

</div>

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](mailto:Balajisriram15@gmail.com).