Re: [petsc-users] Scalable Solver for Incompressible Flow

[Jed Brown]

I would email the author. He's been helpful in the past and that newer paper 
may have been extensions that didn't make it into the upstream example.

Alexander Lindsay  writes:

> Maybe that example was a jumping point for some of those studies, but it
> looks to me like that example has been around since ~2012 and initially
> only touched on SIMPLE, as opposed to addition of SIMPLE into an
> augmented lagrange scheme.
>
> But it does seem that at some point Carola Kruse added Golub-Kahan
> bidiagonalization tests to ex70. I don't know very much about that although
> it seems to be related to AL methods ... but requires that the matrix be
> symmetric?
>
> On Fri, Jul 28, 2023 at 7:04 PM Jed Brown  wrote:
>
>> See src/snes/tutorials/ex70.c for the code that I think was used for that
>> paper.
>>
>> Alexander Lindsay  writes:
>>
>> > Sorry for the spam. Looks like these authors have published multiple
>> papers on the subject
>> >
>> > cover.jpg
>> > Combining the Augmented Lagrangian Preconditioner with the Simple Schur
>> Complement Approximation | SIAM Journal on
>> > Scientific Computingdoi.org
>> > cover.jpg
>> >
>> >  On Jul 28, 2023, at 12:59 PM, Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>> >
>> >  Do you know of anyone who has applied the augmented Lagrange
>> methodology to a finite volume discretization?
>> >
>> >  On Jul 6, 2023, at 6:25 PM, Matthew Knepley  wrote:
>> >
>> >  On Thu, Jul 6, 2023 at 8:30 PM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>> >
>> >  This is an interesting article that compares a multi-level ILU
>> algorithm to approximate commutator and augmented
>> >  Lagrange methods: https://doi.org/10.1002/fld.5039
>> >
>> >  That is for incompressible NS. The results are not better than
>> https://arxiv.org/abs/1810.03315, and that PC is considerably
>> >  simpler and already implemented in PETSc. There is an update in to this
>> >
>> >
>> >
>> https://epubs.siam.org/doi/abs/10.1137/21M1430698?casa_token=Fp_XhuZStZ0A:YDhnkW9XvAom_b8KocWz-hBEI7FAt46aw3ICa0FvCrOVCtYr9bwvtqJ4aBOTkDSvANKh6YTQEw
>> >
>> >
>> >  which removes the need for complicated elements.
>> >
>> >  You might need stuff like ILU for compressible flow, but I think
>> incompressible is solved.
>> >
>> >Thanks,
>> >
>> >   Matt
>> >
>> >  On Wed, Jun 28, 2023 at 11:37 AM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>> >
>> >  I do believe that based off the results in
>> https://doi.org/10.1137/040608817 we should be able to make LSC, with
>> >  proper scaling, compare very favorably with PCD
>> >
>> >  On Tue, Jun 27, 2023 at 10:41 AM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>> >
>> >  I've opened https://gitlab.com/petsc/petsc/-/merge_requests/6642 which
>> adds a couple more scaling
>> >  applications of the inverse of the diagonal of A
>> >
>> >  On Mon, Jun 26, 2023 at 6:06 PM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>> >
>> >  I guess that similar to the discussions about selfp, the approximation
>> of the velocity mass matrix by the
>> >  diagonal of the velocity sub-matrix will improve when running a
>> transient as opposed to a steady
>> >  calculation, especially if the time derivative is lumped Just
>> thinking while typing
>> >
>> >  On Mon, Jun 26, 2023 at 6:03 PM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>> >
>> >  Returning to Sebastian's question about the correctness of the current
>> LSC implementation: in the
>> >  taxonomy paper that Jed linked to (which talks about SIMPLE, PCD, and
>> LSC), equation 21 shows four
>> >  applications of the inverse of the velocity mass matrix. In the PETSc
>> implementation there are at
>> >  most two applications of the reciprocal of the diagonal of A (an
>> approximation to the velocity mass
>> >  matrix without more plumbing, as already pointed out). It seems like
>> for code implementations in
>> >  which there are possible scaling differences between the velocity and
>> pressure equations, that this
>> >  difference in the number of inverse applications could be significant?
>> I know Jed said that these
>> >  scalings wouldn't really matter if you have a uniform grid, but I'm not
>> 100% convinced yet.
>> >
>> >  I might try fiddling around with adding two more reciprocal
>> applications.
>> >
>> >  On Fri, Jun 23, 2023 at 1:09 PM Pierre Jolivet 
>> wrote:
>> >
>> >  On 23 Jun 2023, at 10:06 PM, Pierre Jolivet 
>> wrote:
>> >
>> >  On 23 Jun 2023, at 9:39 PM, Alexander Lindsay 
>> wrote:
>> >
>> >  Ah, I see that if I use Pierre's new 'full' option for
>> -mat_schur_complement_ainv_type
>> >
>> >  That was not initially done by me
>> >
>> >  Oops, sorry for the noise, looks like it was done by me indeed in
>> >  9399e4fd88c6621aad8fe9558ce84df37bd6fada…
>> >
>> >  Thanks,
>> >  Pierre
>> >
>> >  (though I recently tweaked MatSchurComplementComputeExplicitOperator()
>> a bit to use
>> >  KSPMatSolve(), so that if you have a small Schur complement — which is
>> 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

Maybe that example was a jumping point for some of those studies, but it
looks to me like that example has been around since ~2012 and initially
only touched on SIMPLE, as opposed to addition of SIMPLE into an
augmented lagrange scheme.

But it does seem that at some point Carola Kruse added Golub-Kahan
bidiagonalization tests to ex70. I don't know very much about that although
it seems to be related to AL methods ... but requires that the matrix be
symmetric?

On Fri, Jul 28, 2023 at 7:04 PM Jed Brown  wrote:

> See src/snes/tutorials/ex70.c for the code that I think was used for that
> paper.
>
> Alexander Lindsay  writes:
>
> > Sorry for the spam. Looks like these authors have published multiple
> papers on the subject
> >
> > cover.jpg
> > Combining the Augmented Lagrangian Preconditioner with the Simple Schur
> Complement Approximation | SIAM Journal on
> > Scientific Computingdoi.org
> > cover.jpg
> >
> >  On Jul 28, 2023, at 12:59 PM, Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
> >
> >  Do you know of anyone who has applied the augmented Lagrange
> methodology to a finite volume discretization?
> >
> >  On Jul 6, 2023, at 6:25 PM, Matthew Knepley  wrote:
> >
> >  On Thu, Jul 6, 2023 at 8:30 PM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
> >
> >  This is an interesting article that compares a multi-level ILU
> algorithm to approximate commutator and augmented
> >  Lagrange methods: https://doi.org/10.1002/fld.5039
> >
> >  That is for incompressible NS. The results are not better than
> https://arxiv.org/abs/1810.03315, and that PC is considerably
> >  simpler and already implemented in PETSc. There is an update in to this
> >
> >
> >
> https://epubs.siam.org/doi/abs/10.1137/21M1430698?casa_token=Fp_XhuZStZ0A:YDhnkW9XvAom_b8KocWz-hBEI7FAt46aw3ICa0FvCrOVCtYr9bwvtqJ4aBOTkDSvANKh6YTQEw
> >
> >
> >  which removes the need for complicated elements.
> >
> >  You might need stuff like ILU for compressible flow, but I think
> incompressible is solved.
> >
> >Thanks,
> >
> >   Matt
> >
> >  On Wed, Jun 28, 2023 at 11:37 AM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
> >
> >  I do believe that based off the results in
> https://doi.org/10.1137/040608817 we should be able to make LSC, with
> >  proper scaling, compare very favorably with PCD
> >
> >  On Tue, Jun 27, 2023 at 10:41 AM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
> >
> >  I've opened https://gitlab.com/petsc/petsc/-/merge_requests/6642 which
> adds a couple more scaling
> >  applications of the inverse of the diagonal of A
> >
> >  On Mon, Jun 26, 2023 at 6:06 PM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
> >
> >  I guess that similar to the discussions about selfp, the approximation
> of the velocity mass matrix by the
> >  diagonal of the velocity sub-matrix will improve when running a
> transient as opposed to a steady
> >  calculation, especially if the time derivative is lumped Just
> thinking while typing
> >
> >  On Mon, Jun 26, 2023 at 6:03 PM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
> >
> >  Returning to Sebastian's question about the correctness of the current
> LSC implementation: in the
> >  taxonomy paper that Jed linked to (which talks about SIMPLE, PCD, and
> LSC), equation 21 shows four
> >  applications of the inverse of the velocity mass matrix. In the PETSc
> implementation there are at
> >  most two applications of the reciprocal of the diagonal of A (an
> approximation to the velocity mass
> >  matrix without more plumbing, as already pointed out). It seems like
> for code implementations in
> >  which there are possible scaling differences between the velocity and
> pressure equations, that this
> >  difference in the number of inverse applications could be significant?
> I know Jed said that these
> >  scalings wouldn't really matter if you have a uniform grid, but I'm not
> 100% convinced yet.
> >
> >  I might try fiddling around with adding two more reciprocal
> applications.
> >
> >  On Fri, Jun 23, 2023 at 1:09 PM Pierre Jolivet 
> wrote:
> >
> >  On 23 Jun 2023, at 10:06 PM, Pierre Jolivet 
> wrote:
> >
> >  On 23 Jun 2023, at 9:39 PM, Alexander Lindsay 
> wrote:
> >
> >  Ah, I see that if I use Pierre's new 'full' option for
> -mat_schur_complement_ainv_type
> >
> >  That was not initially done by me
> >
> >  Oops, sorry for the noise, looks like it was done by me indeed in
> >  9399e4fd88c6621aad8fe9558ce84df37bd6fada…
> >
> >  Thanks,
> >  Pierre
> >
> >  (though I recently tweaked MatSchurComplementComputeExplicitOperator()
> a bit to use
> >  KSPMatSolve(), so that if you have a small Schur complement — which is
> not really the case
> >  for NS — this could be a viable option, it was previously painfully
> slow).
> >
> >  Thanks,
> >  Pierre
> >
> >  that I get a single iteration for the Schur complement solve with LU.
> That's a nice testing
> >  option
> >
> >  On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay <
> 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Jed Brown]

See src/snes/tutorials/ex70.c for the code that I think was used for that paper.

Alexander Lindsay  writes:

