> On Oct 6, 2020, at 6:57 AM, Olivier Jamond <[email protected]> wrote:
>
>
> On 03/10/2020 00:23, Barry Smith wrote:
>> I think what Jed is saying is that you should just actually build your
>> preconditioner for your Ct*C + Qt*S*Q operator with S. Because Ct is tall
>> and skinny the eigenstructure of Ct*C + Qt*S*Q is just the eigenstructure of
>> S with a low rank "modification" and Krylov methods (GMRES) are good at
>> solving problems where the eigenstructure of the preconditioner is only a
>> small rank modification of the eigenstructure of the operator you are supply
>> to GMRES. In the best situation each new iteration of GMRES corrects one
>> more of the "rogue" eigen directions. I would first use a direct solver with
>> S just to test how well it works as a preconditioner and then switch to GAMG
>> or whatever should work efficiently for solving your particular S matrix.
>>
>> I'd be interested in hearing how well the Ainsworth Formula works, it is
>> something that might be worth adding to PCFIELDSPLIT.
>>
>>
>> Barry
>
> Hi Barry,
>
> Thanks for these clarifications.
>
> To give some context, the test I am working on is a traction on an
> elastoplastic cube in large strain on which I apply 2% of strain at the first
> loading increment. The cube has 14739 dofs, and the number of rows of the C
> matrix is 867.
>
> In this simple case, the C matrix just refers to simple dirichlet conditions.
> Then Q is diagonal with 1. on dofs without dirichlet on 0. for dofs with
> dirichlets. Q'*S*Q is like S with zeros on lines/columns referring to dofs
> with dirichlet, and then C'*C just re-add non null value on the diagonal for
> the dofs with dirichlet. In the end, I feel that in this case the ainsworth
> method just do exactly the same as row/column elimination that can be done
> with MatZeroRowsColumns and the x and b optional vectors provided.
>
> On this test, with '-ksp_rtol 1.e-9' and '-ksp_type gmres', using S as a
> preconditionning matrix and a direct solver gives 65 iterations of the gmres
> for my first newton iteration (where S is SPD) and between 170 and 290 for
> the next ones (S is still symmetric but has negative eigenvalues). If I use
> '-pc_type gamg', the number of iterations of the gmres for the first (SPD)
> newton iteration is (14 with Sp / 23 with S), and for the next ones (not SPD)
> it is (~45 with Sp / ~180 with S).
Given the structure of C it seems you should just explicitly construct Sp
and use GAMG (or other preconditioners, even a direct solver) directly on Sp.
Trying to avoid explicitly forming Sp will give you a much slower performing
solving for what benefit? If C was just some generic monster than forming Sp
might be unrealistic but in your case CCt is is block diagonal with tiny blocks
which means (C*Ct)^(-1) is block diagonal with tiny blocks (the blocks are the
inverses of the blocks of (C*Ct)).
Sp = Ct*C + Qt * S * Q = Ct*C + [I - Ct * (C*Ct)^(-1)*C] S [I - Ct *
(C*Ct)^(-1)*C]
[Ct * (C*Ct)^(-1)*C] will again be block diagonal with slightly larger blocks.
You can do D = (C*Ct) with MatMatMult() then write custom code that zips
through the diagonal blocks of D inverting all of them to get iD then use
MatPtAP applied to C and iD to get Ct * (C*Ct)^(-1)*C then MatShift() to
include the I then MatPtAP or MatRAR to get [I - Ct * (C*Ct)^(-1)*C] S [I - Ct
* (C*Ct)^(-1)*C] then finally MatAXPY() to get Sp. The complexity of each of
the Mat operations is very low because of the absurdly simple structure of C
and its descendants. You might even be able to just use MUMPS to give you the
explicit inv(C*Ct) without writing custom code to get iD.
Perhaps I am missing something? Note you can also prototype this process in
Matlab very quickly to find any glitches. Hopefully you will find [I - Ct *
(C*Ct)^(-1)*C] S [I - Ct * (C*Ct)^(-1)*C] has only a slightly "bulked out"
sparsity of S.
>
> In this case with only simple dirichlets, I think I would like that the
> PCApply does something like: (I-Q)*jacobi(Ct*C)*(I-Q) + Q*precond(S)*Q. BUt I
> am not sure how to do that (I am quite newbie with petsc)... With a PCShell
> and PCShellSetApply?
I don't think there is any reason to think that using I-Q)*jacobi(Ct*C)*(I-Q)
+ Q*precond(S)*Q. would be a particularly good preconditioner for Sp. That is
much better than using a preconditioner built from S. But you can use
PCCOMPOSITE and PCSHELL with KSP inside for the two jacobi(Ct*C) and precond(S).
Barry
>
> In the end, if we found something that works well with the ainsworth formula,
> it would be nice to have it natively with PCFIELDSPLIT!
>
> Many thanks,
> Olivier
>
