Re: [petsc-dev] MatNest and FieldSplit

[Pierre Jolivet via petsc-dev]

> On 15 Apr 2019, at 2:17 PM, Matthew Knepley  > wrote:
> 
> On Mon, Apr 15, 2019 at 5:03 AM Pierre Jolivet via petsc-dev 
> mailto:petsc-dev@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 
>> mailto:petsc-dev@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
>>  
>> 
>> 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
>>  
>> 
>>  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
>>  
>> 
>>  and 
>> 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 
>>> mailto:petsc-dev@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 >>> > 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.1]
 
 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 
 mailto:petsc-dev@mcs.anl.gov>>

Re: [petsc-dev] MatNest and FieldSplit

[Pierre Jolivet via petsc-dev]

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…?
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.

Thanks,
Pierre

> On 14 Apr 2019, at 10:54 PM, Pierre Jolivet via petsc-dev 
>  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
>  
> 
> 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
>  
> 
>  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
>  
> 
>  and 
> 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 
>> mailto:petsc-dev@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 >> > 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.1]
>>> 
>>> 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 
>>> mailto:petsc-dev@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.1]]
>>> 
>>> 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 
>>> mailto:petsc-dev@mcs.anl.gov>> wrote:
>>> On Sun, Mar 24, 2019 at 10:21 AM Pierre Jolivet >> > 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

Re: [petsc-dev] MatNest and FieldSplit

[Pierre Jolivet via petsc-dev]

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
 

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
 

 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
 

 and 
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 
>  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 > > 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.1]
>> 
>> 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 
>> mailto:petsc-dev@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.1]]
>> 
>> 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 
>> mailto:petsc-dev@mcs.anl.gov>> wrote:
>> On Sun, Mar 24, 2019 at 10:21 AM Pierre Jolivet > > 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 >> > wrote:
>>> 
>>> On Sat, Mar 23, 2019 at 9:12 PM Pierre Jolivet via petsc-dev 
>>> mailto:petsc-dev@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 interes

Re: [petsc-dev] MatNest and FieldSplit

[Pierre Jolivet via petsc-dev]

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  
> 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.1]
> 
> 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  > 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.1]]
> 
> 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 
> mailto:petsc-dev@mcs.anl.gov>> wrote:
> On Sun, Mar 24, 2019 at 10:21 AM Pierre Jolivet  > 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 > > wrote:
>> 
>> On Sat, Mar 23, 2019 at 9:12 PM Pierre Jolivet via petsc-dev 
>> mailto:petsc-dev@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/ 
> 
> 
> 
> -- 
> 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/ 

Re: [petsc-dev] MatNest and FieldSplit

[Mark Adams via petsc-dev]

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.1]]


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@mcs.anl.gov> wrote:

> On Sun, Mar 24, 2019 at 10:21 AM Pierre Jolivet <
> pierre.joli...@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  wrote:
>>
>> On Sat, Mar 23, 2019 at 9:12 PM Pierre Jolivet via petsc-dev <
>> petsc-dev@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/
>> 
>>
>>
>>
>
> --
> 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/
> 
>

Re: [petsc-dev] MatNest and FieldSplit

[Pierre Jolivet via petsc-dev]

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 don’t know if this makes more or less sense… sorry :\
Thanks,
Pierre

> On 24 Mar 2019, at 8:42 PM, Matthew Knepley  wrote:
> 
> On Sat, Mar 23, 2019 at 9:12 PM Pierre Jolivet via petsc-dev 
> mailto:petsc-dev@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/ 

[petsc-dev] MatNest and FieldSplit

[Pierre Jolivet via petsc-dev]

I’m trying to figure out why both solutions are not consistent in the following 
example.
Is what I’m doing complete nonsense?

Thanks in advance for your help,
Pierre



nest.c
Description: Binary data