> Sorry for the spam. Looks like these authors have published multiple papers 
> on the subject 
>
> cover.jpg 
> Combining the Augmented Lagrangian Preconditioner with the Simple Schur 
> Complement Approximation | SIAM Journal on  
> Scientific Computingdoi.org  
> cover.jpg
>
>  On Jul 28, 2023, at 12:59 PM, Alexander Lindsay  
> wrote:
>
>  Do you know of anyone who has applied the augmented Lagrange methodology to 
> a finite volume discretization?
>
>  On Jul 6, 2023, at 6:25 PM, Matthew Knepley  wrote:
>
>  On Thu, Jul 6, 2023 at 8:30 PM Alexander Lindsay  
> wrote:
>
>  This is an interesting article that compares a multi-level ILU algorithm to 
> approximate commutator and augmented
>  Lagrange methods: https://doi.org/10.1002/fld.5039
>
>  That is for incompressible NS. The results are not better than 
> https://arxiv.org/abs/1810.03315, and that PC is considerably
>  simpler and already implemented in PETSc. There is an update in to this
>
>   
>  
> https://epubs.siam.org/doi/abs/10.1137/21M1430698?casa_token=Fp_XhuZStZ0A:YDhnkW9XvAom_b8KocWz-hBEI7FAt46aw3ICa0FvCrOVCtYr9bwvtqJ4aBOTkDSvANKh6YTQEw
>  
>
>  which removes the need for complicated elements.
>
>  You might need stuff like ILU for compressible flow, but I think 
> incompressible is solved.
>
>Thanks,
>
>   Matt
>   
>  On Wed, Jun 28, 2023 at 11:37 AM Alexander Lindsay 
>  wrote:
>
>  I do believe that based off the results in https://doi.org/10.1137/040608817 
> we should be able to make LSC, with
>  proper scaling, compare very favorably with PCD
>
>  On Tue, Jun 27, 2023 at 10:41 AM Alexander Lindsay 
>  wrote:
>
>  I've opened https://gitlab.com/petsc/petsc/-/merge_requests/6642 which adds 
> a couple more scaling
>  applications of the inverse of the diagonal of A
>
>  On Mon, Jun 26, 2023 at 6:06 PM Alexander Lindsay  
> wrote:
>
>  I guess that similar to the discussions about selfp, the approximation of 
> the velocity mass matrix by the
>  diagonal of the velocity sub-matrix will improve when running a transient as 
> opposed to a steady
>  calculation, especially if the time derivative is lumped Just thinking 
> while typing
>
>  On Mon, Jun 26, 2023 at 6:03 PM Alexander Lindsay  
> wrote:
>
>  Returning to Sebastian's question about the correctness of the current LSC 
> implementation: in the
>  taxonomy paper that Jed linked to (which talks about SIMPLE, PCD, and LSC), 
> equation 21 shows four
>  applications of the inverse of the velocity mass matrix. In the PETSc 
> implementation there are at
>  most two applications of the reciprocal of the diagonal of A (an 
> approximation to the velocity mass
>  matrix without more plumbing, as already pointed out). It seems like for 
> code implementations in
>  which there are possible scaling differences between the velocity and 
> pressure equations, that this
>  difference in the number of inverse applications could be significant? I 
> know Jed said that these
>  scalings wouldn't really matter if you have a uniform grid, but I'm not 100% 
> convinced yet.
>
>  I might try fiddling around with adding two more reciprocal applications.
>
>  On Fri, Jun 23, 2023 at 1:09 PM Pierre Jolivet  
> wrote:
>
>  On 23 Jun 2023, at 10:06 PM, Pierre Jolivet  wrote:
>
>  On 23 Jun 2023, at 9:39 PM, Alexander Lindsay  
> wrote:
>
>  Ah, I see that if I use Pierre's new 'full' option for 
> -mat_schur_complement_ainv_type
>
>  That was not initially done by me
>
>  Oops, sorry for the noise, looks like it was done by me indeed in
>  9399e4fd88c6621aad8fe9558ce84df37bd6fada…
>
>  Thanks,
>  Pierre
>
>  (though I recently tweaked MatSchurComplementComputeExplicitOperator() a bit 
> to use
>  KSPMatSolve(), so that if you have a small Schur complement — which is not 
> really the case
>  for NS — this could be a viable option, it was previously painfully slow).
>
>  Thanks,
>  Pierre
>
>  that I get a single iteration for the Schur complement solve with LU. That's 
> a nice testing
>  option
>
>  On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay 
>  wrote:
>
>  I guess it is because the inverse of the diagonal form of A00 becomes a poor
>  representation of the inverse of A00? I guess naively I would have thought 
> that the
>  blockdiag form of A00 is A00
>
>  On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay 
>  wrote:
>
>  Hi Jed, I will come back with answers to all of your questions at some 
> point. I
>  mostly just deal with MOOSE users who come to me and tell me their solve is
>  converging slowly, asking me how to fix it. So I generally assume they have
>  built an appropriate mesh and problem size for the problem they want to solve
>  and added appropriate turbulence modeling (although my general assumption
>  is often violated).
>
>  > And to confirm, are you doing a nonlinearly implicit velocity-pressure 
> solve?
>
>  Yes, this is our default.
>
>  A 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

Do you know of anyone who has applied the augmented Lagrange methodology to a finite volume discretization?On Jul 6, 2023, at 6:25 PM, Matthew Knepley  wrote:On Thu, Jul 6, 2023 at 8:30 PM Alexander Lindsay  wrote:This is an interesting article that compares a multi-level ILU algorithm to approximate commutator and augmented Lagrange methods: https://doi.org/10.1002/fld.5039That is for incompressible NS. The results are not better than https://arxiv.org/abs/1810.03315, and that PC is considerably simpler and already implemented in PETSc. There is an update in to this  https://epubs.siam.org/doi/abs/10.1137/21M1430698?casa_token=Fp_XhuZStZ0A:YDhnkW9XvAom_b8KocWz-hBEI7FAt46aw3ICa0FvCrOVCtYr9bwvtqJ4aBOTkDSvANKh6YTQEwwhich removes the need for complicated elements.You might need stuff like ILU for compressible flow, but I think incompressible is solved.  Thanks,     Matt On Wed, Jun 28, 2023 at 11:37 AM Alexander Lindsay  wrote:I do believe that based off the results in https://doi.org/10.1137/040608817 we should be able to make LSC, with proper scaling, compare very favorably with PCDOn Tue, Jun 27, 2023 at 10:41 AM Alexander Lindsay  wrote:I've opened https://gitlab.com/petsc/petsc/-/merge_requests/6642 which adds a couple more scaling applications of the inverse of the diagonal of AOn Mon, Jun 26, 2023 at 6:06 PM Alexander Lindsay  wrote:I guess that similar to the discussions about selfp, the approximation of the velocity mass matrix by the diagonal of the velocity sub-matrix will improve when running a transient as opposed to a steady calculation, especially if the time derivative is lumped Just thinking while typingOn Mon, Jun 26, 2023 at 6:03 PM Alexander Lindsay  wrote:Returning to Sebastian's question about the correctness of the current LSC implementation: in the taxonomy paper that Jed linked to (which talks about SIMPLE, PCD, and LSC), equation 21 shows four applications of the inverse of the velocity mass matrix. In the PETSc implementation there are at most two applications of the reciprocal of the diagonal of A (an approximation to the velocity mass matrix without more plumbing, as already pointed out). It seems like for code implementations in which there are possible scaling differences between the velocity and pressure equations, that this difference in the number of inverse applications could be significant? I know Jed said that these scalings wouldn't really matter if you have a uniform grid, but I'm not 100% convinced yet.I might try fiddling around with adding two more reciprocal applications.On Fri, Jun 23, 2023 at 1:09 PM Pierre Jolivet  wrote:On 23 Jun 2023, at 10:06 PM, Pierre Jolivet  wrote:On 23 Jun 2023, at 9:39 PM, Alexander Lindsay  wrote:Ah, I see that if I use Pierre's new 'full' option for -mat_schur_complement_ainv_typeThat was not initially done by meOops, sorry for the noise, looks like it was done by me indeed in 9399e4fd88c6621aad8fe9558ce84df37bd6fada…Thanks,Pierre (though I recently tweaked MatSchurComplementComputeExplicitOperator() a bit to use KSPMatSolve(), so that if you have a small Schur complement — which is not really the case for NS — this could be a viable option, it was previously painfully slow).Thanks,Pierre that I get a single iteration for the Schur complement solve with LU. That's a nice testing optionOn Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay  wrote:I guess it is because the inverse of the diagonal form of A00 becomes a poor representation of the inverse of A00? I guess naively I would have thought that the blockdiag form of A00 is A00On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay  wrote:Hi Jed, I will come back with answers to all of your questions at some point. I mostly just deal with MOOSE users who come to me and tell me their solve is converging slowly, asking me how to fix it. So I generally assume they have built an appropriate mesh and problem size for the problem they want to solve and added appropriate turbulence modeling (although my general assumption is often violated).> And to confirm, are you doing a nonlinearly implicit velocity-pressure solve?Yes, this is our default.A general question: it seems that it is well known that the quality of selfp degrades with increasing advection. Why is that?On Wed, Jun 7, 2023 at 8:01 PM Jed Brown  wrote:Alexander Lindsay  writes:

> This has been a great discussion to follow. Regarding
>
>> when time stepping, you have enough mass matrix that cheaper preconditioners are good enough
>
> I'm curious what some algebraic recommendations might be for high Re in
> transients. 

What mesh aspect ratio and streamline CFL number? Assuming your model is turbulent, can you say anything about 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Matthew Knepley]

On Thu, Jul 6, 2023 at 8:30 PM Alexander Lindsay 
wrote:

> This is an interesting article that compares a multi-level ILU algorithm
> to approximate commutator and augmented Lagrange methods:
> https://doi.org/10.1002/fld.5039
>

That is for incompressible NS. The results are not better than
https://arxiv.org/abs/1810.03315, and that PC is considerably simpler and
already implemented in PETSc. There is an update in to this


https://epubs.siam.org/doi/abs/10.1137/21M1430698?casa_token=Fp_XhuZStZ0A:YDhnkW9XvAom_b8KocWz-hBEI7FAt46aw3ICa0FvCrOVCtYr9bwvtqJ4aBOTkDSvANKh6YTQEw

which removes the need for complicated elements.

You might need stuff like ILU for compressible flow, but I think
incompressible is solved.

  Thanks,

 Matt


> On Wed, Jun 28, 2023 at 11:37 AM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
>
>> I do believe that based off the results in
>> https://doi.org/10.1137/040608817 we should be able to make LSC, with
>> proper scaling, compare very favorably with PCD
>>
>> On Tue, Jun 27, 2023 at 10:41 AM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>>
>>> I've opened https://gitlab.com/petsc/petsc/-/merge_requests/6642 which
>>> adds a couple more scaling applications of the inverse of the diagonal of A
>>>
>>> On Mon, Jun 26, 2023 at 6:06 PM Alexander Lindsay <
>>> alexlindsay...@gmail.com> wrote:
>>>
 I guess that similar to the discussions about selfp, the approximation
 of the velocity mass matrix by the diagonal of the velocity sub-matrix will
 improve when running a transient as opposed to a steady calculation,
 especially if the time derivative is lumped Just thinking while typing

 On Mon, Jun 26, 2023 at 6:03 PM Alexander Lindsay <
 alexlindsay...@gmail.com> wrote:

> Returning to Sebastian's question about the correctness of the current
> LSC implementation: in the taxonomy paper that Jed linked to (which talks
> about SIMPLE, PCD, and LSC), equation 21 shows four applications of the
> inverse of the velocity mass matrix. In the PETSc implementation there are
> at most two applications of the reciprocal of the diagonal of A (an
> approximation to the velocity mass matrix without more plumbing, as 
> already
> pointed out). It seems like for code implementations in which there are
> possible scaling differences between the velocity and pressure equations,
> that this difference in the number of inverse applications could be
> significant? I know Jed said that these scalings wouldn't really matter if
> you have a uniform grid, but I'm not 100% convinced yet.
>
> I might try fiddling around with adding two more reciprocal
> applications.
>
> On Fri, Jun 23, 2023 at 1:09 PM Pierre Jolivet 
> wrote:
>
>>
>> On 23 Jun 2023, at 10:06 PM, Pierre Jolivet 
>> wrote:
>>
>>
>> On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>>
>> Ah, I see that if I use Pierre's new 'full' option for
>> -mat_schur_complement_ainv_type
>>
>>
>> That was not initially done by me
>>
>>
>> Oops, sorry for the noise, looks like it was done by me indeed
>> in 9399e4fd88c6621aad8fe9558ce84df37bd6fada…
>>
>> Thanks,
>> Pierre
>>
>> (though I recently
>> tweaked MatSchurComplementComputeExplicitOperator() a bit to use
>> KSPMatSolve(), so that if you have a small Schur complement — which is 
>> not
>> really the case for NS — this could be a viable option, it was previously
>> painfully slow).
>>
>> Thanks,
>> Pierre
>>
>> that I get a single iteration for the Schur complement solve with LU.
>> That's a nice testing option
>>
>> On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>>
>>> I guess it is because the inverse of the diagonal form of A00
>>> becomes a poor representation of the inverse of A00? I guess naively I
>>> would have thought that the blockdiag form of A00 is A00
>>>
>>> On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay <
>>> alexlindsay...@gmail.com> wrote:
>>>
 Hi Jed, I will come back with answers to all of your questions at
 some point. I mostly just deal with MOOSE users who come to me and 
 tell me
 their solve is converging slowly, asking me how to fix it. So I 
 generally
 assume they have built an appropriate mesh and problem size for the 
 problem
 they want to solve and added appropriate turbulence modeling (although 
 my
 general assumption is often violated).

 > And to confirm, are you doing a nonlinearly implicit
 velocity-pressure solve?

 Yes, this is our default.

 A general question: it seems that it is well known that the quality
 of selfp degrades with 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

This is an interesting article that compares a multi-level ILU algorithm to
approximate commutator and augmented Lagrange methods:
https://doi.org/10.1002/fld.5039

On Wed, Jun 28, 2023 at 11:37 AM Alexander Lindsay 
wrote:

