MatNest with Shell blocks for multipysics
Hi. We are designing a PETSc application that will employ a SNES solver on a multiphysics problem whose jacobian will have a 2x2 block form, say A=[A00,A01;A10,A11]. We already have code for the top left block A_00 (a MatShell and a related Shell preconditioner) that we wish to reuse. We could implement the other blocks as Shells or assembled matrices. We'd like also to compare our method with existing ones, so we'd like to be quite flexible in the choice of KSP and PC within the SNES. (To this end, implementing an assembled version of the A00 and the other blocks would be easy) I am assuming that, in order to have one or more shell blocks, the full jacobian should be a nested matrix, and I am wondering what is the best way to design the code. We are going to use DMDA's to manipulate Vecs for both variable sets, so the DMComposite approach of https://www.mcs.anl.gov/petsc/petsc-current/src/snes/tutorials/ex28.c.html is intriguing, but I have read in the comments that it has issues with MatNest type. My next guess would be to create the four submatrices ahead and then insert them in a MatNest, like in the Stokes example of https://www.mcs.anl.gov/petsc/petsc-current/src/snes/tutorials/ex70.c.html. However, in order to have shell blocks I guess it is almost mandatory to have the matrix partitioned among cpus as the Vecs are and I don't understand how Vecs end up being partitioned in ex70. We could - create a DMComposite and create the Vecs with it - get the local sizes of the Vecs and subVecs for the two variable groups - create the matrix as in ex70, using the shell type where/when needed, but instead of MatSetSizes(matblock, NULL, NULL, globalRows, globalCols) call MatSetSizes(matblock, localRows, localCols, NULL, NULL) using the local sizes of the subvectors. Does this sound a viable approach? Or do you have some different suggestions? Thanks Matteo
Matteo Semplice <[email protected]> writes:
Hi.
We are designing a PETSc application that will employ a SNES solver on a multiphysics problem whose jacobian will have a 2x2 block form, say A=[A00,A01;A10,A11]. We already have code for the top left block A_00 (a MatShell and a related Shell preconditioner) that we wish to reuse. We could implement the other blocks as Shells or assembled matrices. We'd like also to compare our method with existing ones, so we'd like to be quite flexible in the choice of KSP and PC within the SNES. (To this end, implementing an assembled version of the A00 and the other blocks would be easy)
I am assuming that, in order to have one or more shell blocks, the full jacobian should be a nested matrix, and I am wondering what is the best way to design the code.
We are going to use DMDA's to manipulate Vecs for both variable sets, so the DMComposite approach of https://www.mcs.anl.gov/petsc/petsc-current/src/snes/tutorials/ex28.c.html is intriguing, but I have read in the comments that it has issues with MatNest type.
I think ex28 is better organization of code. You can DMCreateMatrix() and then set types/preallocation for off-diagonal blocks of the MatNest. I think the comment is unclear and not quite what was intended and originally worked (which was to assemble the off-diagonal blocks despite bad preallocation). https://gitlab.com/petsc/petsc/-/commit/6bdeb4dbc27a59cf9af4930e08bd1f9937e4... https://www.mcs.anl.gov/petsc/petsc-current/src/snes/tutorials/ex28.c.html#l... Note that if you're using DMDA and have collocated fields, you can skip all this complexity. And if you have a scattered discretization, consider DMStag. ex28 is showing how to solve a coupled problem where there is no suitable structure to convey the relation between discretizations.
My next guess would be to create the four submatrices ahead and then insert them in a MatNest, like in the Stokes example of https://www.mcs.anl.gov/petsc/petsc-current/src/snes/tutorials/ex70.c.html. However, in order to have shell blocks I guess it is almost mandatory to have the matrix partitioned among cpus as the Vecs are and I don't understand how Vecs end up being partitioned in ex70.
Il 01/07/21 17:52, Jed Brown ha scritto:
I think ex28 is better organization of code. You can DMCreateMatrix() and then set types/preallocation for off-diagonal blocks of the MatNest. I think the comment is unclear and not quite what was intended and originally worked (which was to assemble the off-diagonal blocks despite bad preallocation).
https://gitlab.com/petsc/petsc/-/commit/6bdeb4dbc27a59cf9af4930e08bd1f9937e4...
https://www.mcs.anl.gov/petsc/petsc-current/src/snes/tutorials/ex28.c.html#l...
Thanks! Yesterday I was unable to make it work, but I'll have another go with ex28 then...
Note that if you're using DMDA and have collocated fields, you can skip all this complexity. And if you have a scattered discretization, consider DMStag. ex28 is showing how to solve a coupled problem where there is no suitable structure to convey the relation between discretizations.
The current discretization is in fact colocated with the variables on the same grid, and we might stick to that for a while. However, one key point in the design is that the jacobian will be [A00,A01;A10,A11] and we already have a taylor-made shell preconditioner for A00 and A00 implemented as a shell matrix; the preconditioner for the full Jacobian will be to neglect the A10 block and do a block-triangular solve inverting A00 approximately with the shell preconditioner. I do not understand how creating a DMDA with n0+n1 dofs will let me easily reuse my shell preconditioner code on the top-left block. Matteo
On Thu, Jul 1, 2021 at 11:44 AM Matteo Semplice < [email protected]> wrote:
Il 01/07/21 17:52, Jed Brown ha scritto:
I think ex28 is better organization of code. You can DMCreateMatrix() and then set types/preallocation for off-diagonal blocks of the MatNest. I think the comment is unclear and not quite what was intended and originally worked (which was to assemble the off-diagonal blocks despite bad preallocation).
https://gitlab.com/petsc/petsc/-/commit/6bdeb4dbc27a59cf9af4930e08bd1f9937e4...
https://www.mcs.anl.gov/petsc/petsc-current/src/snes/tutorials/ex28.c.html#l...
Thanks! Yesterday I was unable to make it work, but I'll have another go with ex28 then...
Note that if you're using DMDA and have collocated fields, you can skip all this complexity. And if you have a scattered discretization, consider DMStag. ex28 is showing how to solve a coupled problem where there is no suitable structure to convey the relation between discretizations.
The current discretization is in fact colocated with the variables on the same grid, and we might stick to that for a while.
However, one key point in the design is that the jacobian will be [A00,A01;A10,A11] and we already have a taylor-made shell preconditioner for A00 and A00 implemented as a shell matrix; the preconditioner for the full Jacobian will be to neglect the A10 block and do a block-triangular solve inverting A00 approximately with the shell preconditioner.
Okay, if that is the case, then you should use DMDA to layout the Jacobian, with DMCreateMatrix(), and then PCFIELDSPLIT to preconditioner, since it will automatically split things into your two pieces, and you can use your custom PC for A00 and multiplicative to get your upper triangular PC, or you could use a Schur complement to see if stronger coupling was more effective. Thanks, Matt
I do not understand how creating a DMDA with n0+n1 dofs will let me easily reuse my shell preconditioner code on the top-left block.
Matteo
-- 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/>
On Jul 1, 2021, at 11:44 AM, Matteo Semplice <[email protected]> wrote:
Il 01/07/21 17:52, Jed Brown ha scritto:
I think ex28 is better organization of code. You can DMCreateMatrix() and then set types/preallocation for off-diagonal blocks of the MatNest. I think the comment is unclear and not quite what was intended and originally worked (which was to assemble the off-diagonal blocks despite bad preallocation).
https://gitlab.com/petsc/petsc/-/commit/6bdeb4dbc27a59cf9af4930e08bd1f9937e4...
https://www.mcs.anl.gov/petsc/petsc-current/src/snes/tutorials/ex28.c.html#l...
Thanks! Yesterday I was unable to make it work, but I'll have another go with ex28 then...
Note that if you're using DMDA and have collocated fields, you can skip all this complexity. And if you have a scattered discretization, consider DMStag. ex28 is showing how to solve a coupled problem where there is no suitable structure to convey the relation between discretizations.
The current discretization is in fact colocated with the variables on the same grid, and we might stick to that for a while.
However, one key point in the design is that the jacobian will be [A00,A01;A10,A11] and we already have a taylor-made shell preconditioner for A00 and A00 implemented as a shell matrix; the preconditioner for the full Jacobian will be to neglect the A10 block and do a block-triangular solve inverting A00 approximately with the shell preconditioner.
I do not understand how creating a DMDA with n0+n1 dofs will let me easily reuse my shell preconditioner code on the top-left block.
PCFIELDSPLIT (and friends) do not order the dof by block, rather they "pull out" the required pieces of the vector (using IS's) when needed. Your shell preconditioner will just operate on the "pulled out" vectors. If you use DMDAVecGetArray etc in your shell preconditioner you can create an auxiliary DMDA of that smaller dof to still be able to use the DMDAVecGetArray constructs. You definitely should use an "all dof" DMDA to define your entire problem and not try to "glue" together vectors using DMComposites or other such things. Taking apart is much easier in parallel computing then putting together. Barry
Matteo
Thank you, Matthew and Barry! I can now see a way forward. Il 01/07/21 21:42, Barry Smith ha scritto:
I do not understand how creating a DMDA with n0+n1 dofs will let me easily reuse my shell preconditioner code on the top-left block. PCFIELDSPLIT (and friends) do not order the dof by block, rather they "pull out" the required pieces of the vector (using IS's) when needed. Your shell preconditioner will just operate on the "pulled out" vectors. If you use DMDAVecGetArray etc in your shell preconditioner you can create an auxiliary DMDA of that smaller dof to still be able to use the DMDAVecGetArray constructs. Just to be sure:
- I create a DMDA with n0+n1 dof per node - the jacobian will be associated to this DMDA. (It is not crucial, but can this be a shell matrix?) - I create a multiplicative PCfieldsplit, assign the correct n0 and n1 fields to each split (and get the IS for the splits via PCFieldSplitGetIS, shuld I need them) - the routine A00PCApply for the shell preconditioner of the A00 block, will see a Vec which is really a subvector with n0 dofs per node. In order to use DMDA semantics on this one, I create a DMDA with n0 dofs using DMDACreateCompatibleDMDA and then VecGetArrayDOFS using the smaller DMDA? Best Matteo
Sounds good. Yes the outer most Jacobian can be shell (or even nest but then I think you need to "build" it yourself, I don't think the DMDA will give back an appropriate nest matrix.) > On Jul 1, 2021, at 4:10 PM, Matteo Semplice <[email protected]> wrote: > > Thank you, Matthew and Barry! > > I can now see a way forward. > > Il 01/07/21 21:42, Barry Smith ha scritto: >>> I do not understand how creating a DMDA with n0+n1 dofs will let me easily reuse my shell preconditioner code on the top-left block. >> PCFIELDSPLIT (and friends) do not order the dof by block, rather they "pull out" the required pieces of the vector (using IS's) when needed. Your shell preconditioner will just operate on the "pulled out" vectors. If you use DMDAVecGetArray etc in your shell preconditioner you can create an auxiliary DMDA of that smaller dof to still be able to use the DMDAVecGetArray constructs. > Just to be sure: > - I create a DMDA with n0+n1 dof per node > > - the jacobian will be associated to this DMDA. (It is not crucial, but can this be a shell matrix?) > > - I create a multiplicative PCfieldsplit, assign the correct n0 and n1 fields to each split (and get the IS for the splits via PCFieldSplitGetIS, shuld I need them) > > - the routine A00PCApply for the shell preconditioner of the A00 block, will see a Vec which is really a subvector with n0 dofs per node. In order to use DMDA semantics on this one, I create a DMDA with n0 dofs using DMDACreateCompatibleDMDA and then VecGetArrayDOFS using the smaller DMDA? > > Best > > Matteo > > >
participants (4)
-
Barry Smith -
Jed Brown -
Matteo Semplice -
Matthew Knepley