<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Tue, Mar 21, 2017 at 12:58 PM, Natacha BEREUX <span dir="ltr"><<a href="mailto:natacha.bereux@gmail.com" target="_blank">natacha.bereux@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div>Dear PETSc user's, <br></div>I am trying to solve a poroelasticity problem with an additional temperature equation. The problem is a 3 fields problem involving a displacement field (u), a pressure field (p) and a temperature field (t).<br></div><div>I have seen similar examples in <br><a href="http://www.mcs.anl.gov/papers/P2017-0112.pdf" target="_blank">http://www.mcs.anl.gov/papers/<wbr>P2017-0112.pdf</a><br></div><div>or in Matt's talk <a href="http://www.caam.rice.edu/~mk51/presentations/SIAMCSE13.pdf" target="_blank">http://www.caam.rice.edu/~<wbr>mk51/presentations/SIAMCSE13.<wbr>pdf</a><br></div><div>I would like to reproduce them, but I am encountering troubles whem trying to do so.<br></div></div></blockquote><div><br></div><div>Yes, my talk is not clear on this point. The option you want, -pc_fieldsplit_<num>_fields,</div><div>only works if you use a DM to describe the splits. Here is the code for DMDA, which assumes</div><div>a regular division:</div><div><br></div><div> <a href="https://bitbucket.org/petsc/petsc/src/d4e0040555dc16f2fb1f7e2e0304e363dcc11328/src/ksp/pc/impls/fieldsplit/fieldsplit.c?at=master&fileviewer=file-view-default#fieldsplit.c-287">https://bitbucket.org/petsc/petsc/src/d4e0040555dc16f2fb1f7e2e0304e363dcc11328/src/ksp/pc/impls/fieldsplit/fieldsplit.c?at=master&fileviewer=file-view-default#fieldsplit.c-287</a></div><div><br></div><div>and for DMPlex and DMShell,</div><div><br></div><div> <a href="https://bitbucket.org/petsc/petsc/src/d4e0040555dc16f2fb1f7e2e0304e363dcc11328/src/ksp/pc/impls/fieldsplit/fieldsplit.c?at=master&fileviewer=file-view-default#fieldsplit.c-330">https://bitbucket.org/petsc/petsc/src/d4e0040555dc16f2fb1f7e2e0304e363dcc11328/src/ksp/pc/impls/fieldsplit/fieldsplit.c?at=master&fileviewer=file-view-default#fieldsplit.c-330</a></div><div><br></div><div>We could, of course, do something like this by just merging the ISes, but then we lose the ability to do it recursively since</div><div>we would not know how to split them again. So we would have to carry along something that tells us that, and this is exactly</div><div>the job that the DM is doing.</div><div><br></div><div>I think the remedy is as easy as specifying a DMShell that has a PetscSection (DMSetDefaultSection) with your ordering, and</div><div>I think this is how Firedrake (<a href="http://www.firedrakeproject.org/">http://www.firedrakeproject.org/</a>) does it. However, I usually use a DMPlex which knows about my</div><div>mesh, so I am not sure if this strategy has any holes.</div><div><br></div><div>The PetscSection is just a model of your ordering. You have a domain of "points", which are usually pieces of the mesh, which</div><div>map to a # of dofs. For example, linear elements associate 1 dof to every vertex. The Section allows dofs for several fields to</div><div>be added. This (I think) is described in the manual.</div><div><br></div><div>Does this make sense?</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div></div><div>Here is how I proceed:<br><br></div><div>I have a monolithic matrix A stemming . <br></div><div>I build 3 index sets for u,p, and t in A. <br></div><div>Then I set up the KSP context :<br><br></div><div>call KSPCreate(PETSC_COMM_WORLD,<wbr>ksp,ierr)<br>call KSPSetOperators(ksp,A,A,ierr)<br>call KSPGetPC(ksp, pc, ierr)<br>call PCSetType(pc, PCFIELDSPLIT, ierr)<br>call PCFieldSplitSetIS(pc,'u',is_u, ierr)<br>call PCFieldSplitSetIS(pc,'p',is_p, ierr)<br>call PCFieldSplitSetIS(pc,'t',is_t, ierr)<br>call PCSetFromOptions(pc,ierr)<br>call KSPSetFromOptions(ksp,ierr)<br>call KSPSolve(ksp,b,x,ierr)<br><br></div><div>I run the code with the following options<br>-ksp_view<br>-log_view<br>-ksp_rtol 1.0e-5<br>-ksp_type fgmres<br>-pc_type fieldsplit<br>-pc_fieldsplit_type additive<br>-pc_fieldsplit_0_fields 0,1<br>-pc_fieldsplit_1_fields 2<br>-pc_fieldsplit_0_pc_type lu<br>-pc_fieldsplit_0_ksp_type preonly<br>-pc_fieldsplit_1_pc_type lu<br>-pc_fieldsplit_1_ksp_type preonly<br><br></div><div>I would like to group u and p fields in a "0" block and then temperature remains in "1" block. <br></div><div>(I start with direct solves in each blocks to check the block decomposition but I intend to do use iterative methods later, and more precisely to use a Schur fieldsplit preconditionner for the "0" block) <br></div><div>The output follows :<br><br>KSP Object: 1 MPI processes<br> type: fgmres<br> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement<br> GMRES: happy breakdown tolerance 1e-30<br> maximum iterations=10000, initial guess is zero<br> tolerances: relative=1e-05, absolute=1e-50, divergence=10000.<br> right preconditioning<br> using UNPRECONDITIONED norm type for convergence test<br>PC Object: 1 MPI processes<br> type: fieldsplit<br> FieldSplit with ADDITIVE composition: total splits = 3<br> Solver info for each split is in the following KSP objects:<br> Split number 0 Defined by IS<br> KSP Object: (fieldsplit_u_) 1 MPI processes<br> type: preonly<br> maximum iterations=10000, initial guess is zero<br> tolerances: relative=1e-05, absolute=1e-50, divergence=10000.<br> left preconditioning<br> using DEFAULT norm type for convergence test<br> PC Object: (fieldsplit_u_) 1 MPI processes<br> type: ilu<br> PC has not been set up so information may be incomplete<br> ILU: out-of-place factorization<br> 0 levels of fill<br> tolerance for zero pivot 2.22045e-14<br> matrix ordering: natural<br> linear system matrix = precond matrix:<br> Mat Object: (fieldsplit_u_) 1 MPI processes<br> type: seqaij<br> rows=60, cols=60<br> total: nonzeros=3600, allocated nonzeros=3600<br> total number of mallocs used during MatSetValues calls =0<br> using I-node routines: found 12 nodes, limit used is 5<br> Split number 1 Defined by IS<br> KSP Object: (fieldsplit_p_) 1 MPI processes<br> type: preonly<br> maximum iterations=10000, initial guess is zero<br> tolerances: relative=1e-05, absolute=1e-50, divergence=10000.<br> left preconditioning<br> using DEFAULT norm type for convergence test<br> PC Object: (fieldsplit_p_) 1 MPI processes<br> type: ilu<br> PC has not been set up so information may be incomplete<br> ILU: out-of-place factorization<br> 0 levels of fill<br> tolerance for zero pivot 2.22045e-14<br> matrix ordering: natural<br> linear system matrix = precond matrix:<br> Mat Object: (fieldsplit_p_) 1 MPI processes<br> type: seqaij<br> rows=8, cols=8<br> total: nonzeros=64, allocated nonzeros=64<br> total number of mallocs used during MatSetValues calls =0<br> using I-node routines: found 2 nodes, limit used is 5<br> Split number 2 Defined by IS<br> KSP Object: (fieldsplit_t_) 1 MPI processes<br> type: preonly<br> maximum iterations=10000, initial guess is zero<br> tolerances: relative=1e-05, absolute=1e-50, divergence=10000.<br> left preconditioning<br> using DEFAULT norm type for convergence test<br> PC Object: (fieldsplit_t_) 1 MPI processes<br> type: ilu<br> PC has not been set up so information may be incomplete<br> ILU: out-of-place factorization<br> 0 levels of fill<br> tolerance for zero pivot 2.22045e-14<br> matrix ordering: natural<br> linear system matrix = precond matrix:<br> Mat Object: (fieldsplit_t_) 1 MPI processes<br> type: seqaij<br> rows=8, cols=8<br> total: nonzeros=64, allocated nonzeros=64<br> total number of mallocs used during MatSetValues calls =0<br> using I-node routines: found 2 nodes, limit used is 5<br> linear system matrix = precond matrix:<br> Mat Object: 1 MPI processes<br> type: seqaij<br> rows=76, cols=76<br> total: nonzeros=5776, allocated nonzeros=5776<br> total number of mallocs used during MatSetValues calls =0<br> using I-node routines: found 16 nodes, limit used is 5<br><br></div><div>The preconditionner has 3 splits, whereas I would like to group (u,p) together and see 2 splits <br></div><div>I suspect that <br></div><div>-pc_fieldsplit_0_fields 0,1<br>-pc_fieldsplit_1_fields 2<br></div><div>are not the appropriate options.<br></div><div>Am I correct ? <br></div><div>What is the right way for grouping two fields defined by index sets ? <br></div><div><br></div><div>Any help would be greatly appreciated !<span class="gmail-HOEnZb"><font color="#888888"><br><br></font></span></div><span class="gmail-HOEnZb"><font color="#888888"><div>Natacha <br></div><br></font></span></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">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></div>