0
$\begingroup$

I implemented a FEM solver in MATLAB for Poisson's equation in 3D, using hexahedron and sparse matrix for the Laplacian. I was using the backslash but now I have to use a few iterative methods (GMRES with ilu for example). With $N_x$, $N_y$ and $N_z$ (in general, $N$) being the number of elements (50, 100 and 200 as test cases) I have some problems. For $N=50$ it takes more than $600$ seconds using cgs with incomplete LU preconditioner. For $N=100$ it takes more than 60 minutes and it is still running. I did not set the tolerance or the maximum number of iterations since it does not seem to affect the timing problem. Is this possible or it should not take that long? Is there a way to speed it up?

Here's the code:

[L,U] = ilu(Kb);
[u_h,flag,relres,iter,resvec] = cgs(Kb,fb,[],[],L,U);

where Kb is the Laplacian and fb the right-hand side.

Thank you.

Improve this question
$\endgroup$
3
  • $\begingroup$ Why do you think that tolerance and maximum number of iterations don't affect the timing? What is the time for 1 iteration? Do you have a consistent righthand-side (orthogonal to the kernel if the problem is singular)? Is the solution good for N=50? $\endgroup$
    – VorKir
    Dec 9, 2016 at 20:21
  • $\begingroup$ Is the time spent in computing the ILU, or in running the CG iteration? $\endgroup$
    – Wolfgang Bangerth
    Dec 9, 2016 at 22:18
  • $\begingroup$ @WolfgangBangerth Thank you for the reply. ichol is fast (less than one second for $N=50$), while LU decomposition is the most time consuming task (more than $25$ seconds for $N=50$, while to actually solve it with, for example, cgs it takes only $0.2$ seconds). $\endgroup$
    – wrong_path
    Dec 10, 2016 at 6:53

1 Answer 1

Reset to default
2
$\begingroup$

Your stiffness matrix is symmetric, positive definite so the Krylov-based solver of choice for this problem is conjugate gradients. You want code similar to this

L = ichol(K);
u = pcg(K, b, tol, maxIter, L, L');

assuming you want to use an incomplete factorization as a preconditioner. The stiffness matrix for Poisson's equation on a cube is relatively well-conditioned so you might want to experiment with no preconditioner or a very simple preconditioner like diagonal scaling. But, since you say the MATLAB implementation of ichol is quite fast, that is likely to be best choice.

Other Krylov-based iterative algorithms, like cgs and GMRES that you mention, are primarily for unsymmetric or non-positive definite systems.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.