Scalable Solver for Incompressible Flow
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 [email protected] https://www.itwm.fraunhofer.de
On Tue, May 2, 2023 at 9:07 AM Blauth, Sebastian < [email protected]> 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
Techno- und Wirtschaftsmathematik ITWM
Abteilung Transportvorgänge
Fraunhofer-Platz 1, 67663 Kaiserslautern
Telefon: +49 631 31600-4968
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
On 02.05.2023 15:12, Matthew Knepley wrote:
On Tue, May 2, 2023 at 9:07 AM Blauth, Sebastian <[email protected] <mailto:[email protected]>> 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 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____
[email protected] <mailto:[email protected]>____
https://www.itwm.fraunhofer.de <https://www.itwm.fraunhofer.de>____
__ __
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-- Dr. Sebastian Blauth Fraunhofer-Institut für Techno- und Wirtschaftsmathematik ITWM Abteilung Transportvorgänge Fraunhofer-Platz 1, 67663 Kaiserslautern Telefon: +49 631 31600-4968 [email protected] www.itwm.fraunhofer.de
Sebastian Blauth <[email protected]> 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
On Tue, May 2, 2023 at 2:29 PM Jed Brown <[email protected]> wrote:
Sebastian Blauth <[email protected]> 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.
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/ <http://www.cse.buffalo.edu/~knepley/>
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 <[email protected] <mailto:[email protected]>> wrote:
Sebastian Blauth <[email protected] <mailto:[email protected]>> 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 <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 <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/ <http://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 iterations 38 Newton solver: 1, 1.313e+03 (1.00e-30), 1.283e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 76 Newton solver: 2, 1.198e+02 (1.00e-30), 1.171e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 3, 7.249e-01 (1.00e-30), 7.084e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 4, 3.883e-05 (1.00e-30), 3.795e-08 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 5, 2.778e-12 (1.00e-30), 2.714e-15 (1.00e-10) 16x16 discretization Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol) Newton solver: 0, 1.113e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 62 Newton solver: 1, 8.316e+02 (1.00e-30), 7.475e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 141 Newton solver: 2, 5.806e+01 (1.00e-30), 5.218e-02 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 119 Newton solver: 3, 3.309e-01 (1.00e-30), 2.974e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 118 Newton solver: 4, 9.085e-06 (1.00e-30), 8.166e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 120 Newton solver: 5, 3.475e-12 (1.00e-30), 3.124e-15 (1.00e-10) 32x32 discretization Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol) Newton solver: 0, 1.330e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 98 Newton solver: 1, 5.913e+02 (1.00e-30), 4.445e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 183 Newton solver: 2, 3.214e+01 (1.00e-30), 2.416e-02 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 152 Newton solver: 3, 2.059e-01 (1.00e-30), 1.547e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 151 Newton solver: 4, 6.949e-06 (1.00e-30), 5.223e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 149 Newton solver: 5, 5.300e-12 (1.00e-30), 3.983e-15 (1.00e-10) 64x64 discretization Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol) Newton solver: 0, 1.707e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 198 Newton solver: 1, 4.259e+02 (1.00e-30), 2.494e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 357 Newton solver: 2, 1.706e+01 (1.00e-30), 9.993e-03 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 266 Newton solver: 3, 1.134e-01 (1.00e-30), 6.639e-05 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 261 Newton solver: 4, 4.285e-06 (1.00e-30), 2.510e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 263 Newton solver: 5, 9.650e-12 (1.00e-30), 5.652e-15 (1.00e-10) -- Dr. Sebastian Blauth Fraunhofer-Institut für Techno- und Wirtschaftsmathematik ITWM Abteilung Transportvorgänge Fraunhofer-Platz 1, 67663 Kaiserslautern Telefon: +49 631 31600-4968 [email protected] www.itwm.fraunhofer.de
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 <[email protected] <mailto:[email protected]>> wrote:
Sebastian Blauth <[email protected] <mailto:[email protected]>> 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 <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 <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/ <http://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 iterations 38 Newton solver: 1, 1.313e+03 (1.00e-30), 1.283e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 76 Newton solver: 2, 1.198e+02 (1.00e-30), 1.171e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 3, 7.249e-01 (1.00e-30), 7.084e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 4, 3.883e-05 (1.00e-30), 3.795e-08 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 5, 2.778e-12 (1.00e-30), 2.714e-15 (1.00e-10)
16x16 discretization
Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol)
Newton solver: 0, 1.113e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 62 Newton solver: 1, 8.316e+02 (1.00e-30), 7.475e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 141 Newton solver: 2, 5.806e+01 (1.00e-30), 5.218e-02 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 119 Newton solver: 3, 3.309e-01 (1.00e-30), 2.974e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 118 Newton solver: 4, 9.085e-06 (1.00e-30), 8.166e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 120 Newton solver: 5, 3.475e-12 (1.00e-30), 3.124e-15 (1.00e-10)
32x32 discretization
Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol)
Newton solver: 0, 1.330e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 98 Newton solver: 1, 5.913e+02 (1.00e-30), 4.445e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 183 Newton solver: 2, 3.214e+01 (1.00e-30), 2.416e-02 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 152 Newton solver: 3, 2.059e-01 (1.00e-30), 1.547e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 151 Newton solver: 4, 6.949e-06 (1.00e-30), 5.223e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 149 Newton solver: 5, 5.300e-12 (1.00e-30), 3.983e-15 (1.00e-10)
64x64 discretization
Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol)
Newton solver: 0, 1.707e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 198 Newton solver: 1, 4.259e+02 (1.00e-30), 2.494e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 357 Newton solver: 2, 1.706e+01 (1.00e-30), 9.993e-03 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 266 Newton solver: 3, 1.134e-01 (1.00e-30), 6.639e-05 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 261 Newton solver: 4, 4.285e-06 (1.00e-30), 2.510e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 263 Newton solver: 5, 9.650e-12 (1.00e-30), 5.652e-15 (1.00e-10)
-- Dr. Sebastian Blauth Fraunhofer-Institut für Techno- und Wirtschaftsmathematik ITWM Abteilung Transportvorgänge Fraunhofer-Platz 1, 67663 Kaiserslautern Telefon: +49 631 31600-4968 [email protected] www.itwm.fraunhofer.de
On Mon, May 8, 2023 at 2:32 AM Sebastian Blauth < [email protected]> 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 <[email protected] <mailto:[email protected]>> wrote:
Sebastian Blauth <[email protected] <mailto:[email protected]>> 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 <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 <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/ <http://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 iterations 38 Newton solver: 1, 1.313e+03 (1.00e-30), 1.283e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 76 Newton solver: 2, 1.198e+02 (1.00e-30), 1.171e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 3, 7.249e-01 (1.00e-30), 7.084e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 4, 3.883e-05 (1.00e-30), 3.795e-08 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 5, 2.778e-12 (1.00e-30), 2.714e-15 (1.00e-10)
16x16 discretization
Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol)
Newton solver: 0, 1.113e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 62 Newton solver: 1, 8.316e+02 (1.00e-30), 7.475e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 141 Newton solver: 2, 5.806e+01 (1.00e-30), 5.218e-02 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 119 Newton solver: 3, 3.309e-01 (1.00e-30), 2.974e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 118 Newton solver: 4, 9.085e-06 (1.00e-30), 8.166e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 120 Newton solver: 5, 3.475e-12 (1.00e-30), 3.124e-15 (1.00e-10)
32x32 discretization
Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol)
Newton solver: 0, 1.330e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 98 Newton solver: 1, 5.913e+02 (1.00e-30), 4.445e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 183 Newton solver: 2, 3.214e+01 (1.00e-30), 2.416e-02 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 152 Newton solver: 3, 2.059e-01 (1.00e-30), 1.547e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 151 Newton solver: 4, 6.949e-06 (1.00e-30), 5.223e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 149 Newton solver: 5, 5.300e-12 (1.00e-30), 3.983e-15 (1.00e-10)
64x64 discretization
Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol)
Newton solver: 0, 1.707e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 198 Newton solver: 1, 4.259e+02 (1.00e-30), 2.494e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 357 Newton solver: 2, 1.706e+01 (1.00e-30), 9.993e-03 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 266 Newton solver: 3, 1.134e-01 (1.00e-30), 6.639e-05 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 261 Newton solver: 4, 4.285e-06 (1.00e-30), 2.510e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 263 Newton solver: 5, 9.650e-12 (1.00e-30), 5.652e-15 (1.00e-10)
-- Dr. Sebastian Blauth Fraunhofer-Institut für Techno- und Wirtschaftsmathematik ITWM Abteilung Transportvorgänge Fraunhofer-Platz 1, 67663 Kaiserslautern Telefon: +49 631 31600-4968 [email protected] www.itwm.fraunhofer.de
On Mon, May 8, 2023 at 7:27 AM Mark Adams <[email protected]> wrote:
On Mon, May 8, 2023 at 2:32 AM Sebastian Blauth < [email protected]> 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 <[email protected] <mailto:[email protected]>> wrote:
Sebastian Blauth <[email protected] <mailto:[email protected]>> 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 <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 <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/ <http://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 iterations 38 Newton solver: 1, 1.313e+03 (1.00e-30), 1.283e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 76 Newton solver: 2, 1.198e+02 (1.00e-30), 1.171e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 3, 7.249e-01 (1.00e-30), 7.084e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 4, 3.883e-05 (1.00e-30), 3.795e-08 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 74 Newton solver: 5, 2.778e-12 (1.00e-30), 2.714e-15 (1.00e-10)
16x16 discretization
Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol)
Newton solver: 0, 1.113e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 62 Newton solver: 1, 8.316e+02 (1.00e-30), 7.475e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 141 Newton solver: 2, 5.806e+01 (1.00e-30), 5.218e-02 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 119 Newton solver: 3, 3.309e-01 (1.00e-30), 2.974e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 118 Newton solver: 4, 9.085e-06 (1.00e-30), 8.166e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 120 Newton solver: 5, 3.475e-12 (1.00e-30), 3.124e-15 (1.00e-10)
32x32 discretization
Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol)
Newton solver: 0, 1.330e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 98 Newton solver: 1, 5.913e+02 (1.00e-30), 4.445e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 183 Newton solver: 2, 3.214e+01 (1.00e-30), 2.416e-02 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 152 Newton solver: 3, 2.059e-01 (1.00e-30), 1.547e-04 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 151 Newton solver: 4, 6.949e-06 (1.00e-30), 5.223e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 149 Newton solver: 5, 5.300e-12 (1.00e-30), 3.983e-15 (1.00e-10)
64x64 discretization
Newton solver: iter, abs. residual (abs. tol), rel. residual (rel. tol)
Newton solver: 0, 1.707e+03 (1.00e-30), 1.000e+00 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 198 Newton solver: 1, 4.259e+02 (1.00e-30), 2.494e-01 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 357 Newton solver: 2, 1.706e+01 (1.00e-30), 9.993e-03 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 266 Newton solver: 3, 1.134e-01 (1.00e-30), 6.639e-05 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 261 Newton solver: 4, 4.285e-06 (1.00e-30), 2.510e-09 (1.00e-10) Linear fieldsplit_1_ solve converged due to CONVERGED_RTOL iterations 263 Newton solver: 5, 9.650e-12 (1.00e-30), 5.652e-15 (1.00e-10)
-- Dr. Sebastian Blauth Fraunhofer-Institut für Techno- und Wirtschaftsmathematik ITWM Abteilung Transportvorgänge Fraunhofer-Platz 1, 67663 Kaiserslautern Telefon: +49 631 31600-4968 [email protected] www.itwm.fraunhofer.de
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
Sebastian Blauth <[email protected]> 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.
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 <[email protected]> wrote:
Sebastian Blauth <[email protected]> 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.
Alexander Lindsay <[email protected]> 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.
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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
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 <[email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
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 <[email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <[email protected]> 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 <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> wrote:
Alexander Lindsay <[email protected] <mailto:[email protected]>> 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.
On 23 Jun 2023, at 10:06 PM, Pierre Jolivet <[email protected]> wrote:
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <[email protected]> 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 <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> wrote:
Alexander Lindsay <[email protected] <mailto:[email protected]>> 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.
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 <[email protected]> wrote:
On 23 Jun 2023, at 10:06 PM, Pierre Jolivet <[email protected]> wrote:
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <[email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
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 <[email protected]> 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 <[email protected]> wrote:
On 23 Jun 2023, at 10:06 PM, Pierre Jolivet <[email protected]> wrote:
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <[email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
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 <[email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
On 23 Jun 2023, at 10:06 PM, Pierre Jolivet <[email protected]> wrote:
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <[email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
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 <[email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
On 23 Jun 2023, at 10:06 PM, Pierre Jolivet <[email protected]> wrote:
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <[email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
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 <[email protected]> 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 < [email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
On 23 Jun 2023, at 10:06 PM, Pierre Jolivet <[email protected]> wrote:
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay < [email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
On Thu, Jul 6, 2023 at 8:30 PM Alexander Lindsay <[email protected]> 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_XhuZStZ0AAAA... 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 < [email protected]> 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 < [email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
On 23 Jun 2023, at 10:06 PM, Pierre Jolivet <[email protected]> wrote:
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay < [email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
See src/snes/tutorials/ex70.c for the code that I think was used for that paper. Alexander Lindsay <[email protected]> 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 <[email protected]> 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 <[email protected]> wrote:
On Thu, Jul 6, 2023 at 8:30 PM Alexander Lindsay <[email protected]> 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_XhuZStZ0AAAA...
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 <[email protected]> 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 <[email protected]> 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 <[email protected]> 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 <[email protected]> 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 <[email protected]> wrote:
On 23 Jun 2023, at 10:06 PM, Pierre Jolivet <[email protected]> wrote:
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <[email protected]> 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 <[email protected]> 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 <[email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
-- 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
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 <[email protected]> wrote:
See src/snes/tutorials/ex70.c for the code that I think was used for that paper.
Alexander Lindsay <[email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
On Thu, Jul 6, 2023 at 8:30 PM Alexander Lindsay < [email protected]> 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_XhuZStZ0AAAA...
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 <
[email protected]> 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 < [email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
On 23 Jun 2023, at 10:06 PM, Pierre Jolivet <[email protected]> wrote:
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <[email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
-- 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
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 <[email protected]> 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 <[email protected]> wrote:
See src/snes/tutorials/ex70.c for the code that I think was used for that paper.
Alexander Lindsay <[email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
On Thu, Jul 6, 2023 at 8:30 PM Alexander Lindsay < [email protected]> 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_XhuZStZ0AAAA...
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 <
[email protected]> 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 < [email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
On 23 Jun 2023, at 10:06 PM, Pierre Jolivet <[email protected]> wrote:
On 23 Jun 2023, at 9:39 PM, Alexander Lindsay <[email protected]> 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 < [email protected]> 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 < [email protected]> 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 <[email protected]> wrote:
Alexander Lindsay <[email protected]> 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.
-- 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
participants (7)
-
Alexander Lindsay -
Blauth, Sebastian -
Jed Brown -
Mark Adams -
Matthew Knepley -
Pierre Jolivet -
Sebastian Blauth