Commits

Anonymous committed 208c3aa

SNES ex70.c: compare iteration counts to custom SIMPLE

From [petsc-maint #155997]:

When I compare my version of SIMPLE with the out-of-the-box
version, I see that the latter requires too many iterations for
the fieldsplit_1 solve because it only uses the jacobi
preconditioner, see below. In order to reproduce the results from
Elman e.a. Taxonomy paper, we will need to
replace "-fieldsplit_1_pc_type jacobi" by something more powerful
such as "-fieldsplit_1_pc_type gamg".

In an attempt to quantify the behaviour I've added a routine
"iterationCount" which was meant to give the total number of iterations
of the "fieldsplit_0_ksp" and "fieldsplit_1_ksp" but it only gives the
number for the last coupled iteration.

my version of SIMPLE:
-----------------------------
$ mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_ksp
iterations coupled system = 57
iterations sub-system 0 = 1
iterations sub-system 1 = 48
residual u = 3.03871e-05
residual p = 1.4044e-07
residual [u,p] = 3.03874e-05
discretization error u = 0.000183699
discretization error p = 0.249002
discretization error [u,p] = 0.249002

Out-of-the box version of SIMPLE:
--------------------------------------------
$ mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_inner_pc_type jacobi -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi
iterations coupled system = 43
iterations sub-system 0 = 1
iterations sub-system 1 = 351
residual u = 3.20371e-05
residual p = 1.5949e-07
residual [u,p] = 3.20375e-05
discretization error u = 0.000184736
discretization error p = 0.248922
discretization error [u,p] = 0.248922

  • Participants
  • Parent commits a050ed3
  • Branches klaij/fieldsplit-simple-ex70

Comments (0)

Files changed (1)

File src/snes/examples/tutorials/ex70.c

 /* Discretized with the cell-centered finite-volume method on a                */
 /* Cartesian grid with co-located variables. Variables ordered as              */
 /* [u1...uN v1...vN p1...pN]^T. Matrix [A00 A01; A10, A11] solved with         */
-/* PCFIELDSPLIT.                                                               */
+/* PCFIELDSPLIT. Lower factorization is used to mimick the Semi-Implicit       */
+/* Method for Pressure Linked Equations (SIMPLE) used as preconditioner        */
+/* instead of solver.                                                          */
 /*                                                                             */
 /* Disclaimer: does not contain the pressure-weighed interpolation             */
 /* method needed to suppress spurious pressure modes in real-life              */
 /*                                                                             */
 /* usage:                                                                      */
 /*                                                                             */
-/* mpiexec -n 2 ./stokes -nx 32 -ny 48                                         */
+/* mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_1_pc_type none
 /*                                                                             */
-/*   Runs with PETSc defaults on 32x48 grid, no PC for the Schur               */
-/*   complement because A11 is zero.                                           */
+/*   Runs with PCFIELDSPLIT on 32x48 grid, no PC for the Schur                 */
+/*   complement because A11 is zero. FGMRES is needed because                  */
+/*   PCFIELDSPLIT is a variable preconditioner.                                */
 /*                                                                             */
-/* mpiexec -n 2 ./stokes -nx 32 -ny 48 -fieldsplit_1_user_pc                   */
+/* mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
 /*                                                                             */
 /*   Same as above but with user defined PC for the true Schur                 */
 /*   complement. PC based on the SIMPLE-type approximation (inverse of         */
 /*   A00 approximated by inverse of its diagonal).                             */
 /*                                                                             */
-/* mpiexec -n 2 ./stokes -nx 32 -ny 48 -fieldsplit_1_user_ksp                  */
+/* mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_ksp
 /*                                                                             */
 /*   Replace the true Schur complement with a user defined Schur               */
 /*   complement based on the SIMPLE-type approximation. Same matrix is         */
 /*   used as PC.                                                               */
 /*                                                                             */
-/* mpiexec -n 2 ./stokes -nx 32 -ny 48 -fieldsplit_1_user_ksp -fieldsplit_1_ksp_rtol 0.01 -fieldsplit_0_ksp_rtol 0.01 */
+/* mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_inner_pc_type jacobi -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi
 /*                                                                             */
-/*   SIMPLE-type approximations are crude, there's no benefit in               */
-/*   solving the subsystems in the preconditioner very accurately.             */
+/*   Out-of-the-box SIMPLE-type preconditioning. The major advantage           */
+/*   is that the user neither needs to provide the approximation of            */
+/*   the Schur complement, nor the corresponding preconditioner.               */
 /*                                                                             */
 /*---------------------------------------------------------------------------- */
 
   PetscFunctionReturn(0);
 }
 
+PetscErrorCode countIterations(KSP ksp)
+{
+  PetscErrorCode ierr;
+  PetscInt       n=1,its;
+  KSP            *subksp;
+  PC             pc;
+
+  PetscFunctionBeginUser;
+  /* nb of iterations of coupled system */
+  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
+  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD," iterations coupled system = %d\n", its);CHKERRQ(ierr);
+  /* TODO: this should somehow become the total of all sub-system
+     solver, not just the total of the last fgmres iteration*/
+  /* nb of iterations of sub-system 0 */
+  ierr = PCFieldSplitGetSubKSP(pc,&n,&subksp);CHKERRQ(ierr);
+  ierr = KSPGetIterationNumber(subksp[0],&its);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD," iterations sub-system 0   = %d\n", its);CHKERRQ(ierr);
+  /* nb of iterations of sub-system 1 */
+  ierr = KSPGetIterationNumber(subksp[1],&its);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD," iterations sub-system 1   = %d\n", its);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
 int main(int argc, char **argv)
 {
   Stokes         s;
   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   ierr = StokesSetupPC(&s, ksp);CHKERRQ(ierr);
   ierr = KSPSolve(ksp, s.b, s.x);CHKERRQ(ierr);
+  ierr = countIterations(ksp);CHKERRQ(ierr);
 
   /* don't trust, verify! */
   ierr = StokesCalcResidual(&s);CHKERRQ(ierr);