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 <[email protected]> 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 <[email protected]> > 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 <[email protected]> > 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 >> [email protected], 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 <[email protected]> >> 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 >>> [email protected], I can have a look. >>> >>> Thanks, >>> Pierre >>> >>> On 6 Sep 2023, at 8:14 AM, Carl-Johan Thore <[email protected]> >>> 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 <[email protected]> >>> 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 <[email protected]> >>>> wrote: >>>> >>>> >>>> On 28 Aug 2023, at 6:50 PM, Carl-Johan Thore <[email protected]> >>>> 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 < >>>> [email protected]> wrote: >>>> >>>>> >>>>> >>>>> On 27 Aug 2023, at 12:14 AM, Carl-Johan Thore <[email protected]> >>>>> 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 <[email protected]> >>>>> *Sent:* Saturday, August 26, 2023 4:36 PM >>>>> *To:* Carl-Johan Thore <[email protected]> >>>>> *Cc:* Carl-Johan Thore <[email protected]>; >>>>> [email protected] >>>>> *Subject:* Re: [petsc-users] PCFIELDSPLIT with MATSBAIJ >>>>> >>>>> >>>>> >>>>> >>>>> On 26 Aug 2023, at 11:16 PM, Carl-Johan Thore <[email protected]> >>>>> 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 < >>>>> [email protected]> 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> >>> >>> >
