[petsc-users] Saddle point problem with nested matrix and a relatively small number of Lagrange multipliers

[Matthew Knepley]

On Thu, Sep 3, 2020 at 11:43 AM Olivier Jamond <olivier.jamond at cea.fr>
wrote:

> Hello,
>
> I am working on a finite-elements/finite-volumes code, whose distributed
> solver is based on petsc. For FE, it relies on Lagrange multipliers for
> the imposition of various boundary conditions or interactions (simple
> dirichlet, contact, ...). This results in saddle point problems:
>
> [K C^t][U]=[F]
> [C  0 ][L] [D]
>
> Most of the time, the relations related to the matrix C are applied to
> dofs on the boundary of the domain. Then the size of L is much smaller
> than the size of U, which becomes more and more true as  the mesh is
> refined.
>
> The code construct this matrix as a nested matrix (one of the reason is
> that for some interactions such as contact, whereas being quite small,
> the size of the matrix C change constantly, and having it 'blended' into
> a monolithic 'big' matrix would require to recompute its profile/
> reallocate / ... each time), and use fieldsplit preconditioner of type
> PC_COMPOSITE_SCHUR. I would like to solve the system using iterative
> methods to access good extensibility on a large number of subdomains.
>
> Simple BC such as Dirichlet can be eliminated into K (and must be in
> order to make K invertible).
>
> My problem is the computational cost of these constraints treated with
> Lagrange multipliers, whereas their number becomes more and more
> neglectable as the mesh is refined. To give an idea, let's consider a
> simple elastic cube with dirichlet BCs which are all eliminated (to
> ensure invertibility of K) but one on a single dof.
>
> -ksp_type preonly
> -pc_type fieldsplit
> -pc_fieldsplit_type schur
> -pc_fieldsplit_schur_factorization_type full
> -pc_fieldsplit_schur_precondition selfp
>
> -fieldsplit_u_ksp_type cg
> -fieldsplit_u_pc_type bjacobi
>
> -fieldsplit_l_ksp_type cg
> -fieldsplit_l_pc_type bjacobi
>
> it seems that my computation time is multiplied by a factor 3: 3 ksp
> solves of the big block 'K' are needed to apply the schur preconditioner
> (assuming that the ksp(S,Sp) converges in 1 iteration). It seems
> expensive for a single dof dirichlet!
>

I am not sure you can get around this cost. In this case, it reduces to the
well-known
Sherman-Morrison formula (
https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula),
which Woodbury generalized. It seems to have the same number of solves.

  Thanks,

     Matt


> And for some relations treated by Lagrange multipliers which involve
> many dofs, the number of ksp solve of the big block 'K' is ( 2 + number
> of iteration of ksp(S,Sp)). To reduce this, one can think about solving
> the ksp(S,Sp) with a direct solver, but then one must use
> "-pc_fieldsplit_schur_precondition self" which is advised against in the
> documentation...
>
> To illustrate this, on a small elasticity case: 32x32x32 cube on 8
> processors, dirichlet on the top and bottom faces:
> * if all the dirichlet are eliminated (no C matrix, monolithic solve of
> the K bloc)
>    - computation time for the solve: ~400ms
> * if only the dirichlet of the bottom face are eliminated
>    - computation time for the solve: ~35000ms
>    - number of iteration of ksp(S,Sp): 37
>    - total number of iterations of ksp(K): 4939
> * only the dirichlet of the bottom face are eliminated with these options:
> -ksp_type fgmres
> -pc_type fieldsplit
> -pc_fieldsplit_type schur
> -pc_fieldsplit_schur_factorization_type full
> -pc_fieldsplit_schur_precondition selfp
>
> -fieldsplit_u_ksp_type cg
> -fieldsplit_u_pc_type bjacobi
>
> -fieldsplit_l_ksp_type cg
> -fieldsplit_l_pc_type bjacobi
> -fieldsplit_l_ksp_rtol 1e-10
> -fieldsplit_l_inner_ksp_type preonly
> -fieldsplit_l_inner_pc_type jacobi
> -fieldsplit_l_upper_ksp_type preonly
> -fieldsplit_l_upper_pc_type jacobi
>
>    - computation time for the solve: ~50000ms
>    - total number of iterations of ksp(K): 7424
>    - 'main' ksp number of iterations: 7424
>
> Then in the end, my question is: is there a smarter way to handle such
> 'small' constraint matrices C, with the (maybe wrong) idea that a small
> number of extra dofs (the lagrange multipliers) should result in a small
> extra computation time ?
>
> Thanks!
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200907/aee474ba/attachment.html>