> I do believe that based off the results in
> https://doi.org/10.1137/040608817 we should be able to make LSC, with
> proper scaling, compare very favorably with PCD
>
> On Tue, Jun 27, 2023 at 10:41 AM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
>
>> I've opened https://gitlab.com/petsc/petsc/-/merge_requests/6642 which
>> adds a couple more scaling applications of the inverse of the diagonal of A
>>
>> On Mon, Jun 26, 2023 at 6:06 PM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>>
>>> I guess that similar to the discussions about selfp, the approximation
>>> of the velocity mass matrix by the diagonal of the velocity sub-matrix will
>>> improve when running a transient as opposed to a steady calculation,
>>> especially if the time derivative is lumped Just thinking while typing
>>>
>>> On Mon, Jun 26, 2023 at 6:03 PM Alexander Lindsay <
>>> alexlindsay...@gmail.com> wrote:
>>>
 Returning to Sebastian's question about the correctness of the current
 LSC implementation: in the taxonomy paper that Jed linked to (which talks
 about SIMPLE, PCD, and LSC), equation 21 shows four applications of the
 inverse of the velocity mass matrix. In the PETSc implementation there are
 at most two applications of the reciprocal of the diagonal of A (an
 approximation to the velocity mass matrix without more plumbing, as already
 pointed out). It seems like for code implementations in which there are
 possible scaling differences between the velocity and pressure equations,
 that this difference in the number of inverse applications could be
 significant? I know Jed said that these scalings wouldn't really matter if
 you have a uniform grid, but I'm not 100% convinced yet.

 I might try fiddling around with adding two more reciprocal
 applications.

 On Fri, Jun 23, 2023 at 1:09 PM Pierre Jolivet 
 wrote:

>
> On 23 Jun 2023, at 10:06 PM, Pierre Jolivet 
> wrote:
>
>
> On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
>
> Ah, I see that if I use Pierre's new 'full' option for
> -mat_schur_complement_ainv_type
>
>
> That was not initially done by me
>
>
> Oops, sorry for the noise, looks like it was done by me indeed
> in 9399e4fd88c6621aad8fe9558ce84df37bd6fada…
>
> Thanks,
> Pierre
>
> (though I recently tweaked MatSchurComplementComputeExplicitOperator()
> a bit to use KSPMatSolve(), so that if you have a small Schur complement —
> which is not really the case for NS — this could be a viable option, it 
> was
> previously painfully slow).
>
> Thanks,
> Pierre
>
> that I get a single iteration for the Schur complement solve with LU.
> That's a nice testing option
>
> On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
>
>> I guess it is because the inverse of the diagonal form of A00 becomes
>> a poor representation of the inverse of A00? I guess naively I would have
>> thought that the blockdiag form of A00 is A00
>>
>> On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>>
>>> Hi Jed, I will come back with answers to all of your questions at
>>> some point. I mostly just deal with MOOSE users who come to me and tell 
>>> me
>>> their solve is converging slowly, asking me how to fix it. So I 
>>> generally
>>> assume they have built an appropriate mesh and problem size for the 
>>> problem
>>> they want to solve and added appropriate turbulence modeling (although 
>>> my
>>> general assumption is often violated).
>>>
>>> > And to confirm, are you doing a nonlinearly implicit
>>> velocity-pressure solve?
>>>
>>> Yes, this is our default.
>>>
>>> A general question: it seems that it is well known that the quality
>>> of selfp degrades with increasing advection. Why is that?
>>>
>>> On Wed, Jun 7, 2023 at 8:01 PM Jed Brown  wrote:
>>>
 Alexander Lindsay  writes:

 > This has been a great discussion to follow. Regarding
 >
 >> when time stepping, you have enough mass matrix that cheaper
 preconditioners are good enough
 >
 > I'm curious what some algebraic recommendations might be for high
 Re in
 > transients.

 What mesh aspect ratio and streamline CFL number? Assuming your
 model is turbulent, can you say anything about momentum thickness 
 Reynolds
 number Re_θ? What is your wall normal spacing in plus units? (Wall 
 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

I do believe that based off the results in https://doi.org/10.1137/040608817 we
should be able to make LSC, with proper scaling, compare very favorably
with PCD

On Tue, Jun 27, 2023 at 10:41 AM Alexander Lindsay 
wrote:

> I've opened https://gitlab.com/petsc/petsc/-/merge_requests/6642 which
> adds a couple more scaling applications of the inverse of the diagonal of A
>
> On Mon, Jun 26, 2023 at 6:06 PM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
>
>> I guess that similar to the discussions about selfp, the approximation of
>> the velocity mass matrix by the diagonal of the velocity sub-matrix will
>> improve when running a transient as opposed to a steady calculation,
>> especially if the time derivative is lumped Just thinking while typing
>>
>> On Mon, Jun 26, 2023 at 6:03 PM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>>
>>> Returning to Sebastian's question about the correctness of the current
>>> LSC implementation: in the taxonomy paper that Jed linked to (which talks
>>> about SIMPLE, PCD, and LSC), equation 21 shows four applications of the
>>> inverse of the velocity mass matrix. In the PETSc implementation there are
>>> at most two applications of the reciprocal of the diagonal of A (an
>>> approximation to the velocity mass matrix without more plumbing, as already
>>> pointed out). It seems like for code implementations in which there are
>>> possible scaling differences between the velocity and pressure equations,
>>> that this difference in the number of inverse applications could be
>>> significant? I know Jed said that these scalings wouldn't really matter if
>>> you have a uniform grid, but I'm not 100% convinced yet.
>>>
>>> I might try fiddling around with adding two more reciprocal applications.
>>>
>>> On Fri, Jun 23, 2023 at 1:09 PM Pierre Jolivet 
>>> wrote:
>>>

 On 23 Jun 2023, at 10:06 PM, Pierre Jolivet 
 wrote:


 On 23 Jun 2023, at 9:39 PM, Alexander Lindsay 
 wrote:

 Ah, I see that if I use Pierre's new 'full' option for
 -mat_schur_complement_ainv_type


 That was not initially done by me


 Oops, sorry for the noise, looks like it was done by me indeed
 in 9399e4fd88c6621aad8fe9558ce84df37bd6fada…

 Thanks,
 Pierre

 (though I recently tweaked MatSchurComplementComputeExplicitOperator()
 a bit to use KSPMatSolve(), so that if you have a small Schur complement —
 which is not really the case for NS — this could be a viable option, it was
 previously painfully slow).

 Thanks,
 Pierre

 that I get a single iteration for the Schur complement solve with LU.
 That's a nice testing option

 On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay <
 alexlindsay...@gmail.com> wrote:

> I guess it is because the inverse of the diagonal form of A00 becomes
> a poor representation of the inverse of A00? I guess naively I would have
> thought that the blockdiag form of A00 is A00
>
> On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
>
>> Hi Jed, I will come back with answers to all of your questions at
>> some point. I mostly just deal with MOOSE users who come to me and tell 
>> me
>> their solve is converging slowly, asking me how to fix it. So I generally
>> assume they have built an appropriate mesh and problem size for the 
>> problem
>> they want to solve and added appropriate turbulence modeling (although my
>> general assumption is often violated).
>>
>> > And to confirm, are you doing a nonlinearly implicit
>> velocity-pressure solve?
>>
>> Yes, this is our default.
>>
>> A general question: it seems that it is well known that the quality
>> of selfp degrades with increasing advection. Why is that?
>>
>> On Wed, Jun 7, 2023 at 8:01 PM Jed Brown  wrote:
>>
>>> Alexander Lindsay  writes:
>>>
>>> > This has been a great discussion to follow. Regarding
>>> >
>>> >> when time stepping, you have enough mass matrix that cheaper
>>> preconditioners are good enough
>>> >
>>> > I'm curious what some algebraic recommendations might be for high
>>> Re in
>>> > transients.
>>>
>>> What mesh aspect ratio and streamline CFL number? Assuming your
>>> model is turbulent, can you say anything about momentum thickness 
>>> Reynolds
>>> number Re_θ? What is your wall normal spacing in plus units? (Wall 
>>> resolved
>>> or wall modeled?)
>>>
>>> And to confirm, are you doing a nonlinearly implicit
>>> velocity-pressure solve?
>>>
>>> > I've found one-level DD to be ineffective when applied
>>> monolithically or to the momentum block of a split, as it scales with 
>>> the
>>> mesh size.
>>>
>>> I wouldn't put too much weight on "scaling with mesh size" per se.
>>> You want an efficient 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

I've opened https://gitlab.com/petsc/petsc/-/merge_requests/6642 which adds
a couple more scaling applications of the inverse of the diagonal of A

On Mon, Jun 26, 2023 at 6:06 PM Alexander Lindsay 
wrote:

> I guess that similar to the discussions about selfp, the approximation of
> the velocity mass matrix by the diagonal of the velocity sub-matrix will
> improve when running a transient as opposed to a steady calculation,
> especially if the time derivative is lumped Just thinking while typing
>
> On Mon, Jun 26, 2023 at 6:03 PM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
>
>> Returning to Sebastian's question about the correctness of the current
>> LSC implementation: in the taxonomy paper that Jed linked to (which talks
>> about SIMPLE, PCD, and LSC), equation 21 shows four applications of the
>> inverse of the velocity mass matrix. In the PETSc implementation there are
>> at most two applications of the reciprocal of the diagonal of A (an
>> approximation to the velocity mass matrix without more plumbing, as already
>> pointed out). It seems like for code implementations in which there are
>> possible scaling differences between the velocity and pressure equations,
>> that this difference in the number of inverse applications could be
>> significant? I know Jed said that these scalings wouldn't really matter if
>> you have a uniform grid, but I'm not 100% convinced yet.
>>
>> I might try fiddling around with adding two more reciprocal applications.
>>
>> On Fri, Jun 23, 2023 at 1:09 PM Pierre Jolivet 
>> wrote:
>>
>>>
>>> On 23 Jun 2023, at 10:06 PM, Pierre Jolivet 
>>> wrote:
>>>
>>>
>>> On 23 Jun 2023, at 9:39 PM, Alexander Lindsay 
>>> wrote:
>>>
>>> Ah, I see that if I use Pierre's new 'full' option for
>>> -mat_schur_complement_ainv_type
>>>
>>>
>>> That was not initially done by me
>>>
>>>
>>> Oops, sorry for the noise, looks like it was done by me indeed
>>> in 9399e4fd88c6621aad8fe9558ce84df37bd6fada…
>>>
>>> Thanks,
>>> Pierre
>>>
>>> (though I recently tweaked MatSchurComplementComputeExplicitOperator() a
>>> bit to use KSPMatSolve(), so that if you have a small Schur complement —
>>> which is not really the case for NS — this could be a viable option, it was
>>> previously painfully slow).
>>>
>>> Thanks,
>>> Pierre
>>>
>>> that I get a single iteration for the Schur complement solve with LU.
>>> That's a nice testing option
>>>
>>> On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay <
>>> alexlindsay...@gmail.com> wrote:
>>>
 I guess it is because the inverse of the diagonal form of A00 becomes a
 poor representation of the inverse of A00? I guess naively I would have
 thought that the blockdiag form of A00 is A00

 On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay <
 alexlindsay...@gmail.com> wrote:

> Hi Jed, I will come back with answers to all of your questions at some
> point. I mostly just deal with MOOSE users who come to me and tell me 
> their
> solve is converging slowly, asking me how to fix it. So I generally assume
> they have built an appropriate mesh and problem size for the problem they
> want to solve and added appropriate turbulence modeling (although my
> general assumption is often violated).
>
> > And to confirm, are you doing a nonlinearly implicit
> velocity-pressure solve?
>
> Yes, this is our default.
>
> A general question: it seems that it is well known that the quality of
> selfp degrades with increasing advection. Why is that?
>
> On Wed, Jun 7, 2023 at 8:01 PM Jed Brown  wrote:
>
>> Alexander Lindsay  writes:
>>
>> > This has been a great discussion to follow. Regarding
>> >
>> >> when time stepping, you have enough mass matrix that cheaper
>> preconditioners are good enough
>> >
>> > I'm curious what some algebraic recommendations might be for high
>> Re in
>> > transients.
>>
>> What mesh aspect ratio and streamline CFL number? Assuming your model
>> is turbulent, can you say anything about momentum thickness Reynolds 
>> number
>> Re_θ? What is your wall normal spacing in plus units? (Wall resolved or
>> wall modeled?)
>>
>> And to confirm, are you doing a nonlinearly implicit
>> velocity-pressure solve?
>>
>> > I've found one-level DD to be ineffective when applied
>> monolithically or to the momentum block of a split, as it scales with the
>> mesh size.
>>
>> I wouldn't put too much weight on "scaling with mesh size" per se.
>> You want an efficient solver for the coarsest mesh that delivers 
>> sufficient
>> accuracy in your flow regime. Constants matter.
>>
>> Refining the mesh while holding time steps constant changes the
>> advective CFL number as well as cell Peclet/cell Reynolds numbers. A
>> meaningful scaling study is to increase Reynolds number (e.g., by growing
>> the domain) while keeping mesh size 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

I guess that similar to the discussions about selfp, the approximation of
the velocity mass matrix by the diagonal of the velocity sub-matrix will
improve when running a transient as opposed to a steady calculation,
especially if the time derivative is lumped Just thinking while typing

On Mon, Jun 26, 2023 at 6:03 PM Alexander Lindsay 
wrote:

> Returning to Sebastian's question about the correctness of the current LSC
> implementation: in the taxonomy paper that Jed linked to (which talks about
> SIMPLE, PCD, and LSC), equation 21 shows four applications of the inverse
> of the velocity mass matrix. In the PETSc implementation there are at most
> two applications of the reciprocal of the diagonal of A (an approximation
> to the velocity mass matrix without more plumbing, as already pointed out).
> It seems like for code implementations in which there are possible scaling
> differences between the velocity and pressure equations, that this
> difference in the number of inverse applications could be significant? I
> know Jed said that these scalings wouldn't really matter if you have a
> uniform grid, but I'm not 100% convinced yet.
>
> I might try fiddling around with adding two more reciprocal applications.
>
> On Fri, Jun 23, 2023 at 1:09 PM Pierre Jolivet 
> wrote:
>
>>
>> On 23 Jun 2023, at 10:06 PM, Pierre Jolivet 
>> wrote:
>>
>>
>> On 23 Jun 2023, at 9:39 PM, Alexander Lindsay 
>> wrote:
>>
>> Ah, I see that if I use Pierre's new 'full' option for
>> -mat_schur_complement_ainv_type
>>
>>
>> That was not initially done by me
>>
>>
>> Oops, sorry for the noise, looks like it was done by me indeed
>> in 9399e4fd88c6621aad8fe9558ce84df37bd6fada…
>>
>> Thanks,
>> Pierre
>>
>> (though I recently tweaked MatSchurComplementComputeExplicitOperator() a
>> bit to use KSPMatSolve(), so that if you have a small Schur complement —
>> which is not really the case for NS — this could be a viable option, it was
>> previously painfully slow).
>>
>> Thanks,
>> Pierre
>>
>> that I get a single iteration for the Schur complement solve with LU.
>> That's a nice testing option
>>
>> On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>>
>>> I guess it is because the inverse of the diagonal form of A00 becomes a
>>> poor representation of the inverse of A00? I guess naively I would have
>>> thought that the blockdiag form of A00 is A00
>>>
>>> On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay <
>>> alexlindsay...@gmail.com> wrote:
>>>
 Hi Jed, I will come back with answers to all of your questions at some
 point. I mostly just deal with MOOSE users who come to me and tell me their
 solve is converging slowly, asking me how to fix it. So I generally assume
 they have built an appropriate mesh and problem size for the problem they
 want to solve and added appropriate turbulence modeling (although my
 general assumption is often violated).

 > And to confirm, are you doing a nonlinearly implicit
 velocity-pressure solve?

 Yes, this is our default.

 A general question: it seems that it is well known that the quality of
 selfp degrades with increasing advection. Why is that?

 On Wed, Jun 7, 2023 at 8:01 PM Jed Brown  wrote:

> Alexander Lindsay  writes:
>
> > This has been a great discussion to follow. Regarding
> >
> >> when time stepping, you have enough mass matrix that cheaper
> preconditioners are good enough
> >
> > I'm curious what some algebraic recommendations might be for high Re
> in
> > transients.
>
> What mesh aspect ratio and streamline CFL number? Assuming your model
> is turbulent, can you say anything about momentum thickness Reynolds 
> number
> Re_θ? What is your wall normal spacing in plus units? (Wall resolved or
> wall modeled?)
>
> And to confirm, are you doing a nonlinearly implicit velocity-pressure
> solve?
>
> > I've found one-level DD to be ineffective when applied
> monolithically or to the momentum block of a split, as it scales with the
> mesh size.
>
> I wouldn't put too much weight on "scaling with mesh size" per se. You
> want an efficient solver for the coarsest mesh that delivers sufficient
> accuracy in your flow regime. Constants matter.
>
> Refining the mesh while holding time steps constant changes the
> advective CFL number as well as cell Peclet/cell Reynolds numbers. A
> meaningful scaling study is to increase Reynolds number (e.g., by growing
> the domain) while keeping mesh size matched in terms of plus units in the
> viscous sublayer and Kolmogorov length in the outer boundary layer. That
> turns out to not be a very automatic study to do, but it's what matters 
> and
> you can spend a lot of time chasing ghosts with naive scaling studies.
>

>>
>>

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

Returning to Sebastian's question about the correctness of the current LSC
implementation: in the taxonomy paper that Jed linked to (which talks about
SIMPLE, PCD, and LSC), equation 21 shows four applications of the inverse
of the velocity mass matrix. In the PETSc implementation there are at most
two applications of the reciprocal of the diagonal of A (an approximation
to the velocity mass matrix without more plumbing, as already pointed out).
It seems like for code implementations in which there are possible scaling
differences between the velocity and pressure equations, that this
difference in the number of inverse applications could be significant? I
know Jed said that these scalings wouldn't really matter if you have a
uniform grid, but I'm not 100% convinced yet.

I might try fiddling around with adding two more reciprocal applications.

On Fri, Jun 23, 2023 at 1:09 PM Pierre Jolivet 
wrote:

>
> On 23 Jun 2023, at 10:06 PM, Pierre Jolivet 
> wrote:
>
>
> On 23 Jun 2023, at 9:39 PM, Alexander Lindsay 
> wrote:
>
> Ah, I see that if I use Pierre's new 'full' option for
> -mat_schur_complement_ainv_type
>
>
> That was not initially done by me
>
>
> Oops, sorry for the noise, looks like it was done by me indeed
> in 9399e4fd88c6621aad8fe9558ce84df37bd6fada…
>
> Thanks,
> Pierre
>
> (though I recently tweaked MatSchurComplementComputeExplicitOperator() a
> bit to use KSPMatSolve(), so that if you have a small Schur complement —
> which is not really the case for NS — this could be a viable option, it was
> previously painfully slow).
>
> Thanks,
> Pierre
>
> that I get a single iteration for the Schur complement solve with LU.
> That's a nice testing option
>
> On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
>
>> I guess it is because the inverse of the diagonal form of A00 becomes a
>> poor representation of the inverse of A00? I guess naively I would have
>> thought that the blockdiag form of A00 is A00
>>
>> On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay <
>> alexlindsay...@gmail.com> wrote:
>>
>>> Hi Jed, I will come back with answers to all of your questions at some
>>> point. I mostly just deal with MOOSE users who come to me and tell me their
>>> solve is converging slowly, asking me how to fix it. So I generally assume
>>> they have built an appropriate mesh and problem size for the problem they
>>> want to solve and added appropriate turbulence modeling (although my
>>> general assumption is often violated).
>>>
>>> > And to confirm, are you doing a nonlinearly implicit velocity-pressure
>>> solve?
>>>
>>> Yes, this is our default.
>>>
>>> A general question: it seems that it is well known that the quality of
>>> selfp degrades with increasing advection. Why is that?
>>>
>>> On Wed, Jun 7, 2023 at 8:01 PM Jed Brown  wrote:
>>>
 Alexander Lindsay  writes:

 > This has been a great discussion to follow. Regarding
 >
 >> when time stepping, you have enough mass matrix that cheaper
 preconditioners are good enough
 >
 > I'm curious what some algebraic recommendations might be for high Re
 in
 > transients.

 What mesh aspect ratio and streamline CFL number? Assuming your model
 is turbulent, can you say anything about momentum thickness Reynolds number
 Re_θ? What is your wall normal spacing in plus units? (Wall resolved or
 wall modeled?)

 And to confirm, are you doing a nonlinearly implicit velocity-pressure
 solve?

 > I've found one-level DD to be ineffective when applied monolithically
 or to the momentum block of a split, as it scales with the mesh size.

 I wouldn't put too much weight on "scaling with mesh size" per se. You
 want an efficient solver for the coarsest mesh that delivers sufficient
 accuracy in your flow regime. Constants matter.

 Refining the mesh while holding time steps constant changes the
 advective CFL number as well as cell Peclet/cell Reynolds numbers. A
 meaningful scaling study is to increase Reynolds number (e.g., by growing
 the domain) while keeping mesh size matched in terms of plus units in the
 viscous sublayer and Kolmogorov length in the outer boundary layer. That
 turns out to not be a very automatic study to do, but it's what matters and
 you can spend a lot of time chasing ghosts with naive scaling studies.

>>>
>
>

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Pierre Jolivet]

