<div dir="ltr"><div dir="ltr">On Mon, Jan 23, 2023 at 9:18 AM Jonas Lundgren via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>
<div lang="EN-US">
<div>
<p class="MsoNormal"><span lang="SV">Hi!<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="SV"><u></u> <u></u></span></p>
<p class="MsoNormal">(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.)<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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:<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">// Preliminaries<u></u><u></u></p>
<p class="MsoNormal">KSP ksp;<u></u><u></u></p>
<p class="MsoNormal">PC pc;<u></u><u></u></p>
<p class="MsoNormal">KSPCreate(PETSC_COMM_WORLD, ksp);<u></u><u></u></p>
<p class="MsoNormal">KSPSetType(ksp, KSPBCGS);<u></u><u></u></p>
<p class="MsoNormal">KSPSetFromOptions(ksp); // here, “-pc_type redistribute” is read from command line<u></u><u></u></p>
<p class="MsoNormal">KSPSetOperators(ksp, K, K);<u></u><u></u></p>
<p class="MsoNormal">KSPGetPC(ksp, &pc);<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">// Preconditioner<u></u><u></u></p>
<p class="MsoNormal">PCSetType(pc, PCFIELDSPLIT);<u></u><u></u></p>
<p class="MsoNormal">PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);<u></u><u></u></p>
<p class="MsoNormal">PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_LOWER);<u></u><u></u></p>
<p class="MsoNormal">PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELFP, K);<u></u><u></u></p>
<p class="MsoNormal">PCFieldSplitSetIS(pc,"0",isu);<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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.)<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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.</p></div></div></div></blockquote><div><br></div><div>Let me reply to the substance of the question. I will restate it to make sure I understand. You would like to eliminate the diagonal rows/cols from your operator before the solve, which PCREDISTRIBUTE can do. So you could have</div><div><br></div><div> -ksp_type preonly -pc_type redistribute</div><div> -redistribute_ksp_type bcgs -redistribute_pc_type fieldsplit</div><div> -redistribute_pc_fieldsplit_detect_saddle_point</div><div> -redistribute_pc_fieldsplit_type schur -redistribute_pc_fieldsplit_schur_fact_type lower</div><div> -redistribute_pc_fieldsplit_schur_precondition selfp</div><div><br></div><div>Using options solves the ordering problems you are having in setup and this is what we recommend. Satisfying the ordering constraints by hand, rather than programmatically, is very hard.</div><div><br></div><div>Also note the selfp is not a good preconditioner for Stokes. The Schur complement looks like</div><div><br></div><div> div (lap)^{-1} grad ~ id</div><div><br></div><div>so you want the mass matrix as the preconditioning matrix, not the diagonally weighted pressure Laplacian, which is what you get with selfp. What you can do is make a preconditioning matrix for the whole problem that looks like</div><div><br></div><div> / A B \</div><div> \ C M /</div><div><br></div><div>where M is the mass matrix, or leave out B and C and use -pc_amat to cleverly get them from the system matrix. You can see this in SNES ex69.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div lang="EN-US"><div><p class="MsoNormal"> 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:<u></u><u></u></p>
<p class="MsoNormal">PC ipc;<u></u><u></u></p>
<p class="MsoNormal">PCRedistributeGetKSP(pc, &iksp);<u></u><u></u></p>
<p class="MsoNormal">KSPGetPC(iksp, &ipc);<u></u><u></u></p>
<p class="MsoNormal">and letting “ipc” replace “pc” in the second code block (“Preconditioner”).<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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).
<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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.<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">I think I can solve the matrix problem by inserting the following (instead of the second to last row above):<u></u><u></u></p>
<p class="MsoNormal">PCFieldSplitSetSchurPre(ipc,PC_FIELDSPLIT_SCHUR_PRE_SELFP, ipc->mat);<u></u><u></u></p>
<p class="MsoNormal">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?<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">(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…)<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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:<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: Object is in wrong state<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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):<u></u><u></u></p>
<p class="MsoNormal">[1]PETSC ERROR: Argument out of range<u></u><u></u></p>
<p class="MsoNormal">[1]PETSC ERROR: Index 0's value 3864 is larger than maximum given 2334<u></u><u></u></p>
<p class="MsoNormal">[1]PETSC ERROR: #1 ISComplement() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/vec/is/is/utils/iscoloring.c:804<u></u><u></u></p>
<p class="MsoNormal">[1]PETSC ERROR: #2 PCFieldSplitSetDefaults() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/fieldsplit/fieldsplit.c:544<u></u><u></u></p>
<p class="MsoNormal">[1]PETSC ERROR: #3 PCSetUp_FieldSplit() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/fieldsplit/fieldsplit.c:588<u></u><u></u></p>
<p class="MsoNormal">[1]PETSC ERROR: #4 PCSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/interface/precon.c:994<u></u><u></u></p>
<p class="MsoNormal">[1]PETSC ERROR: #5 KSPSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/ksp/interface/itfunc.c:406<u></u><u></u></p>
<p class="MsoNormal">[1]PETSC ERROR: #6 PCSetUp_Redistribute() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/redistribute/redistribute.c:240<u></u><u></u></p>
<p class="MsoNormal">[1]PETSC ERROR: #7 PCSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/interface/precon.c:994<u></u><u></u></p>
<p class="MsoNormal">[1]PETSC ERROR: #8 KSPSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/ksp/interface/itfunc.c:406<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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 <a href="https://petsc.org/main/src/ksp/pc/impls/redistribute/redistribute.c.html" target="_blank">
https://petsc.org/main/src/ksp/pc/impls/redistribute/redistribute.c.html</a> <u></u>
<u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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:<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: Petsc has generated inconsistent data<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: Unhandled case, must have at least two fields, not 1<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: #1 PCFieldSplitSetDefaults() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/fieldsplit/fieldsplit.c:550<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: #2 PCSetUp_FieldSplit() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/fieldsplit/fieldsplit.c:588<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: #3 PCSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/interface/precon.c:994<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: #4 KSPSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/ksp/interface/itfunc.c:406<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: #5 PCSetUp_Redistribute() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/impls/redistribute/redistribute.c:240<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: #6 PCSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/pc/interface/precon.c:994<u></u><u></u></p>
<p class="MsoNormal">[0]PETSC ERROR: #7 KSPSetUp() at /mnt/c/Users/jonlu89/install/petsc-v3-18-1/src/ksp/ksp/interface/itfunc.c:406<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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”).<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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?<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">Is this even a sane way of looking at this problem?<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">Best regards,<u></u><u></u></p>
<p class="MsoNormal">Jonas Lundgren<u></u><u></u></p>
</div>
</div>
</div></blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>