[petsc-dev] MatNest and FieldSplit

Pierre Jolivet pierre.jolivet at enseeiht.fr
Mon Apr 15 07:57:15 CDT 2019 


> On 15 Apr 2019, at 2:17 PM, Matthew Knepley <knepley at gmail.com <mailto:knepley at gmail.com>> wrote:
> 
> On Mon, Apr 15, 2019 at 5:03 AM Pierre Jolivet via petsc-dev <petsc-dev at mcs.anl.gov <mailto:petsc-dev at mcs.anl.gov>> wrote:
> OK, my solver is now working properly with a call to
> PetscObjectCompose((PetscObject)isS, "pmat", (PetscObject) myS);
> 
> I have two follow-up questions:
> 1) am I supposed to call this, or is it the sign of something done wrong in my sequence of SNESSolve/KSPSetUp/KSPSetFromOptions/KSPSetOperators…?
> 
> Yes, unfortunately. We want to let the user pass in preconditioning matrices for arbitrarily nested splits without pulling out
> the objects explicitly, so we need a way to locate a particular split. I guess you could imagine using the prefix, which would
> be an XML way of doing things. We are doing it by attaching things to the IS, since the Solver/Mat objects are created on the
> fly. In fact, the IS is created on the fly by a DM sometimes, so you can attach to the DM (I do this for Plex).

Indeed, I found — to my surprise one of the very few occurence of PetscObjectCompose((PetscObject)isS, “pmat”, …) — in dm/interface/dmi.c.
And actually, I realised later on that even for a “simple” example like snes/examples/tutorials/ex70.c with -user_ksp turned on, the call to KSPSetOperators(subksp[1], s->myS, s->myS) is not sufficient to avoid the call to MatCreateSubMatrix on the Schur complement. Which means (at least) this part of the solver is behaving like it should, thanks for making it clear to me.

> 2) I have currently hardwired the split name I’m using when calling PCFieldSplitGetIS to get “isS” (from above) for debugging purposes. Could I create a PR with a new function like PCFieldSplitGetISByIndex(PC pc,const PetscInt n,IS *is) that will return the nth IS? Right now, I would need to get the KSP prefix followed up by some string comparison to get the actual IS prefix, whereas I know the position of the KSP in the PC_FieldSplitLink.
> 
> That sounds fine to me.