> On 23 Jun 2023, at 10:06 PM, Pierre Jolivet  wrote:
> 
> 
>> On 23 Jun 2023, at 9:39 PM, Alexander Lindsay  
>> wrote:
>> 
>> Ah, I see that if I use Pierre's new 'full' option for 
>> -mat_schur_complement_ainv_type
> 
> That was not initially done by me

Oops, sorry for the noise, looks like it was done by me indeed in 
9399e4fd88c6621aad8fe9558ce84df37bd6fada…

Thanks,
Pierre

> (though I recently tweaked MatSchurComplementComputeExplicitOperator() a bit 
> to use KSPMatSolve(), so that if you have a small Schur complement — which is 
> not really the case for NS — this could be a viable option, it was previously 
> painfully slow).
> 
> Thanks,
> Pierre
> 
>> that I get a single iteration for the Schur complement solve with LU. That's 
>> a nice testing option
>> 
>> On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay > > wrote:
>>> I guess it is because the inverse of the diagonal form of A00 becomes a 
>>> poor representation of the inverse of A00? I guess naively I would have 
>>> thought that the blockdiag form of A00 is A00
>>> 
>>> On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay 
>>> mailto:alexlindsay...@gmail.com>> wrote:
 Hi Jed, I will come back with answers to all of your questions at some 
 point. I mostly just deal with MOOSE users who come to me and tell me 
 their solve is converging slowly, asking me how to fix it. So I generally 
 assume they have built an appropriate mesh and problem size for the 
 problem they want to solve and added appropriate turbulence modeling 
 (although my general assumption is often violated).
 
 > And to confirm, are you doing a nonlinearly implicit velocity-pressure 
 > solve?
 
 Yes, this is our default.
 
 A general question: it seems that it is well known that the quality of 
 selfp degrades with increasing advection. Why is that?
 
 On Wed, Jun 7, 2023 at 8:01 PM Jed Brown >>> > wrote:
> Alexander Lindsay  > writes:
> 
> > This has been a great discussion to follow. Regarding
> >
> >> when time stepping, you have enough mass matrix that cheaper 
> >> preconditioners are good enough
> >
> > I'm curious what some algebraic recommendations might be for high Re in
> > transients. 
> 
> What mesh aspect ratio and streamline CFL number? Assuming your model is 
> turbulent, can you say anything about momentum thickness Reynolds number 
> Re_θ? What is your wall normal spacing in plus units? (Wall resolved or 
> wall modeled?)
> 
> And to confirm, are you doing a nonlinearly implicit velocity-pressure 
> solve?
> 
> > I've found one-level DD to be ineffective when applied monolithically 
> > or to the momentum block of a split, as it scales with the mesh size. 
> 
> I wouldn't put too much weight on "scaling with mesh size" per se. You 
> want an efficient solver for the coarsest mesh that delivers sufficient 
> accuracy in your flow regime. Constants matter.
> 
> Refining the mesh while holding time steps constant changes the advective 
> CFL number as well as cell Peclet/cell Reynolds numbers. A meaningful 
> scaling study is to increase Reynolds number (e.g., by growing the 
> domain) while keeping mesh size matched in terms of plus units in the 
> viscous sublayer and Kolmogorov length in the outer boundary layer. That 
> turns out to not be a very automatic study to do, but it's what matters 
> and you can spend a lot of time chasing ghosts with naive scaling studies.
> 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Pierre Jolivet]

> On 23 Jun 2023, at 9:39 PM, Alexander Lindsay  
> wrote:
> 
> Ah, I see that if I use Pierre's new 'full' option for 
> -mat_schur_complement_ainv_type

That was not initially done by me (though I recently tweaked 
MatSchurComplementComputeExplicitOperator() a bit to use KSPMatSolve(), so that 
if you have a small Schur complement — which is not really the case for NS — 
this could be a viable option, it was previously painfully slow).

Thanks,
Pierre

