[petsc-users] Direct solver for 3D-FEM

[Thomas Witkowski]

Jed Brown wrote:
> On Fri, May 20, 2011 at 14:55, Thomas Witkowski 
> <thomas.witkowski at tu-dresden.de 
> <mailto:thomas.witkowski at tu-dresden.de>> wrote:
>
>     Could one of you explain me why direct solvers (I make use of
>     UMFPACK) seems to work quite bad for 3D-FEM? For a small test
>     case, I take a unite square/box and refine it globally
>     (bisectioning of triangle/tetrahedrons). I solve a six order PDE
>      (that leads to symmetric but indefinite matrices) on it. In 2D
>     the resulting system has 49.923 rows with 808.199 non zeros,
>     UMFPACK solves the system within 2.8 seconds. For 3D I've refine
>     the box such that the resulting matrix has more or less the same
>     number of non zeros, in this case 898.807 on 27.027 rows. UMFPACK
>     needs 32 seconds to solve it, so more than 10 time as for the 2D
>     case. Even if I make use of fourth order Lagrange functions for
>     the 2D case, which leads to much denser matrices (49.923 rows and
>     2.705.927 non zeros), UMFPACK solves it within 3.2 seconds. Is
>     there any good reason for it?
>
>
> It has to do with maximal independent sets, or (more roughly) that the 
> ratio of surface area/volume (3D) is larger than perimeter/area (2D). 
> This also affects the amount of communication. You can try using MUMPS 
> instead of Umfpack, but it won't change the asymptotics.
>
> For PDE problems in 2D, direct solvers with optimal ordering need O(n 
> log n) memory and O(n^{1.5}) flops. In 3D, they need O(n^{4/3}) memory 
> and O(n^2) flops.
Thanks for the information. One related question: The system that is 
solved here is a 3x3 block matrix: Two of these blocks are empty, two 
are the discretized (FEM) identity, three are the discretized Laplace 
operator and one block is the discretization of a complex sum of 
operators. At all, there are only three different blocks. Can I make 
somehow advantages of this structure? I could create the three inverses 
of the three different blocks, but exists than a cheap/simple way to 
produce the action of the inverse of the whole 3x3 block matrix?

Thomas