Done (PR #1544).

Thanks,
Pierre

>    Matt
>  
> Thanks,
> Pierre
> 
>> On 14 Apr 2019, at 10:54 PM, Pierre Jolivet via petsc-dev <petsc-dev at mcs.anl.gov <mailto:petsc-dev at mcs.anl.gov>> wrote:
>> 
>> I think I figured out what my problem is _exactly_.
>> The Mat inside the MATNEST on which I’m using a PCFIELDSPLIT is unassembled before the first KSPSolve, except for the last field.
>> Matt, you nailed it, when I call KSPSetFromOptions on the global PCFIELDSPLIT, then KSPSetUp explicitly on the inner PCFIELDSPLIT, the matrices associated to all fields are created here: https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/pc/impls/fieldsplit/fieldsplit.c.html#line656 <https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/pc/impls/fieldsplit/fieldsplit.c.html#line656>
>> Then, whatever I do with the matrix from the last field, currently trying MatCopy(myAssembledS, generatedS), before the first KSPSolve is reset by https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/pc/impls/fieldsplit/fieldsplit.c.html#line688 <https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/pc/impls/fieldsplit/fieldsplit.c.html#line688> and my solver fails…
>> 
>> So pretty simple question, how do I set a “pmat” for my last assembled field so that https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/pc/impls/fieldsplit/fieldsplit.c.html#line646 <https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/pc/impls/fieldsplit/fieldsplit.c.html#line646> and https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/pc/impls/fieldsplit/fieldsplit.c.html#line686 <https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/pc/impls/fieldsplit/fieldsplit.c.html#line686> will return a non null pmat.
>> 
>> This may sound really trivial but I’m lost in limbo right now. When everything is not wrapped inside an outer PCFIELDSPLIT, everything just work.
>> 
>> Thanks,
>> Pierre
>> 
>>> On 25 Mar 2019, at 6:57 PM, Pierre Jolivet via petsc-dev <petsc-dev at mcs.anl.gov <mailto:petsc-dev at mcs.anl.gov>> wrote:
>>> 
>>> Thanks, this makes (slightly) more sense to me know.
>>> For some reason my application is still not acting properly but I must be screwing somewhere else the nested FieldSplit…
>>> 
>>> Thank you,
>>> Pierre
>>> 
>>>> On 24 Mar 2019, at 11:42 PM, Dave May via petsc-dev <petsc-dev at mcs.anl.gov <mailto:petsc-dev at mcs.anl.gov>> wrote:
>>>> 
>>>> Matt is right.
>>>> 
>>>> When you defined the operator S, you basically invalidate the operator N (in the sense that they are no longer consistent). Hence when you use KSP nest to solve your problem your A matrix looks like 
>>>>   A = diag[1, 2, 4, 0, 8]
>>>> but the B matrix you have defined looks like
>>>>   B = diag[1, 2, 4, 0.00001]
>>>> 
>>>> The only way to obtain the correct answer with your code is thus to use the option
>>>> -ksp_type preonly
>>>> 
>>>> Thanks
>>>> Dave
>>>> 
>>>> 
>>>> 
>>>> On Sun, 24 Mar 2019 at 22:09, Mark Adams via petsc-dev <petsc-dev at mcs.anl.gov <mailto:petsc-dev at mcs.anl.gov>> wrote:
>>>> I think he is saying that this line seems to have no effect (and the comment is hence wrong):
>>>> 
>>>> KSPSetOperators(subksp[nsplits - 1], S, S);
>>>> // J2 = [[4, 0] ; [0, 0.00001]]
>>>> 
>>>> J2 is a 2x2 but this block has been changed into two single equation fields. Does this KSPSetOperators supposed to copy this 1x1 S matrix into the (1,1) block of the "J2", or do some sort of correct mixing internally, to get what he wants?
>>>> 
>>>> BTW, this line does not seem necessary to me so maybe I'm missing something.
>>>> 
>>>> KSPSetOperators(sub, J2, J2);
>>>> 
>>>> 
>>>> On Sun, Mar 24, 2019 at 4:33 PM Matthew Knepley via petsc-dev <petsc-dev at mcs.anl.gov <mailto:petsc-dev at mcs.anl.gov>> wrote:
>>>> On Sun, Mar 24, 2019 at 10:21 AM Pierre Jolivet <pierre.jolivet at enseeiht.fr <mailto:pierre.jolivet at enseeiht.fr>> wrote:
>>>> It’s a 4x4 matrix.
>>>> The first 2x2 diagonal matrix is a field.
>>>> The second 2x2 diagonal matrix is another field.
>>>> In the second field, the first diagonal coefficient is a subfield.
>>>> In the second field, the second diagonal coefficient is another subfield.
>>>> I’m changing the operators from the second subfield (last diagonal coefficient of the matrix).
>>>> When I solve a system with the complete matrix (2 fields), I get a different “partial solution" than when I solve the “partial system” on just the second field (with the two subfields in which I modified the operators from the second one).
>>>> 
>>>> I may understand waht you are doing.
>>>> Fieldsplit calls MatGetSubMatrix() which can copy values, depending on the implementation,
>>>> so changing values in the original matrix may or may not change it in the PC.
>>>>  
>>>>    Matt
>>>> 
>>>> I don’t know if this makes more or less sense… sorry :\
>>>> Thanks,
>>>> Pierre
>>>> 
>>>>> On 24 Mar 2019, at 8:42 PM, Matthew Knepley <knepley at gmail.com <mailto:knepley at gmail.com>> wrote:
>>>>> 
>>>>> On Sat, Mar 23, 2019 at 9:12 PM Pierre Jolivet via petsc-dev <petsc-dev at mcs.anl.gov <mailto:petsc-dev at mcs.anl.gov>> wrote:
>>>>> I’m trying to figure out why both solutions are not consistent in the following example.
>>>>> Is what I’m doing complete nonsense?
>>>>> 
>>>>> The code does not make clear what you are asking. I can see its a nested fieldsplit.
>>>>> 
>>>>>   Thanks,
>>>>> 
>>>>>      Matt
>>>>>  
>>>>> Thanks in advance for your help,
>>>>> Pierre
>>>>> 
>>>>> 
>>>>> 
>>>>> -- 
>>>>> 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/>
>>>> 
>>>> 
>>>> 
>>>> -- 
>>>> 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/>
>>> 
>> 
> 
> 
> 
> -- 
> 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20190415/89a93c4f/attachment.html>

More information about the petsc-dev mailing list