> that I get a single iteration for the Schur complement solve with LU. That's 
> a nice testing option
> 
> On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay  > wrote:
>> I guess it is because the inverse of the diagonal form of A00 becomes a poor 
>> representation of the inverse of A00? I guess naively I would have thought 
>> that the blockdiag form of A00 is A00
>> 
>> On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay > > wrote:
>>> Hi Jed, I will come back with answers to all of your questions at some 
>>> point. I mostly just deal with MOOSE users who come to me and tell me their 
>>> solve is converging slowly, asking me how to fix it. So I generally assume 
>>> they have built an appropriate mesh and problem size for the problem they 
>>> want to solve and added appropriate turbulence modeling (although my 
>>> general assumption is often violated).
>>> 
>>> > And to confirm, are you doing a nonlinearly implicit velocity-pressure 
>>> > solve?
>>> 
>>> Yes, this is our default.
>>> 
>>> A general question: it seems that it is well known that the quality of 
>>> selfp degrades with increasing advection. Why is that?
>>> 
>>> On Wed, Jun 7, 2023 at 8:01 PM Jed Brown >> > wrote:
 Alexander Lindsay >>> > writes:
 
 > This has been a great discussion to follow. Regarding
 >
 >> when time stepping, you have enough mass matrix that cheaper 
 >> preconditioners are good enough
 >
 > I'm curious what some algebraic recommendations might be for high Re in
 > transients. 
 
 What mesh aspect ratio and streamline CFL number? Assuming your model is 
 turbulent, can you say anything about momentum thickness Reynolds number 
 Re_θ? What is your wall normal spacing in plus units? (Wall resolved or 
 wall modeled?)
 
 And to confirm, are you doing a nonlinearly implicit velocity-pressure 
 solve?
 
 > I've found one-level DD to be ineffective when applied monolithically or 
 > to the momentum block of a split, as it scales with the mesh size. 
 
 I wouldn't put too much weight on "scaling with mesh size" per se. You 
 want an efficient solver for the coarsest mesh that delivers sufficient 
 accuracy in your flow regime. Constants matter.
 
 Refining the mesh while holding time steps constant changes the advective 
 CFL number as well as cell Peclet/cell Reynolds numbers. A meaningful 
 scaling study is to increase Reynolds number (e.g., by growing the domain) 
 while keeping mesh size matched in terms of plus units in the viscous 
 sublayer and Kolmogorov length in the outer boundary layer. That turns out 
 to not be a very automatic study to do, but it's what matters and you can 
 spend a lot of time chasing ghosts with naive scaling studies.

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

Ah, I see that if I use Pierre's new 'full' option for
-mat_schur_complement_ainv_type that I get a single iteration for the Schur
complement solve with LU. That's a nice testing option

On Fri, Jun 23, 2023 at 12:02 PM Alexander Lindsay 
wrote:

> I guess it is because the inverse of the diagonal form of A00 becomes a
> poor representation of the inverse of A00? I guess naively I would have
> thought that the blockdiag form of A00 is A00
>
> On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay <
> alexlindsay...@gmail.com> wrote:
>
>> Hi Jed, I will come back with answers to all of your questions at some
>> point. I mostly just deal with MOOSE users who come to me and tell me their
>> solve is converging slowly, asking me how to fix it. So I generally assume
>> they have built an appropriate mesh and problem size for the problem they
>> want to solve and added appropriate turbulence modeling (although my
>> general assumption is often violated).
>>
>> > And to confirm, are you doing a nonlinearly implicit velocity-pressure
>> solve?
>>
>> Yes, this is our default.
>>
>> A general question: it seems that it is well known that the quality of
>> selfp degrades with increasing advection. Why is that?
>>
>> On Wed, Jun 7, 2023 at 8:01 PM Jed Brown  wrote:
>>
>>> Alexander Lindsay  writes:
>>>
>>> > This has been a great discussion to follow. Regarding
>>> >
>>> >> when time stepping, you have enough mass matrix that cheaper
>>> preconditioners are good enough
>>> >
>>> > I'm curious what some algebraic recommendations might be for high Re in
>>> > transients.
>>>
>>> What mesh aspect ratio and streamline CFL number? Assuming your model is
>>> turbulent, can you say anything about momentum thickness Reynolds number
>>> Re_θ? What is your wall normal spacing in plus units? (Wall resolved or
>>> wall modeled?)
>>>
>>> And to confirm, are you doing a nonlinearly implicit velocity-pressure
>>> solve?
>>>
>>> > I've found one-level DD to be ineffective when applied monolithically
>>> or to the momentum block of a split, as it scales with the mesh size.
>>>
>>> I wouldn't put too much weight on "scaling with mesh size" per se. You
>>> want an efficient solver for the coarsest mesh that delivers sufficient
>>> accuracy in your flow regime. Constants matter.
>>>
>>> Refining the mesh while holding time steps constant changes the
>>> advective CFL number as well as cell Peclet/cell Reynolds numbers. A
>>> meaningful scaling study is to increase Reynolds number (e.g., by growing
>>> the domain) while keeping mesh size matched in terms of plus units in the
>>> viscous sublayer and Kolmogorov length in the outer boundary layer. That
>>> turns out to not be a very automatic study to do, but it's what matters and
>>> you can spend a lot of time chasing ghosts with naive scaling studies.
>>>
>>

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

I guess it is because the inverse of the diagonal form of A00 becomes a
poor representation of the inverse of A00? I guess naively I would have
thought that the blockdiag form of A00 is A00

On Fri, Jun 23, 2023 at 10:18 AM Alexander Lindsay 
wrote:

> Hi Jed, I will come back with answers to all of your questions at some
> point. I mostly just deal with MOOSE users who come to me and tell me their
> solve is converging slowly, asking me how to fix it. So I generally assume
> they have built an appropriate mesh and problem size for the problem they
> want to solve and added appropriate turbulence modeling (although my
> general assumption is often violated).
>
> > And to confirm, are you doing a nonlinearly implicit velocity-pressure
> solve?
>
> Yes, this is our default.
>
> A general question: it seems that it is well known that the quality of
> selfp degrades with increasing advection. Why is that?
>
> On Wed, Jun 7, 2023 at 8:01 PM Jed Brown  wrote:
>
>> Alexander Lindsay  writes:
>>
>> > This has been a great discussion to follow. Regarding
>> >
>> >> when time stepping, you have enough mass matrix that cheaper
>> preconditioners are good enough
>> >
>> > I'm curious what some algebraic recommendations might be for high Re in
>> > transients.
>>
>> What mesh aspect ratio and streamline CFL number? Assuming your model is
>> turbulent, can you say anything about momentum thickness Reynolds number
>> Re_θ? What is your wall normal spacing in plus units? (Wall resolved or
>> wall modeled?)
>>
>> And to confirm, are you doing a nonlinearly implicit velocity-pressure
>> solve?
>>
>> > I've found one-level DD to be ineffective when applied monolithically
>> or to the momentum block of a split, as it scales with the mesh size.
>>
>> I wouldn't put too much weight on "scaling with mesh size" per se. You
>> want an efficient solver for the coarsest mesh that delivers sufficient
>> accuracy in your flow regime. Constants matter.
>>
>> Refining the mesh while holding time steps constant changes the advective
>> CFL number as well as cell Peclet/cell Reynolds numbers. A meaningful
>> scaling study is to increase Reynolds number (e.g., by growing the domain)
>> while keeping mesh size matched in terms of plus units in the viscous
>> sublayer and Kolmogorov length in the outer boundary layer. That turns out
>> to not be a very automatic study to do, but it's what matters and you can
>> spend a lot of time chasing ghosts with naive scaling studies.
>>
>

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

Hi Jed, I will come back with answers to all of your questions at some
point. I mostly just deal with MOOSE users who come to me and tell me their
solve is converging slowly, asking me how to fix it. So I generally assume
they have built an appropriate mesh and problem size for the problem they
want to solve and added appropriate turbulence modeling (although my
general assumption is often violated).

> And to confirm, are you doing a nonlinearly implicit velocity-pressure
solve?

Yes, this is our default.

A general question: it seems that it is well known that the quality of
selfp degrades with increasing advection. Why is that?

On Wed, Jun 7, 2023 at 8:01 PM Jed Brown  wrote:

> Alexander Lindsay  writes:
>
> > This has been a great discussion to follow. Regarding
> >
> >> when time stepping, you have enough mass matrix that cheaper
> preconditioners are good enough
> >
> > I'm curious what some algebraic recommendations might be for high Re in
> > transients.
>
> What mesh aspect ratio and streamline CFL number? Assuming your model is
> turbulent, can you say anything about momentum thickness Reynolds number
> Re_θ? What is your wall normal spacing in plus units? (Wall resolved or
> wall modeled?)
>
> And to confirm, are you doing a nonlinearly implicit velocity-pressure
> solve?
>
> > I've found one-level DD to be ineffective when applied monolithically or
> to the momentum block of a split, as it scales with the mesh size.
>
> I wouldn't put too much weight on "scaling with mesh size" per se. You
> want an efficient solver for the coarsest mesh that delivers sufficient
> accuracy in your flow regime. Constants matter.
>
> Refining the mesh while holding time steps constant changes the advective
> CFL number as well as cell Peclet/cell Reynolds numbers. A meaningful
> scaling study is to increase Reynolds number (e.g., by growing the domain)
> while keeping mesh size matched in terms of plus units in the viscous
> sublayer and Kolmogorov length in the outer boundary layer. That turns out
> to not be a very automatic study to do, but it's what matters and you can
> spend a lot of time chasing ghosts with naive scaling studies.
>

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Jed Brown]

Alexander Lindsay  writes:

> This has been a great discussion to follow. Regarding
>
>> when time stepping, you have enough mass matrix that cheaper preconditioners 
>> are good enough
>
> I'm curious what some algebraic recommendations might be for high Re in
> transients. 

What mesh aspect ratio and streamline CFL number? Assuming your model is 
turbulent, can you say anything about momentum thickness Reynolds number Re_θ? 
What is your wall normal spacing in plus units? (Wall resolved or wall modeled?)

And to confirm, are you doing a nonlinearly implicit velocity-pressure solve?

> I've found one-level DD to be ineffective when applied monolithically or to 
> the momentum block of a split, as it scales with the mesh size. 

I wouldn't put too much weight on "scaling with mesh size" per se. You want an 
efficient solver for the coarsest mesh that delivers sufficient accuracy in 
your flow regime. Constants matter.

Refining the mesh while holding time steps constant changes the advective CFL 
number as well as cell Peclet/cell Reynolds numbers. A meaningful scaling study 
is to increase Reynolds number (e.g., by growing the domain) while keeping mesh 
size matched in terms of plus units in the viscous sublayer and Kolmogorov 
length in the outer boundary layer. That turns out to not be a very automatic 
study to do, but it's what matters and you can spend a lot of time chasing 
ghosts with naive scaling studies.

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Alexander Lindsay]

This has been a great discussion to follow. Regarding

> when time stepping, you have enough mass matrix that cheaper
preconditioners are good enough

I'm curious what some algebraic recommendations might be for high Re in
transients. I've found one-level DD to be ineffective when applied
monolithically or to the momentum block of a split, as it scales with the
mesh size. For high Re boomeramg is ineffective perhaps until
https://gitlab.com/petsc/petsc/-/issues/1362 is resolved. I should try
fiddling around again with Pierre's work in HPDDM, but curious if there are
other PETSc PC recs, or if I need to overcome my inertia/laziness and move
beyond command line options.

On Mon, May 8, 2023 at 6:46 AM Jed Brown  wrote:

> Sebastian Blauth  writes:
>
> > Hello everyone,
> >
> > I wanted to briefly follow up on my question (see my last reply).
> > Does anyone know / have an idea why the LSC preconditioner in PETSc does
> > not seem to scale well with the problem size (the outer fgmres solver I
> > am using nearly scale nearly linearly with the problem size in my
> example).
>
> The implementation was tested on heterogeneous Stokes problems from
> geodynamics, and perhaps not on NS (or not with the discretization you're
> using).
>
> https://doi.org/10.1016/j.pepi.2008.07.036
>
> There is a comment about not having plumbing to provide a mass matrix. A
> few lines earlier there is code using PetscObjectQuery, and that same
> pattern could be applied for the mass matrix. If you're on a roughly
> uniform mesh, including the mass scaling will probably have little effect,
> but it could have a big impact in the presence of highly anistropic
> elements or a broad range of scales.
>
> I don't think LSC has gotten a lot of use since everyone I know who tried
> it has been sort of disappointed relative to other methods (e.g., inverse
> viscosity scaled mass matrix for heterogeneous Stokes, PCD for moderate Re
> Navier-Stokes). Of course there are no steady solutions to high Re so you
> either have a turbulence model or are time stepping. I'm not aware of work
> with LSC with turbulence models, and when time stepping, you have enough
> mass matrix that cheaper preconditioners are good enough. That said, it
> would be a great contribution to support this scaling.
>
> > I have also already tried using -ksp_diagonal_scale but the results are
> > identical.
>
> That's expected, and likely to mess up some MG implementations so I
> wouldn't recommend it.
>

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Jed Brown]

Sebastian Blauth  writes:

