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
