[petsc-users] Using PCREDISTRIBUTE together with PCFIELDSPLIT

Fri Jan 20 18:21:20 CST 2023

Hi!

(Sorry for a long message, I have tried to cover the essentials only. I am happy to provide further details and logs if necessary, but I have tried to keep it as short as possible.)

I am trying to solve a Stokes flow problem with PCFIELDSPLIT as preconditioner (and KSPBCGS as solver). I have successfully set up the preconditioner and solved several examples by applying the following options when setting up the solver:

// Preliminaries
KSP ksp;
PC pc;
KSPCreate(PETSC_COMM_WORLD, ksp);
KSPSetType(ksp, KSPBCGS);
KSPSetFromOptions(ksp);  // here, "-pc_type redistribute" is read from command line
KSPSetOperators(ksp, K, K);
KSPGetPC(ksp, &pc);

// Preconditioner
PCSetType(pc, PCFIELDSPLIT);
PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_LOWER);
PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELFP, K);
PCFieldSplitSetIS(pc,"0",isu);

where "isu" in the last row is an IS containing all flow indices (no pressure indices), and "K" is my system matrix, which is created using DMCreateMatrix(dm_state, &K); and "dm_state" is a DMStag object. (However, K is not assembled yet.)

I want to try to use PCREDISTRIBUTE on top of this, since I believe that I have a lot of degrees-of-freedom (DOF) locked, and if I can reduce the size of my linear problem, the solution time can decrease as well. My idea was to use PCREDISTRIBUTE as the main preconditioner (with KSPPREONLY as the solver), and to move the KSPBCGS and PCFIELDSPLIT down one level, to act as a sub-KSP and sub-PC, by introducing the following between the two code blocks:
PC ipc;
PCRedistributeGetKSP(pc, &iksp);
KSPGetPC(iksp, &ipc);
and letting "ipc" replace "pc" in the second code block ("Preconditioner").

The two last rows of in the "Preconditioner" block refers to "K" and "isu", which are defined in terms of the entire system, not the reduced. I believe that I have to replace these two variables with a Mat ("A") that has locked DOFs removed, and an IS ("isfree") that contains only free flow DOF indices (no locked flow indices and no pressure indices).

And here comes my problem: the Mat "A" and the IS "isfree" will only be available when the original KSP ("ksp") - or possibly the PC "pc" - are set up (using KSPSetUp(ksp); or PCSetUp(pc); ). I have no way of knowing how "A" or "isfree" looks before it is automatically extracted from "K" and "isu" during set up - they will be different depending on the specifics of the example I am running.

I think I can solve the matrix problem by inserting the following (instead of the second to last row above):
PCFieldSplitSetSchurPre(ipc,PC_FIELDSPLIT_SCHUR_PRE_SELFP, ipc->mat);
Even though "ipc->mat" is not created yet, I believe that this might lead to the use of the extracted sub-flow-matrix as (a part of the) preconditioner for the flow block in PCFIELDSPLIT. I do not know for sure, since I have other issues that have prevented me from testing this approach. Does this sound reasonable?

(One might come up with a better alternative than the system matrix to use for this sub-PC, but that work is yet to be done...)

Furthermore, the second issue I have is regarding the IS "isfree". After the "Preconditioner" code block I am trying to set the Fieldsplit sub-PCs by calling PCFieldSplitSchurGetSubKSP(); It appears as if I need to set up the outer KSP "ksp" or PC "pc" before I call PCFieldSplitSchurGetSubKSP(); or otherwise I get an error message looking like:
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()

However, when I try to call KSPSetUp() as suggested above (which is possible thanks to assembling the system matrix before the call to set up the KSP), I just get another error message indicating that I have the wrong IS ("isu") (most likely because the indices "isu" refers to the entire system, not the reduced):
[1]PETSC ERROR: Argument out of range
[1]PETSC ERROR: Index 0's value 3864 is larger than maximum given 2334
[1]PETSC ERROR: #1 ISComplement() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/vec/is/is/utils/iscoloring.c:804
[1]PETSC ERROR: #2 PCFieldSplitSetDefaults() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/fieldsplit/fieldsplit.c:544
[1]PETSC ERROR: #3 PCSetUp_FieldSplit() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/fieldsplit/fieldsplit.c:588
[1]PETSC ERROR: #4 PCSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/interface/precon.c:994
[1]PETSC ERROR: #5 KSPSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/ksp/interface/itfunc.c:406
[1]PETSC ERROR: #6 PCSetUp_Redistribute() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/redistribute/redistribute.c:240
[1]PETSC ERROR: #7 PCSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/interface/precon.c:994
[1]PETSC ERROR: #8 KSPSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/ksp/interface/itfunc.c:406

Okay, so now I want to obtain the reduced IS "isfree" that contains only free flow DOFs (not locked flow DOFs or pressure DOFs). One option is to try to obtain it from the same IS that extracts the sub-matrix when setting up the PCREDISTRIBUTE preconditioner, see red->is in MatCreateSubMatrix() in https://petsc.org/main/src/ksp/pc/impls/redistribute/redistribute.c.html

And annoyingly, I need to have called PCSetUp(pc) or KSPSetUp(ksp) before using red->is (thus, before calling PCFieldSplitSetIS() ), which is something that is impossible it seems:
[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Unhandled case, must have at least two fields, not 1
[0]PETSC ERROR: #1 PCFieldSplitSetDefaults() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/fieldsplit/fieldsplit.c:550
[0]PETSC ERROR: #2 PCSetUp_FieldSplit() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/fieldsplit/fieldsplit.c:588
[0]PETSC ERROR: #3 PCSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/interface/precon.c:994
[0]PETSC ERROR: #4 KSPSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/ksp/interface/itfunc.c:406
[0]PETSC ERROR: #5 PCSetUp_Redistribute() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/redistribute/redistribute.c:240
[0]PETSC ERROR: #6 PCSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/interface/precon.c:994
[0]PETSC ERROR: #7 KSPSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/ksp/interface/itfunc.c:406

So, to summarize, I seem to be stuck in a Catch 22; I cannot set up the KSP ("ksp") before I have the correct IS ("isfree") set to my Fieldsplit precondtitioner, but I cannot obtain the correct IS preconditioner ("isfree") before I set up the (outer) KSP ("ksp").

Is there any possible ideas on how to resolve this issue? Do I have to hard code the IS "isfree" in order to achieve what I want? (I don't know how...) Or maybe this approach is wrong all together. Can I achieve my goal of using PCFIELDSPLIT and PCREDISTRIBUTE together in another way?

Is this even a sane way of looking at this problem?

Best regards,
Jonas Lundgren
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230121/8aaa1777/attachment-0001.html>