> Hello everyone,
>
> I wanted to briefly follow up on my question (see my last reply).
> Does anyone know / have an idea why the LSC preconditioner in PETSc does 
> not seem to scale well with the problem size (the outer fgmres solver I 
> am using nearly scale nearly linearly with the problem size in my example).

The implementation was tested on heterogeneous Stokes problems from 
geodynamics, and perhaps not on NS (or not with the discretization you're 
using).

https://doi.org/10.1016/j.pepi.2008.07.036

There is a comment about not having plumbing to provide a mass matrix. A few 
lines earlier there is code using PetscObjectQuery, and that same pattern could 
be applied for the mass matrix. If you're on a roughly uniform mesh, including 
the mass scaling will probably have little effect, but it could have a big 
impact in the presence of highly anistropic elements or a broad range of scales.

I don't think LSC has gotten a lot of use since everyone I know who tried it 
has been sort of disappointed relative to other methods (e.g., inverse 
viscosity scaled mass matrix for heterogeneous Stokes, PCD for moderate Re 
Navier-Stokes). Of course there are no steady solutions to high Re so you 
either have a turbulence model or are time stepping. I'm not aware of work with 
LSC with turbulence models, and when time stepping, you have enough mass matrix 
that cheaper preconditioners are good enough. That said, it would be a great 
contribution to support this scaling.

> I have also already tried using -ksp_diagonal_scale but the results are 
> identical.

That's expected, and likely to mess up some MG implementations so I wouldn't 
recommend it.

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Matthew Knepley]

On Mon, May 8, 2023 at 7:27 AM Mark Adams  wrote:

>
>
> On Mon, May 8, 2023 at 2:32 AM Sebastian Blauth <
> sebastian.bla...@itwm.fraunhofer.de> wrote:
>
>> Hello everyone,
>>
>> I wanted to briefly follow up on my question (see my last reply).
>> Does anyone know / have an idea why the LSC preconditioner in PETSc does
>> not seem to scale well with the problem size (the outer fgmres solver I
>> am using nearly scale nearly linearly with the problem size in my
>> example).
>> I have also already tried using -ksp_diagonal_scale but the results are
>> identical.
>> Any help is really appreciated.
>>
>
> I would start by finding results in the literature that you like, in that
> they are on a problem similar to yours and you like the performance, and
> reproduce them in your code.
> If you can do that then you have a research problem to see how to get your
> problem to work well.
> If your solver does not reproduce published results then we might be
> able to provide some advice.
>

I have not seen LSC scale well. The only thing I have seen to actually work
is patch smoothing (I gave a paper reference). It could
be that we have a bug in LSC, but I thought we verified it with the
Shuttleworth paper.

  THanks,

 Matt


> Mark
>
>
>>
>> Thanks a lot,
>> Sebastian
>>
>> On 03.05.2023 09:07, Sebastian Blauth wrote:
>> > First of all, yes you are correct that I am trying to solve the
>> > stationary incompressible Navier Stokes equations.
>> >
>> > On 02.05.2023 21:33, Matthew Knepley wrote:
>> >> On Tue, May 2, 2023 at 2:29 PM Jed Brown > >> > wrote:
>> >>
>> >> Sebastian Blauth > >> > writes:
>> >>
>> >>  > I agree with your comment for the Stokes equations - for these,
>> I
>> >> have
>> >>  > already tried and used the pressure mass matrix as part of a
>> >> (additive)
>> >>  > block preconditioner and it gave mesh independent results.
>> >>  >
>> >>  > However, for the Navier Stokes equations, is the Schur
>> complement
>> >> really
>> >>  > spectrally equivalent to the pressure mass matrix?
>> >>
>> >> No, it's not. You'd want something like PCD (better, but not
>> >> algebraic) or LSC.
>> >>
>> >
>> > I would like to take a look at the LSC preconditioner. For this, I did
>> > also not achieve mesh-independent results. I am using the following
>> > options (I know that the tolerances are too high at the moment, but it
>> > should just illustrate the behavior w.r.t. mesh refinement). Again, I
>> am
>> > using a simple 2D channel problem for testing purposes.
>> >
>> > I am using the following options
>> >
>> > -ksp_type fgmres
>> > -ksp_gmres_restart 100
>> > -ksp_gmres_cgs_refinement_type refine_ifneeded
>> > -ksp_max_it 1000
>> > -ksp_rtol 1e-10
>> > -ksp_atol 1e-30
>> > -pc_type fieldsplit
>> > -pc_fieldsplit_type schur
>> > -pc_fieldsplit_schur_fact_type full
>> > -pc_fieldsplit_schur_precondition self
>> > -fieldsplit_0_ksp_type preonly
>> > -fieldsplit_0_pc_type lu
>> > -fieldsplit_1_ksp_type gmres
>> > -fieldsplit_1_ksp_pc_side right
>> > -fieldsplit_1_ksp_gmres_restart 100
>> > -fieldsplit_1_ksp_gmres_cgs_refinement_type refine_ifneeded
>> > -fieldsplit_1_ksp_max_it 1000
>> > -fieldsplit_1_ksp_rtol 1e-10
>> > -fieldsplit_1_ksp_atol 1e-30
>> > -fieldsplit_1_pc_type lsc
>> > -fieldsplit_1_lsc_ksp_type preonly
>> > -fieldsplit_1_lsc_pc_type lu
>> > -fieldsplit_1_ksp_converged_reason
>> >
>> > Again, the direct solvers are used so that only the influence of the
>> LSC
>> > preconditioner is seen. I have suitable preconditioners for all of
>> these
>> > available (using boomeramg).
>> >
>> > At the bottom, I attach the output for different discretizations. As
>> you
>> > can see there, the number of iterations increases nearly linearly with
>> > the problem size.
>> >
>> > I think that the problem could occur due to a wrong scaling. In your
>> > docs https://petsc.org/release/manualpages/PC/PCLSC/ , you write that
>> > the LSC preconditioner is implemented as
>> >
>> > inv(S) \approx inv(A10 A01) A10 A00 A01 inv(A10 A01)
>> >
>> > However, in the book of Elman, Sylvester and Wathen (Finite Elements
>> and
>> > Fast Iterative Solvers), the LSC preconditioner is defined as
>> >
>> >  inv(S) \approx inv(A10 inv(T) A01) A10 inv(T) A00 inv(T) A01
>> > inv(A10 inv(T) A01)
>> >
>> > where T = diag(Q) and Q is the velocity mass matrix.
>> >
>> > There is an options -pc_lsc_scale_diag, which states that it uses the
>> > diagonal of A for scaling. I suppose, that this means, that the
>> diagonal
>> > of the A00 block is used for scaling - however, this is not the right
>> > scaling, is it? Even in the source code for the LSC preconditioner, in
>> > /src/ksp/pc/impls/lsc/lsc.c it is mentioned, that a mass matrix should
>> > be used...
>> > Is there any way to implement this in PETSc? Maybe by supplying the
>> mass
>> > matrix as Pmat?
>> >
>> > Thanks a lot in advance,
>> > Sebastian

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Mark Adams]

On Mon, May 8, 2023 at 2:32 AM Sebastian Blauth <
sebastian.bla...@itwm.fraunhofer.de> wrote:

> Hello everyone,
>
> I wanted to briefly follow up on my question (see my last reply).
> Does anyone know / have an idea why the LSC preconditioner in PETSc does
> not seem to scale well with the problem size (the outer fgmres solver I
> am using nearly scale nearly linearly with the problem size in my example).
> I have also already tried using -ksp_diagonal_scale but the results are
> identical.
> Any help is really appreciated.
>

I would start by finding results in the literature that you like, in that
they are on a problem similar to yours and you like the performance, and
reproduce them in your code.
If you can do that then you have a research problem to see how to get your
problem to work well.
If your solver does not reproduce published results then we might be
able to provide some advice.

Mark


>
> Thanks a lot,
> Sebastian
>
> On 03.05.2023 09:07, Sebastian Blauth wrote:
> > First of all, yes you are correct that I am trying to solve the
> > stationary incompressible Navier Stokes equations.
> >
> > On 02.05.2023 21:33, Matthew Knepley wrote:
> >> On Tue, May 2, 2023 at 2:29 PM Jed Brown  >> > wrote:
> >>
> >> Sebastian Blauth  >> > writes:
> >>
> >>  > I agree with your comment for the Stokes equations - for these, I
> >> have
> >>  > already tried and used the pressure mass matrix as part of a
> >> (additive)
> >>  > block preconditioner and it gave mesh independent results.
> >>  >
> >>  > However, for the Navier Stokes equations, is the Schur complement
> >> really
> >>  > spectrally equivalent to the pressure mass matrix?
> >>
> >> No, it's not. You'd want something like PCD (better, but not
> >> algebraic) or LSC.
> >>
> >
> > I would like to take a look at the LSC preconditioner. For this, I did
> > also not achieve mesh-independent results. I am using the following
> > options (I know that the tolerances are too high at the moment, but it
> > should just illustrate the behavior w.r.t. mesh refinement). Again, I am
> > using a simple 2D channel problem for testing purposes.
> >
> > I am using the following options
> >
> > -ksp_type fgmres
> > -ksp_gmres_restart 100
> > -ksp_gmres_cgs_refinement_type refine_ifneeded
> > -ksp_max_it 1000
> > -ksp_rtol 1e-10
> > -ksp_atol 1e-30
> > -pc_type fieldsplit
> > -pc_fieldsplit_type schur
> > -pc_fieldsplit_schur_fact_type full
> > -pc_fieldsplit_schur_precondition self
> > -fieldsplit_0_ksp_type preonly
> > -fieldsplit_0_pc_type lu
> > -fieldsplit_1_ksp_type gmres
> > -fieldsplit_1_ksp_pc_side right
> > -fieldsplit_1_ksp_gmres_restart 100
> > -fieldsplit_1_ksp_gmres_cgs_refinement_type refine_ifneeded
> > -fieldsplit_1_ksp_max_it 1000
> > -fieldsplit_1_ksp_rtol 1e-10
> > -fieldsplit_1_ksp_atol 1e-30
> > -fieldsplit_1_pc_type lsc
> > -fieldsplit_1_lsc_ksp_type preonly
> > -fieldsplit_1_lsc_pc_type lu
> > -fieldsplit_1_ksp_converged_reason
> >
> > Again, the direct solvers are used so that only the influence of the LSC
> > preconditioner is seen. I have suitable preconditioners for all of these
> > available (using boomeramg).
> >
> > At the bottom, I attach the output for different discretizations. As you
> > can see there, the number of iterations increases nearly linearly with
> > the problem size.
> >
> > I think that the problem could occur due to a wrong scaling. In your
> > docs https://petsc.org/release/manualpages/PC/PCLSC/ , you write that
> > the LSC preconditioner is implemented as
> >
> > inv(S) \approx inv(A10 A01) A10 A00 A01 inv(A10 A01)
> >
> > However, in the book of Elman, Sylvester and Wathen (Finite Elements and
> > Fast Iterative Solvers), the LSC preconditioner is defined as
> >
> >  inv(S) \approx inv(A10 inv(T) A01) A10 inv(T) A00 inv(T) A01
> > inv(A10 inv(T) A01)
> >
> > where T = diag(Q) and Q is the velocity mass matrix.
> >
> > There is an options -pc_lsc_scale_diag, which states that it uses the
> > diagonal of A for scaling. I suppose, that this means, that the diagonal
> > of the A00 block is used for scaling - however, this is not the right
> > scaling, is it? Even in the source code for the LSC preconditioner, in
> > /src/ksp/pc/impls/lsc/lsc.c it is mentioned, that a mass matrix should
> > be used...
> > Is there any way to implement this in PETSc? Maybe by supplying the mass
> > matrix as Pmat?
> >
> > Thanks a lot in advance,
> > Sebastian
> >
> >>
> >> I think you can do a better job than that using something like
> >>
> >> https://arxiv.org/abs/1810.03315 
> >>
> >> Basically, you use an augmented Lagrangian thing to make the Schur
> >> complement well-conditioned,
> >> and then use a special smoother to handle that perturbation.
> >>
> >>  > And even if it is, the convergence is only good for small
> >> Reynolds numbers, for moderately high ones 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Sebastian Blauth]

Hello everyone,

I wanted to briefly follow up on my question (see my last reply).
Does anyone know / have an idea why the LSC preconditioner in PETSc does 
not seem to scale well with the problem size (the outer fgmres solver I 
am using nearly scale nearly linearly with the problem size in my example).
I have also already tried using -ksp_diagonal_scale but the results are 
identical.

Any help is really appreciated.

Thanks a lot,
Sebastian

On 03.05.2023 09:07, Sebastian Blauth wrote:
First of all, yes you are correct that I am trying to solve the 
stationary incompressible Navier Stokes equations.


