[petsc-users] PCFIELDSPLIT with MATSBAIJ

[Carl-Johan Thore]

Great, that's as I suspected then.
The matrix comes from density-based topology optimization using an FE-model
for Stokes-Brinkman with pressure jump stabilization,
so it's a difficult problem I would say.
I'll experiment with the inner solvers and see what works.
Thanks again for the help,
Carl-Johan

On Wed, Sep 6, 2023 at 2:16 PM Pierre Jolivet <pierre.jolivet at lip6.fr>
wrote:

> (switching back to petsc-users with the previous attachment striped)
> I’m assuming your problem is kind of ill-posed?
> PCICC and PCILU have different default parameters, and thus can yield
> different convergence curves, especially for difficult problems.
> When switching from AIJ to SBAIJ, ILU will switch to ICC.
> If I switch to something more “robust” (PCCHOLESKY from MUMPS), then I get
> the same convergence (up to machine precision).
> I think this highlight there is no issue with the new
> PCFIELDSPLIT/MatCreateSubMatri[x,ces]() code, but merely that you need to
> be a bit careful with the inner solvers.
>
> Thanks,
> Pierre
>
>
>
> On 6 Sep 2023, at 11:27 AM, Carl-Johan Thore <carljohanthore at gmail.com>
> wrote:
>
> Ok, here's data from a small instance of my problem. It's in Matlab-format
> which is what I usually use, but let me know if you want something else.
> I can send smaller versions of the problem if you prefer that (really
> small versions will take a bit of time to code so that one doesn't just get
> the solution
> 0 everywhere).
>
> One IS, ISO, is included because in my code PETSc computes the other IS as
> the complement. Also included is the output from PETSc for aij and
> sbaij, and a Matlab-script to check the data.
>
> Kind regards,
> Carl-Johan
>
> On Wed, Sep 6, 2023 at 10:10 AM Carl-Johan Thore <carljohanthore at gmail.com>
> wrote:
>
>> "Naïve question, but is your problem really symmetric?
>> For example, if you do FEM, how do you impose Dirichlet boundary
>> conditions?"
>>
>> The structure of the matrix is K =  [A B; B^T C] with A and C symmetric,
>> so the problem should be symmetric. Numerically, the maximum relative
>> difference between K and K^T is
>> on the order of, let's say, 1e-18 which I guess is a good as can be
>> expected? I'm using FEM with (homogeneous) Dirichlet boundary conditions
>> imposed by zeroing rows and columns
>> for the fixed DOFs and adding ones to the corresponding diagonal
>> elements, i.e. same as what MatZeroRowsColumns does.  MatZeroRowsColumns
>> doesn't work for sbaij but we already have an implementation in which rows
>> and columns are cancelled in the element matrices before assembly which is
>> well-tested for the aij-case. We've also verified numerically that aij and
>> sbaij yields the same upper triangular part, and that the RHS is the same
>> in both cases.
>>
>> "If you can reproduce this on a smaller example and are at liberty to
>> share the Mat, IS, and RHS Vec, feel free to send them at
>> petsc-maint at mcs.anl.gov, I can have a look."
>>
>> Yes, this is reproducible on smaller problems. In case I send you Mat,
>> IS, and RHS, which format is preferable?
>>
>> Kind regards,
>> Carl-Johan
>>
>>
>> On Wed, Sep 6, 2023 at 8:21 AM Pierre Jolivet <pierre.jolivet at lip6.fr>
>> wrote:
>>
>>> Naïve question, but is your problem really symmetric?
>>> For example, if you do FEM, how do you impose Dirichlet boundary
>>> conditions?
>>> If you can reproduce this on a smaller example and are at liberty to
>>> share the Mat, IS, and RHS Vec, feel free to send them at
>>> petsc-maint at mcs.anl.gov, I can have a look.
>>>
>>> Thanks,
>>> Pierre
>>>
>>> On 6 Sep 2023, at 8:14 AM, Carl-Johan Thore <carljohanthore at gmail.com>
>>> wrote:
>>>
>>> 
>>> Thanks for this!
>>> I now get convergence with sbaij and fieldsplit. However, as can be seen
>>> in the attached file, it requires around three times as many outer
>>> iterations as with aij to reach the same point (to within rtol). Switching
>>> from -pc_fieldsplit_schur_fact_type LOWER to FULL (which is the "correct"
>>> choice?) helps quite a bit, but still sbaij takes many more iterations. So
>>> it seems that your sbaij-implementation works but that I'm doing something
>>> wrong when setting up the fieldsplit pc, or that some of the subsolvers
>>> doesn't work properly with sbaij. I've tried bjacobi as you had in your
>>> logs, but not hypre yet. But anyway, one doesn't expect the matrix format
>>> to have this much impact on convergence? Any other suggestions?
>>> /Carl-Johan
>>>
>>> On Mon, Sep 4, 2023 at 9:31 PM Pierre Jolivet <Pierre.Jolivet at lip6.fr>
>>> wrote:
>>>
>>>> The branch should now be good to go (
>>>> https://gitlab.com/petsc/petsc/-/merge_requests/6841).
>>>> Sorry, I made a mistake before, hence the error on PetscObjectQuery().
>>>> I’m not sure the code will be covered by the pipeline, but I have
>>>> tested this on a Raviart—Thomas discretization with PCFIELDSPLIT.
>>>> You’ll see in the attached logs that:
>>>> 1) the numerics match
>>>> 2) in the SBAIJ case, PCFIELDSPLIT extract the (non-symmetric) A_{01}
>>>> block from the global (symmetric) A and we get the A_{10} block cheaply by
>>>> just using MatCreateHermitianTranspose(), instead of calling another time
>>>> MatCreateSubMatrix()
>>>> Please let me know if you have some time to test the branch and whether
>>>> it fails or succeeds on your test cases.
>>>>
>>>> Also, I do not agree with what Hong said.
>>>> Sometimes, the assembly of a coefficient can be more expensive than the
>>>> communication of the said coefficient.
>>>> So they are instances where SBAIJ would be more efficient than AIJ even
>>>> if it would require more communication, it is not a black and white picture.
>>>>
>>>> Thanks,
>>>> Pierre
>>>>
>>>>
>>>> On 28 Aug 2023, at 12:12 PM, Pierre Jolivet <Pierre.Jolivet at lip6.fr>
>>>> wrote:
>>>>
>>>>
>>>> On 28 Aug 2023, at 6:50 PM, Carl-Johan Thore <carljohanthore at gmail.com>
>>>> wrote:
>>>>
>>>> I've tried the new files, and with them, PCFIELDSPLIT now gets set up
>>>> without crashes (but the setup is significantly slower than for MATAIJ)
>>>>
>>>>
>>>> I’ll be back from Japan at the end of this week, my schedule is too
>>>> packed to get anything done in the meantime.
>>>> But I’ll let you know when things are working properly (last I checked,
>>>> I think it was working, but I may have forgotten about a corner case or
>>>> two).
>>>> But yes, though one would except things to be faster and less memory
>>>> intensive with SBAIJ, it’s unfortunately not always the case.
>>>>
>>>> Thanks,
>>>> Pierre
>>>>
>>>> Unfortunately I still get errors later in the process:
>>>>
>>>> [0]PETSC ERROR: Null argument, when expecting valid pointer
>>>> [0]PETSC ERROR: Null Pointer: Parameter # 1
>>>> [0]PETSC ERROR: Petsc Development GIT revision:
>>>> v3.19.4-1023-ga6d78fcba1d  GIT Date: 2023-08-22 20:32:33 -0400
>>>> [0]PETSC ERROR: Configure options -f --with-fortran-bindings=0
>>>> --with-cuda --with-cusp --download-scalapack --download-hdf5
>>>> --download-zlib --download-mumps --download-parmetis --download-metis
>>>> --download-ptscotch --download-hypre --download-spai
>>>> [0]PETSC ERROR: #1 PetscObjectQuery() at
>>>> /mnt/c/mathware/petsc/petsc-v3-19-4/src/sys/objects/inherit.c:742
>>>> [0]PETSC ERROR: #2 MatCreateSubMatrix_MPISBAIJ() at
>>>> /mnt/c/mathware/petsc/petsc-v3-19-4/src/mat/impls/sbaij/mpi/mpisbaij.c:1414
>>>> [0]PETSC ERROR: #3 MatCreateSubMatrix() at
>>>> /mnt/c/mathware/petsc/petsc-v3-19-4/src/mat/interface/matrix.c:8476
>>>> [0]PETSC ERROR: #4 PCSetUp_FieldSplit() at
>>>> /mnt/c/mathware/petsc/petsc-v3-19-4/src/ksp/pc/impls/fieldsplit/fieldsplit.c:826
>>>> [0]PETSC ERROR: #5 PCSetUp() at
>>>> /mnt/c/mathware/petsc/petsc-v3-19-4/src/ksp/pc/interface/precon.c:1069
>>>> [0]PETSC ERROR: #6 KSPSetUp() at
>>>> /mnt/c/mathware/petsc/petsc-v3-19-4/src/ksp/ksp/interface/itfunc.c:415
>>>>
>>>> The code I'm running here works without any problems for MATAIJ. To run
>>>> it with MATSBAIJ I've simply used the command-line option
>>>> -dm_mat_type sbaij
>>>>
>>>>
>>>> Kind regards,
>>>> Carl-Johan
>>>>
>>>>
>>>> On Sat, Aug 26, 2023 at 5:21 PM Pierre Jolivet via petsc-users <
>>>> petsc-users at mcs.anl.gov> wrote:
>>>>
>>>>>
>>>>>
>>>>> On 27 Aug 2023, at 12:14 AM, Carl-Johan Thore <carl-johan.thore at liu.se>
>>>>> wrote:
>>>>>
>>>>> “Well, your A00 and A11 will possibly be SBAIJ also, so you’ll end up
>>>>> with the same issue.”
>>>>> I’m not sure I follow. Does PCFIELDSPLIT extract further submatrices
>>>>> from these blocks, or is there
>>>>> somewhere else in the code that things will go wrong?
>>>>>
>>>>>
>>>>> Ah, no, you are right, in that case it should work.
>>>>>
>>>>> For the MATNEST I was thinking to get some savings from the
>>>>> block-symmetry at least
>>>>> even if symmetry in A00 and A11 cannot be exploited; using SBAIJ for
>>>>> them would just be a
>>>>> (pretty big) bonus.
>>>>>
>>>>> “I’ll rebase on top of main and try to get it integrated if it could
>>>>> be useful to you (but I’m traveling
>>>>> right now so everything gets done more slowly, sorry).”
>>>>> Sound great, Thanks again!
>>>>>
>>>>>
>>>>> The MR is there https://gitlab.com/petsc/petsc/-/merge_requests/6841.
>>>>> I need to add a new code path in MatCreateRedundantMatrix() to make
>>>>> sure the resulting Mat is indeed SBAIJ, but that is orthogonal to the
>>>>> PCFIELDSPLIT issue.
>>>>> The branch should be usable in its current state.
>>>>>
>>>>> Thanks,
>>>>> Pierre
>>>>>
>>>>>
>>>>> *From:* Pierre Jolivet <pierre.jolivet at lip6.fr>
>>>>> *Sent:* Saturday, August 26, 2023 4:36 PM
>>>>> *To:* Carl-Johan Thore <carl-johan.thore at liu.se>
>>>>> *Cc:* Carl-Johan Thore <carljohanthore at gmail.com>;
>>>>> petsc-users at mcs.anl.gov
>>>>> *Subject:* Re: [petsc-users] PCFIELDSPLIT with MATSBAIJ
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On 26 Aug 2023, at 11:16 PM, Carl-Johan Thore <carl-johan.thore at liu.se>
>>>>> wrote:
>>>>>
>>>>> "(Sadly) MATSBAIJ is extremely broken, in particular, it cannot be
>>>>> used to retrieve rectangular blocks in MatCreateSubMatrices, thus you
>>>>> cannot get the A01 and A10 blocks in PCFIELDSPLIT.
>>>>> I have a branch that fixes this, but I haven’t rebased in a while (and
>>>>> I’m AFK right now), would you want me to rebase and give it a go, or must
>>>>> you stick to a release tarball?"
>>>>>
>>>>> Ok, would be great if you could look at this! I don't need to stick to
>>>>> any particular branch.
>>>>>
>>>>> Do you think MATNEST could be an alternative here?
>>>>>
>>>>>
>>>>> Well, your A00 and A11 will possibly be SBAIJ also, so you’ll end up
>>>>> with the same issue.
>>>>> I’m using both approaches (monolithic SBAIJ or Nest + SBAIJ), it was
>>>>> crashing but I think it was thoroughly fixed in
>>>>> https://gitlab.com/petsc/petsc/-/commits/jolivet/feature-matcreatesubmatrices-rectangular-sbaij/
>>>>> It is ugly code on top of ugly code, so I didn’t try to get it
>>>>> integrated and just used the branch locally, and then moved to some other
>>>>> stuff.
>>>>> I’ll rebase on top of main and try to get it integrated if it could be
>>>>> useful to you (but I’m traveling right now so everything gets done more
>>>>> slowly, sorry).
>>>>>
>>>>> Thanks,
>>>>> Pierre
>>>>>
>>>>>
>>>>> My matrix is
>>>>> [A00 A01;
>>>>> A01^t A11]
>>>>> so perhaps with MATNEST I can make use of the block-symmetry at least,
>>>>> and then use MATSBAIJ for
>>>>> A00 and A11 if it's possible to combine matrix types which the manual
>>>>> seems to imply.
>>>>>
>>>>> Kind regards
>>>>> Carl-Johan
>>>>>
>>>>>
>>>>>
>>>>> On 26 Aug 2023, at 10:09 PM, Carl-Johan Thore <
>>>>> carljohanthore at gmail.com> wrote:
>>>>> 
>>>>> Hi,
>>>>>
>>>>> I'm trying to use PCFIELDSPLIT with MATSBAIJ in PETSc 3.19.4.
>>>>> According to the manual "[t]he fieldsplit preconditioner cannot
>>>>> currently be used with the MATBAIJ or MATSBAIJ data formats if the
>>>>> blocksize is larger than 1". Since my blocksize is exactly 1 it would
>>>>> seem that I can use PCFIELDSPLIT. But this fails with "PETSC ERROR: For
>>>>> symmetric format, iscol must equal isrow"
>>>>> from MatCreateSubMatrix_MPISBAIJ. Tracing backwards one ends up in
>>>>> fieldsplit.c at
>>>>>
>>>>> /* extract the A01 and A10 matrices */ ilink = jac->head;
>>>>> PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis)); if
>>>>> (jac->offdiag_use_amat) { PetscCall(MatCreateSubMatrix(pc->mat,
>>>>> ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B)); } else {
>>>>>        PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis,
>>>>> MAT_INITIAL_MATRIX, &jac->B)); }
>>>>>
>>>>> This, since my A01 and A10 are not square, seems to explain why iscol
>>>>> is not equal to isrow.
>>>>> From this I gather that it is in general NOT possible to use
>>>>> PCFIELDSPLIT with MATSBAIJ even with block size 1?
>>>>>
>>>>> Kind regards,
>>>>> Carl-Johan
>>>>>
>>>>>
>>>>>
>>>>
>>>> <experiment_with_sbaij.txt>
>>>
>>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230906/d6084319/attachment-0001.html>