On 02.05.2023 21:33, Matthew Knepley wrote:
On Tue, May 2, 2023 at 2:29 PM Jed Brown > wrote:


    Sebastian Blauth mailto:sebastian.bla...@itwm.fraunhofer.de>> writes:

 > I agree with your comment for the Stokes equations - for these, I
    have
 > already tried and used the pressure mass matrix as part of a
    (additive)
 > block preconditioner and it gave mesh independent results.
 >
 > However, for the Navier Stokes equations, is the Schur complement
    really
 > spectrally equivalent to the pressure mass matrix?

    No, it's not. You'd want something like PCD (better, but not
    algebraic) or LSC.



I would like to take a look at the LSC preconditioner. For this, I did 
also not achieve mesh-independent results. I am using the following 
options (I know that the tolerances are too high at the moment, but it 
should just illustrate the behavior w.r.t. mesh refinement). Again, I am 
using a simple 2D channel problem for testing purposes.


I am using the following options

-ksp_type fgmres
-ksp_gmres_restart 100
-ksp_gmres_cgs_refinement_type refine_ifneeded
-ksp_max_it 1000
-ksp_rtol 1e-10
-ksp_atol 1e-30
-pc_type fieldsplit
-pc_fieldsplit_type schur
-pc_fieldsplit_schur_fact_type full
-pc_fieldsplit_schur_precondition self
-fieldsplit_0_ksp_type preonly
-fieldsplit_0_pc_type lu
-fieldsplit_1_ksp_type gmres
-fieldsplit_1_ksp_pc_side right
-fieldsplit_1_ksp_gmres_restart 100
-fieldsplit_1_ksp_gmres_cgs_refinement_type refine_ifneeded
-fieldsplit_1_ksp_max_it 1000
-fieldsplit_1_ksp_rtol 1e-10
-fieldsplit_1_ksp_atol 1e-30
-fieldsplit_1_pc_type lsc
-fieldsplit_1_lsc_ksp_type preonly
-fieldsplit_1_lsc_pc_type lu
-fieldsplit_1_ksp_converged_reason

Again, the direct solvers are used so that only the influence of the LSC 
preconditioner is seen. I have suitable preconditioners for all of these 
available (using boomeramg).


At the bottom, I attach the output for different discretizations. As you 
can see there, the number of iterations increases nearly linearly with 
the problem size.


I think that the problem could occur due to a wrong scaling. In your 
docs https://petsc.org/release/manualpages/PC/PCLSC/ , you write that 
the LSC preconditioner is implemented as


    inv(S) \approx inv(A10 A01) A10 A00 A01 inv(A10 A01)

However, in the book of Elman, Sylvester and Wathen (Finite Elements and 
Fast Iterative Solvers), the LSC preconditioner is defined as


     inv(S) \approx inv(A10 inv(T) A01) A10 inv(T) A00 inv(T) A01 
inv(A10 inv(T) A01)


where T = diag(Q) and Q is the velocity mass matrix.

There is an options -pc_lsc_scale_diag, which states that it uses the 
diagonal of A for scaling. I suppose, that this means, that the diagonal 
of the A00 block is used for scaling - however, this is not the right 
scaling, is it? Even in the source code for the LSC preconditioner, in 
/src/ksp/pc/impls/lsc/lsc.c it is mentioned, that a mass matrix should 
be used...
Is there any way to implement this in PETSc? Maybe by supplying the mass 
matrix as Pmat?


Thanks a lot in advance,
Sebastian



I think you can do a better job than that using something like

https://arxiv.org/abs/1810.03315 

Basically, you use an augmented Lagrangian thing to make the Schur 
complement well-conditioned,

and then use a special smoother to handle that perturbation.

 > And even if it is, the convergence is only good for small
    Reynolds numbers, for moderately high ones the convergence really
    deteriorates. This is why I am trying to make
    fieldsplit_schur_precondition selfp work better (this is, if I
    understand it correctly, a SIMPLE type preconditioner).

    SIMPLE is for short time steps (not too far from resolving CFL) and
    bad for steady. This taxonomy is useful, though the problems are
    super academic and they don't use high aspect ratio.



Okay, I get that I cannot expect the SIMPLE preconditioning 
(schur_precondition selfp) to work efficiently. I guess the reason it 
works for small time steps (for the instationary equations) is that the 
velocity block becomes diagonally dominant in this case, so that diag(A) 
is a good approximation of A.




    https://doi.org/10.1016/j.jcp.2007.09.026
    


    Thanks,

       Matt

--
What most experimenters take for 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Sebastian Blauth]

First of all, yes you are correct that I am trying to solve the 
stationary incompressible Navier Stokes equations.


On 02.05.2023 21:33, Matthew Knepley wrote:
On Tue, May 2, 2023 at 2:29 PM Jed Brown > wrote:


Sebastian Blauth mailto:sebastian.bla...@itwm.fraunhofer.de>> writes:

 > I agree with your comment for the Stokes equations - for these, I
have
 > already tried and used the pressure mass matrix as part of a
(additive)
 > block preconditioner and it gave mesh independent results.
 >
 > However, for the Navier Stokes equations, is the Schur complement
really
 > spectrally equivalent to the pressure mass matrix?

No, it's not. You'd want something like PCD (better, but not
algebraic) or LSC.



I would like to take a look at the LSC preconditioner. For this, I did 
also not achieve mesh-independent results. I am using the following 
options (I know that the tolerances are too high at the moment, but it 
should just illustrate the behavior w.r.t. mesh refinement). Again, I am 
using a simple 2D channel problem for testing purposes.


I am using the following options

-ksp_type fgmres
-ksp_gmres_restart 100
-ksp_gmres_cgs_refinement_type refine_ifneeded
-ksp_max_it 1000
-ksp_rtol 1e-10
-ksp_atol 1e-30
-pc_type fieldsplit
-pc_fieldsplit_type schur
-pc_fieldsplit_schur_fact_type full
-pc_fieldsplit_schur_precondition self
-fieldsplit_0_ksp_type preonly
-fieldsplit_0_pc_type lu
-fieldsplit_1_ksp_type gmres
-fieldsplit_1_ksp_pc_side right
-fieldsplit_1_ksp_gmres_restart 100
-fieldsplit_1_ksp_gmres_cgs_refinement_type refine_ifneeded
-fieldsplit_1_ksp_max_it 1000
-fieldsplit_1_ksp_rtol 1e-10
-fieldsplit_1_ksp_atol 1e-30
-fieldsplit_1_pc_type lsc
-fieldsplit_1_lsc_ksp_type preonly
-fieldsplit_1_lsc_pc_type lu
-fieldsplit_1_ksp_converged_reason

Again, the direct solvers are used so that only the influence of the LSC 
preconditioner is seen. I have suitable preconditioners for all of these 
available (using boomeramg).


At the bottom, I attach the output for different discretizations. As you 
can see there, the number of iterations increases nearly linearly with 
the problem size.


I think that the problem could occur due to a wrong scaling. In your 
docs https://petsc.org/release/manualpages/PC/PCLSC/ , you write that 
the LSC preconditioner is implemented as


   inv(S) \approx inv(A10 A01) A10 A00 A01 inv(A10 A01)

However, in the book of Elman, Sylvester and Wathen (Finite Elements and 
Fast Iterative Solvers), the LSC preconditioner is defined as


inv(S) \approx inv(A10 inv(T) A01) A10 inv(T) A00 inv(T) A01 
inv(A10 inv(T) A01)


where T = diag(Q) and Q is the velocity mass matrix.

There is an options -pc_lsc_scale_diag, which states that it uses the 
diagonal of A for scaling. I suppose, that this means, that the diagonal 
of the A00 block is used for scaling - however, this is not the right 
scaling, is it? Even in the source code for the LSC preconditioner, in 
/src/ksp/pc/impls/lsc/lsc.c it is mentioned, that a mass matrix should 
be used...
Is there any way to implement this in PETSc? Maybe by supplying the mass 
matrix as Pmat?


Thanks a lot in advance,
Sebastian



I think you can do a better job than that using something like

https://arxiv.org/abs/1810.03315 

Basically, you use an augmented Lagrangian thing to make the Schur 
complement well-conditioned,

and then use a special smoother to handle that perturbation.

 > And even if it is, the convergence is only good for small
Reynolds numbers, for moderately high ones the convergence really
deteriorates. This is why I am trying to make
fieldsplit_schur_precondition selfp work better (this is, if I
understand it correctly, a SIMPLE type preconditioner).

SIMPLE is for short time steps (not too far from resolving CFL) and
bad for steady. This taxonomy is useful, though the problems are
super academic and they don't use high aspect ratio.



Okay, I get that I cannot expect the SIMPLE preconditioning 
(schur_precondition selfp) to work efficiently. I guess the reason it 
works for small time steps (for the instationary equations) is that the 
velocity block becomes diagonally dominant in this case, so that diag(A) 
is a good approximation of A.




https://doi.org/10.1016/j.jcp.2007.09.026



    Thanks,

       Matt

--
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/ 



And here is the output of my scaling tests

8x8 discretization

Newton solver:  iter,  abs. residual (abs. tol),  rel. residual (rel. tol)

Newton solver: 0,  1.023e+03 (1.00e-30),  1.000e+00 (1.00e-10)
  Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Matthew Knepley]

On Tue, May 2, 2023 at 2:29 PM Jed Brown  wrote:

> Sebastian Blauth  writes:
>
> > I agree with your comment for the Stokes equations - for these, I have
> > already tried and used the pressure mass matrix as part of a (additive)
> > block preconditioner and it gave mesh independent results.
> >
> > However, for the Navier Stokes equations, is the Schur complement really
> > spectrally equivalent to the pressure mass matrix?
>
> No, it's not. You'd want something like PCD (better, but not algebraic) or
> LSC.
>

I think you can do a better job than that using something like

  https://arxiv.org/abs/1810.03315

Basically, you use an augmented Lagrangian thing to make the Schur
complement well-conditioned,
and then use a special smoother to handle that perturbation.


> > And even if it is, the convergence is only good for small Reynolds
> numbers, for moderately high ones the convergence really deteriorates. This
> is why I am trying to make fieldsplit_schur_precondition selfp work better
> (this is, if I understand it correctly, a SIMPLE type preconditioner).
>
> SIMPLE is for short time steps (not too far from resolving CFL) and bad
> for steady. This taxonomy is useful, though the problems are super academic
> and they don't use high aspect ratio.
>
> https://doi.org/10.1016/j.jcp.2007.09.026
>

   Thanks,

  Matt

-- 
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-users] Scalable Solver for Incompressible Flow

[Jed Brown]

Sebastian Blauth  writes:

> I agree with your comment for the Stokes equations - for these, I have 
> already tried and used the pressure mass matrix as part of a (additive) 
> block preconditioner and it gave mesh independent results.
>
> However, for the Navier Stokes equations, is the Schur complement really 
> spectrally equivalent to the pressure mass matrix? 

No, it's not. You'd want something like PCD (better, but not algebraic) or LSC.

> And even if it is, the convergence is only good for small Reynolds numbers, 
> for moderately high ones the convergence really deteriorates. This is why I 
> am trying to make fieldsplit_schur_precondition selfp work better (this is, 
> if I understand it correctly, a SIMPLE type preconditioner).

SIMPLE is for short time steps (not too far from resolving CFL) and bad for 
steady. This taxonomy is useful, though the problems are super academic and 
they don't use high aspect ratio.

https://doi.org/10.1016/j.jcp.2007.09.026

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Sebastian Blauth]

On 02.05.2023 15:12, Matthew Knepley wrote:
On Tue, May 2, 2023 at 9:07 AM Blauth, Sebastian 
> wrote:


Hello,

__ __

I am having a problem using / configuring PETSc to obtain a scalable
solver for the incompressible Navier Stokes equations. I am
discretizing the equations using FEM (with the library fenics) and I
am using the stable P2-P1 Taylor-Hood elements. I have read and
tried a lot regarding preconditioners for incompressible Navier
Stokes and I am aware that this is very much an active research
field, but maybe I can get some hints / tips. 

I am interested in solving large-scale 3D problems, but I cannot
even set up a scaleable 2D solver for the problems. All of my
approaches at the moment are trying to use a Schur Complement
approach, but I cannot get a “good” preconditioner for the Schur
complement matrix. For the velocity block, I am using the AMG
provided by hypre (which seems to work fine and is most likely not
the problem).

__ __

To test the solver, I am using a simple 2D channel flow problem with
do-nothing conditions at the outlet.

__ __

I am facing the following difficulties at the moment:

__ __

- First, I am having trouble with using
-pc_fieldsplit_schur_precondition selfp. With this setup, the cost
for solving the Schur complement part in the fieldsplit
preconditioner (approximately) increase when the mesh is refined. I
am using the following options for this setup (note that I am using
exact solves for the velocity part to debug, but using, e.g., gmres
with hypre boomeramg reaches a given tolerance with a number of
iterations that is independent of the mesh)


The diagonal of the momentum block is a bad preconditioner for the Schur 
complement, because S is spectrally equivalent to the mass matrix. You 
should build the mass matrix and use that as the preconditioning matrix 
for the Schur part. The FEniCS people can show you how to do that. This 
will provide mesh-independent convergence (you can see me doing this in 
SNES ex69).


   Thanks,

      Matt


I agree with your comment for the Stokes equations - for these, I have 
already tried and used the pressure mass matrix as part of a (additive) 
block preconditioner and it gave mesh independent results.


However, for the Navier Stokes equations, is the Schur complement really 
spectrally equivalent to the pressure mass matrix? And even if it is, 
the convergence is only good for small Reynolds numbers, for moderately 
high ones the convergence really deteriorates. This is why I am trying 
to make fieldsplit_schur_precondition selfp work better (this is, if I 
understand it correctly, a SIMPLE type preconditioner).


Best regards,
Sebastian




     -ksp_type fgmres

     -ksp_rtol 1e-6

-ksp_atol 1e-30

     -pc_type fieldsplit

     -pc_fieldsplit_type schur

     -pc_fieldsplit_schur_fact_type full

     -pc_fieldsplit_schur_precondition selfp

     -fieldsplit_0_ksp_type preonly

     -fieldsplit_0_pc_type lu

     -fieldsplit_1_ksp_type gmres

     -fieldsplit_1_ksp_pc_side right

     -fieldsplit_1_ksp_max_it 1000

     -fieldsplit_1_ksp_rtol 1e-1

     -fieldsplit_1_ksp_atol 1e-30

     -fieldsplit_1_pc_type lu

     -fieldsplit_1_ksp_converged_reason

     -ksp_converged_reason

__ __

Note, that I use direct solvers for the subproblems to get an
“ideal” convergence. Even if I replace the direct solver with
boomeramg, the behavior is the same and the number of iterations
does not change much. 

In particular, I get the following behavior:

For a 8x8 mesh, I need, on average, 25 iterations to solve
fieldsplit_1

For a 16x16 mesh, I need 40 iterations

For a 32x32 mesh, I need 70 iterations

For a 64x64 mesh, I need 100 iterations

__ __

However, the outer fgmres requires, as expected, always the same
number of iterations to reach convergence (as expected).

I do understand that the selfp preconditioner for the Schur
complement is expected to deteriorate as the Reynolds number
increases and the problem becomes more convective in nature, but I
had hoped that I can at least get a scaleable preconditioner with
respect to the mesh size out of it. Are there any tips on how to
achieve this?

__ __

My second problem is concerning the LSC preconditioner. When I am
using this, again both with exact solves of the linear problems or
when using boomeramg, I do not get a scalable solver with respect to
the mesh size. On the contrary, here the number of solves required
for solving fieldsplit_1 to a fixed relative tolerance seem to
behave linearly w.r.t. the problem size. For 

Re: [petsc-users] Scalable Solver for Incompressible Flow

[Matthew Knepley]

On Tue, May 2, 2023 at 9:07 AM Blauth, Sebastian <
sebastian.bla...@itwm.fraunhofer.de> wrote:

> Hello,
>
>
>
> I am having a problem using / configuring PETSc to obtain a scalable
> solver for the incompressible Navier Stokes equations. I am discretizing
> the equations using FEM (with the library fenics) and I am using the stable
> P2-P1 Taylor-Hood elements. I have read and tried a lot regarding
> preconditioners for incompressible Navier Stokes and I am aware that this
> is very much an active research field, but maybe I can get some hints /
> tips.
>
> I am interested in solving large-scale 3D problems, but I cannot even set
> up a scaleable 2D solver for the problems. All of my approaches at the
> moment are trying to use a Schur Complement approach, but I cannot get a
> “good” preconditioner for the Schur complement matrix. For the velocity
> block, I am using the AMG provided by hypre (which seems to work fine and
> is most likely not the problem).
>
>
>
> To test the solver, I am using a simple 2D channel flow problem with
> do-nothing conditions at the outlet.
>
>
>
> I am facing the following difficulties at the moment:
>
>
>
> - First, I am having trouble with using -pc_fieldsplit_schur_precondition
> selfp. With this setup, the cost for solving the Schur complement part in
> the fieldsplit preconditioner (approximately) increase when the mesh is
> refined. I am using the following options for this setup (note that I am
> using exact solves for the velocity part to debug, but using, e.g., gmres
> with hypre boomeramg reaches a given tolerance with a number of iterations
> that is independent of the mesh)
>

The diagonal of the momentum block is a bad preconditioner for the Schur
complement, because S is spectrally equivalent to the mass matrix. You
should build the mass matrix and use that as the preconditioning matrix for
the Schur part. The FEniCS people can show you how to do that. This will
provide mesh-independent convergence (you can see me doing this in SNES
ex69).

  Thanks,

 Matt


> -ksp_type fgmres
>
> -ksp_rtol 1e-6
>
> -ksp_atol 1e-30
>
> -pc_type fieldsplit
>
> -pc_fieldsplit_type schur
>
> -pc_fieldsplit_schur_fact_type full
>
> -pc_fieldsplit_schur_precondition selfp
>
> -fieldsplit_0_ksp_type preonly
>
> -fieldsplit_0_pc_type lu
>
> -fieldsplit_1_ksp_type gmres
>
> -fieldsplit_1_ksp_pc_side right
>
> -fieldsplit_1_ksp_max_it 1000
>
> -fieldsplit_1_ksp_rtol 1e-1
>
> -fieldsplit_1_ksp_atol 1e-30
>
> -fieldsplit_1_pc_type lu
>
> -fieldsplit_1_ksp_converged_reason
>
> -ksp_converged_reason
>
>
>
> Note, that I use direct solvers for the subproblems to get an “ideal”
> convergence. Even if I replace the direct solver with boomeramg, the
> behavior is the same and the number of iterations does not change much.
>
> In particular, I get the following behavior:
>
> For a 8x8 mesh, I need, on average, 25 iterations to solve fieldsplit_1
>
> For a 16x16 mesh, I need 40 iterations
>
> For a 32x32 mesh, I need 70 iterations
>
> For a 64x64 mesh, I need 100 iterations
>
>
>
> However, the outer fgmres requires, as expected, always the same number of
> iterations to reach convergence (as expected).
>
> I do understand that the selfp preconditioner for the Schur complement is
> expected to deteriorate as the Reynolds number increases and the problem
> becomes more convective in nature, but I had hoped that I can at least get
> a scaleable preconditioner with respect to the mesh size out of it. Are
> there any tips on how to achieve this?
>
>
>
> My second problem is concerning the LSC preconditioner. When I am using
> this, again both with exact solves of the linear problems or when using
> boomeramg, I do not get a scalable solver with respect to the mesh size. On
> the contrary, here the number of solves required for solving fieldsplit_1
> to a fixed relative tolerance seem to behave linearly w.r.t. the problem
> size. For this problem, I suspect that the issue lies in the scaling of the
> LSC preconditioner matrices (in the book of Elman, Sylvester and Wathen,
> the matrices are scaled with the inverse of the diagonal velocity mass
> matrix). Is it possible to achieve this with PETSc? I started experimenting
> with supplying the velocity mass matrix as preconditioner matrix and using
> “use_amat”, but I am not sure where / how to do it this way.
>
>
>
> And finally, more of an observation and question: I noticed that the AMG
> approximations for the velocity block became worse with increase of the
> Reynolds number when using the default options. However, when using
> -pc_hypre_boomeramg_relax_weight_all 0.0 I noticed that boomeramg performed
> way more robustly w.r.t. the Reynolds number. Are there any other ways to
> improve the AMG performance in this regard?
>
>
>
> Thanks a lot in advance and I am looking forward to your reply,
>
> Sebastian
>
>
>
> --
>
> Dr. Sebastian Blauth
>
> Fraunhofer-Institut für
>
> 

[petsc-users] Scalable Solver for Incompressible Flow

[Blauth, Sebastian]

Hello,

 

I am having a problem using / configuring PETSc to obtain a scalable solver
for the incompressible Navier Stokes equations. I am discretizing the
equations using FEM (with the library fenics) and I am using the stable
P2-P1 Taylor-Hood elements. I have read and tried a lot regarding
preconditioners for incompressible Navier Stokes and I am aware that this is
very much an active research field, but maybe I can get some hints / tips. 

I am interested in solving large-scale 3D problems, but I cannot even set up
a scaleable 2D solver for the problems. All of my approaches at the moment
are trying to use a Schur Complement approach, but I cannot get a “good”
preconditioner for the Schur complement matrix. For the velocity block, I am
using the AMG provided by hypre (which seems to work fine and is most likely
not the problem).

 

To test the solver, I am using a simple 2D channel flow problem with
do-nothing conditions at the outlet.

 

I am facing the following difficulties at the moment:

 

- First, I am having trouble with using -pc_fieldsplit_schur_precondition
selfp. With this setup, the cost for solving the Schur complement part in
the fieldsplit preconditioner (approximately) increase when the mesh is
refined. I am using the following options for this setup (note that I am
using exact solves for the velocity part to debug, but using, e.g., gmres
with hypre boomeramg reaches a given tolerance with a number of iterations
that is independent of the mesh)

 

-ksp_type fgmres

-ksp_rtol 1e-6

-ksp_atol 1e-30

-pc_type fieldsplit

-pc_fieldsplit_type schur

-pc_fieldsplit_schur_fact_type full

-pc_fieldsplit_schur_precondition selfp

-fieldsplit_0_ksp_type preonly

-fieldsplit_0_pc_type lu

-fieldsplit_1_ksp_type gmres

-fieldsplit_1_ksp_pc_side right

-fieldsplit_1_ksp_max_it 1000

-fieldsplit_1_ksp_rtol 1e-1

-fieldsplit_1_ksp_atol 1e-30

-fieldsplit_1_pc_type lu

-fieldsplit_1_ksp_converged_reason

-ksp_converged_reason

 

Note, that I use direct solvers for the subproblems to get an “ideal”
convergence. Even if I replace the direct solver with boomeramg, the
behavior is the same and the number of iterations does not change much. 

In particular, I get the following behavior:

For a 8x8 mesh, I need, on average, 25 iterations to solve fieldsplit_1

For a 16x16 mesh, I need 40 iterations

For a 32x32 mesh, I need 70 iterations

For a 64x64 mesh, I need 100 iterations

 

However, the outer fgmres requires, as expected, always the same number of
iterations to reach convergence (as expected).

I do understand that the selfp preconditioner for the Schur complement is
expected to deteriorate as the Reynolds number increases and the problem
becomes more convective in nature, but I had hoped that I can at least get a
scaleable preconditioner with respect to the mesh size out of it. Are there
any tips on how to achieve this?

 

My second problem is concerning the LSC preconditioner. When I am using
this, again both with exact solves of the linear problems or when using
boomeramg, I do not get a scalable solver with respect to the mesh size. On
the contrary, here the number of solves required for solving fieldsplit_1 to
a fixed relative tolerance seem to behave linearly w.r.t. the problem size.
For this problem, I suspect that the issue lies in the scaling of the LSC
preconditioner matrices (in the book of Elman, Sylvester and Wathen, the
matrices are scaled with the inverse of the diagonal velocity mass matrix).
Is it possible to achieve this with PETSc? I started experimenting with
supplying the velocity mass matrix as preconditioner matrix and using
“use_amat”, but I am not sure where / how to do it this way.

 

And finally, more of an observation and question: I noticed that the AMG
approximations for the velocity block became worse with increase of the
Reynolds number when using the default options. However, when using
-pc_hypre_boomeramg_relax_weight_all 0.0 I noticed that boomeramg performed
way more robustly w.r.t. the Reynolds number. Are there any other ways to
improve the AMG performance in this regard?

 

Thanks a lot in advance and I am looking forward to your reply,

Sebastian

 

--

Dr. Sebastian Blauth

Fraunhofer-Institut für

Techno- und Wirtschaftsmathematik ITWM

Abteilung Transportvorgänge

Fraunhofer-Platz 1, 67663 Kaiserslautern

Telefon: +49 631 31600-4968

sebastian.bla...@itwm.fraunhofer.de

https://www.itwm.fraunhofer.de

 



smime.p7s
Description: S/MIME cryptographic signature