feat(par): add convergence check based on petsc L2 norm (MODFLOW-ORG#1970)

* marked this, and petsc preconditioning, as development feature
* modify print summary to reflect new config setting

Author: mjr-deltares

src/Solution/PETSc/PetscConvergence.F90

@@ -3,25 +3,29 @@ module PetscConvergenceModule
   use petscksp
   use KindModule, only: I4B, DP
   use ConstantsModule, only: DPREC, DZERO
+  use SimModule, only: store_error
+  use SimVariablesModule, only: errmsg
   use ListModule
   use ConvergenceSummaryModule
   use ImsLinearSettingsModule
   implicit none
   private
 
-  public :: petsc_check_convergence
+  public :: petsc_cnvg_check
   public :: KSPSetConvergenceTest
 
   ! TODO_MJR: this could be smaller, find a bound
   real(DP), private, parameter :: RNORM_L2_TOL = DPREC
 
-  !< Context for the custom convergence check
+  ! Context for the custom convergence check
   type, public :: PetscCnvgCtxType
     Vec :: x_old !< x vector from the previous iteration
     Vec :: delta_x !< delta in x w.r.t. previous iteration
     Vec :: residual !< the unpreconditoned residual vector (a la IMS)
     integer(I4B) :: icnvg_ims !< IMS convergence number: 1 => converged, -1 => forces next Picard iter
-    integer(I4B) :: icnvgopt !< convergence option from IMS settings
+    integer(I4B) :: icnvgopt !< convergence option:
+                             !! 0,1,2,3,4,.. for equivalent IMS settings,
+                             !! 100,... for PETSc specific settings
     real(DP) :: dvclose !< dep. variable closure criterion
     real(DP) :: rclose !< residual closure criterion
     integer(I4B) :: max_its !< maximum number of inner iterations
@@ -88,9 +92,9 @@ subroutine create(this, mat, settings, summary)
 
   end subroutine create
 
-  !> @brief Routine to check the convergence. This is called
-  !< from within PETSc.
-  subroutine petsc_check_convergence(ksp, n, rnorm_L2, flag, context, ierr)
+  !> @brief Routine to check the convergence following the configuration
+  !< of IMS. (called back from the PETSc solver)
+  subroutine petsc_cnvg_check(ksp, n, rnorm_L2, flag, context, ierr)
     KSP :: ksp !< Iterative context
     PetscInt :: n !< Iteration number
     PetscReal :: rnorm_L2 !< 2-norm (preconditioned) residual value
@@ -99,15 +103,9 @@ subroutine petsc_check_convergence(ksp, n, rnorm_L2, flag, context, ierr)
     PetscErrorCode :: ierr !< error
     ! local
     PetscReal, parameter :: min_one = -1.0
-    PetscReal, dimension(:), pointer :: local_dx, local_res
-    PetscReal :: xnorm_inf_ims, rnorm_inf_ims, rnorm_L2_ims
-    PetscReal :: dvmax_model, rmax_model
-    PetscInt :: idx_dv, idx_r
-    Vec :: x, rhs
-    Mat :: Amat
+    PetscReal :: xnorm_inf, rnorm0, rnorm_inf_ims, rnorm_L2_ims
+    Vec :: x, res
     type(ConvergenceSummaryType), pointer :: summary
-    PetscInt :: iter_cnt
-    PetscInt :: i, j, istart, iend
 
     summary => context%cnvg_summary
 
@@ -116,25 +114,25 @@ subroutine petsc_check_convergence(ksp, n, rnorm_L2, flag, context, ierr)
     call KSPBuildSolution(ksp, PETSC_NULL_VEC, x, ierr)
     CHKERRQ(ierr)
 
-    call KSPGetRhs(ksp, rhs, ierr)
-    CHKERRQ(ierr)
-
-    call KSPGetOperators(ksp, Amat, PETSC_NULL_MAT, ierr)
-    CHKERRQ(ierr)
-
-    call MatMult(Amat, x, context%residual, ierr)
-    CHKERRQ(ierr)
-
-    ! y = x + beta y (i.e. r = b - A*x)
-    call VecAYPX(context%residual, -1.0_DP, rhs, ierr)
+    ! for CG the KSPBuildResidual returns the work vector directly,
+    ! but BCGS (and possibly others) will do the matrix multiplication
+    call KSPBuildResidual(ksp, PETSC_NULL_VEC, PETSC_NULL_VEC, res, ierr)
     CHKERRQ(ierr)
 
-    call VecNorm(context%residual, NORM_2, rnorm_L2_ims, ierr)
-    CHKERRQ(ierr)
+    rnorm0 = huge(rnorm0)
+    if (context%icnvgopt == 2 .or. &
+        context%icnvgopt == 3 .or. &
+        context%icnvgopt == 4) then
+      call VecNorm(res, NORM_2, rnorm_L2_ims, ierr)
+      rnorm0 = rnorm_L2_ims
+      CHKERRQ(ierr)
+    else if (context%icnvgopt == 100) then
+      rnorm0 = rnorm_L2
+    end if
 
     ! n == 0 is before the iteration starts
     if (n == 0) then
-      context%rnorm_L2_init = rnorm_L2_ims
+      context%rnorm_L2_init = rnorm0
       if (rnorm_L2 < RNORM_L2_TOL) then
         ! exact solution found
         flag = KSP_CONVERGED_HAPPY_BREAKDOWN
@@ -144,80 +142,45 @@ subroutine petsc_check_convergence(ksp, n, rnorm_L2, flag, context, ierr)
         flag = KSP_CONVERGED_ITERATING
       end if
       ! early return
+      call VecDestroy(res, ierr)
+      CHKERRQ(ierr)
       return
     end if
 
-    ! increment iteration counter
-    summary%iter_cnt = summary%iter_cnt + 1
-    iter_cnt = summary%iter_cnt
-
-    if (summary%nitermax > 1) then
-      summary%itinner(iter_cnt) = n
-      do i = 1, summary%convnmod
-        summary%convdvmax(i, iter_cnt) = DZERO
-        summary%convlocdv(i, iter_cnt) = 0
-        summary%convrmax(i, iter_cnt) = DZERO
-        summary%convlocr(i, iter_cnt) = 0
-      end do
-    end if
-
     call VecWAXPY(context%delta_x, min_one, context%x_old, x, ierr)
     CHKERRQ(ierr)
 
-    call VecNorm(context%delta_x, NORM_INFINITY, xnorm_inf_ims, ierr)
+    call VecNorm(context%delta_x, NORM_INFINITY, xnorm_inf, ierr)
     CHKERRQ(ierr)
 
-    rnorm_inf_ims = 0.0
+    rnorm_inf_ims = huge(rnorm_inf_ims)
     if (context%icnvgopt == 0 .or. context%icnvgopt == 1) then
-      call VecNorm(context%residual, NORM_INFINITY, rnorm_inf_ims, ierr)
+      call VecNorm(res, NORM_INFINITY, rnorm_inf_ims, ierr)
       CHKERRQ(ierr)
     end if
 
     call VecCopy(x, context%x_old, ierr)
     CHKERRQ(ierr)
 
-    ! get dv and dr per local model (readonly!)
-    call VecGetArrayReadF90(context%delta_x, local_dx, ierr)
-    CHKERRQ(ierr)
-    call VecGetArrayReadF90(context%residual, local_res, ierr)
-    CHKERRQ(ierr)
-    do i = 1, summary%convnmod
-      ! reset
-      dvmax_model = DZERO
-      idx_dv = 0
-      rmax_model = DZERO
-      idx_r = 0
-      ! get first and last model index
-      istart = summary%model_bounds(i)
-      iend = summary%model_bounds(i + 1) - 1
-      do j = istart, iend
-        if (abs(local_dx(j)) > abs(dvmax_model)) then
-          dvmax_model = local_dx(j)
-          idx_dv = j
-        end if
-        if (abs(local_res(j)) > abs(rmax_model)) then
-          rmax_model = local_res(j)
-          idx_r = j
-        end if
-      end do
-      if (summary%nitermax > 1) then
-        summary%convdvmax(i, iter_cnt) = dvmax_model
-        summary%convlocdv(i, iter_cnt) = idx_dv
-        summary%convrmax(i, iter_cnt) = rmax_model
-        summary%convlocr(i, iter_cnt) = idx_r
-      end if
-    end do
-    call VecRestoreArrayF90(context%delta_x, local_dx, ierr)
-    CHKERRQ(ierr)
-    call VecRestoreArrayF90(context%residual, local_res, ierr)
-    CHKERRQ(ierr)
+    ! fill the summary for reporting
+    call fill_cnvg_summary(summary, context%delta_x, res, n)
 
     if (rnorm_L2 < RNORM_L2_TOL) then
       ! exact solution, set to 'converged'
       flag = KSP_CONVERGED_HAPPY_BREAKDOWN
-    else
+    else if (context%icnvgopt < 100) then
       ! IMS check on convergence
-      flag = apply_check(context, n, xnorm_inf_ims, rnorm_inf_ims, rnorm_L2_ims)
+      flag = apply_check(context, n, xnorm_inf, rnorm_inf_ims, rnorm_L2_ims)
+    else if (context%icnvgopt == 100) then
+      ! use PETSc rnorm directly
+      flag = KSP_CONVERGED_ITERATING
+      if (xnorm_inf < context%dvclose .and. rnorm_L2 < context%rclose) then
+        flag = KSP_CONVERGED_HAPPY_BREAKDOWN
+      end if
+    else
+      ! invalid option somehow
+      write (errmsg, '(a,i0)') "Invalid convergence option: ", context%icnvgopt
+      call store_error(errmsg, .true.)
     end if
 
     if (flag == KSP_CONVERGED_ITERATING) then
@@ -227,7 +190,10 @@ subroutine petsc_check_convergence(ksp, n, rnorm_L2, flag, context, ierr)
       end if
     end if
 
-  end subroutine petsc_check_convergence
+    call VecDestroy(res, ierr)
+    CHKERRQ(ierr)
+
+  end subroutine petsc_cnvg_check
 
   !> @brief Apply the IMS convergence check
   !<
@@ -268,6 +234,73 @@ function apply_check(ctx, nit, dvmax, rnorm_inf, rnorm_L2) result(flag)
 
   end function apply_check
 
+  !> @brief Fill the convergence summary from the context
+  !<
+  subroutine fill_cnvg_summary(summary, dx, res, n)
+    type(ConvergenceSummaryType), pointer :: summary !< the convergence summary
+    Vec :: dx !< the vector with changes in x
+    Vec :: res !< the residual vector
+    PetscInt :: n !< the PETSc iteration number
+    ! local
+    PetscReal, dimension(:), pointer :: local_dx, local_res
+    PetscReal :: dvmax_model, rmax_model
+    PetscErrorCode :: ierr
+    PetscInt :: idx_dv, idx_r
+    PetscInt :: i, j, istart, iend
+    PetscInt :: iter_cnt
+
+    ! increment iteration counter
+    summary%iter_cnt = summary%iter_cnt + 1
+    iter_cnt = summary%iter_cnt
+
+    if (summary%nitermax > 1) then
+      summary%itinner(iter_cnt) = n
+      do i = 1, summary%convnmod
+        summary%convdvmax(i, iter_cnt) = DZERO
+        summary%convlocdv(i, iter_cnt) = 0
+        summary%convrmax(i, iter_cnt) = DZERO
+        summary%convlocr(i, iter_cnt) = 0
+      end do
+    end if
+
+    ! get dv and dr per local model (readonly!)
+    call VecGetArrayReadF90(dx, local_dx, ierr)
+    CHKERRQ(ierr)
+    call VecGetArrayReadF90(res, local_res, ierr)
+    CHKERRQ(ierr)
+    do i = 1, summary%convnmod
+      ! reset
+      dvmax_model = DZERO
+      idx_dv = 0
+      rmax_model = DZERO
+      idx_r = 0
+      ! get first and last model index
+      istart = summary%model_bounds(i)
+      iend = summary%model_bounds(i + 1) - 1
+      do j = istart, iend
+        if (abs(local_dx(j)) > abs(dvmax_model)) then
+          dvmax_model = local_dx(j)
+          idx_dv = j
+        end if
+        if (abs(local_res(j)) > abs(rmax_model)) then
+          rmax_model = local_res(j)
+          idx_r = j
+        end if
+      end do
+      if (summary%nitermax > 1) then
+        summary%convdvmax(i, iter_cnt) = dvmax_model
+        summary%convlocdv(i, iter_cnt) = idx_dv
+        summary%convrmax(i, iter_cnt) = rmax_model
+        summary%convlocr(i, iter_cnt) = idx_r
+      end if
+    end do
+    call VecRestoreArrayF90(dx, local_dx, ierr)
+    CHKERRQ(ierr)
+    call VecRestoreArrayF90(res, local_res, ierr)
+    CHKERRQ(ierr)
+
+  end subroutine fill_cnvg_summary
+
   subroutine destroy(this)
     class(PetscCnvgCtxType) :: this
     ! local

src/Solution/PETSc/PetscSolver.F90

@@ -14,6 +14,7 @@ module PetscSolverModule
   use ImsLinearSettingsModule
   use SimVariablesModule, only: iout, simulation_mode
   use SimModule, only: store_error, store_warning
+  use DevFeatureModule, only: dev_feature
 
   implicit none
   private
@@ -28,6 +29,7 @@ module PetscSolverModule
 
     type(ImsLinearSettingsType), pointer :: linear_settings => null() !< pointer to linear settings from IMS
     logical(LGP) :: use_ims_pc !< when true, use custom IMS-style preconditioning
+    logical(LGP) :: use_ims_cnvgopt !< when true, use IMS convergence check in PETSc solve
     KSPType :: ksp_type !< the KSP solver type (CG, BCGS, ...)
     PCType :: pc_type !< the (global) preconditioner type, should be parallel
     PCType :: sub_pc_type !< the block preconditioner type (for the subdomain)
@@ -99,7 +101,8 @@ subroutine petsc_initialize(this, matrix, linear_settings, convergence_summary)
 
     call this%petsc_check_settings(linear_settings)
 
-    this%use_ims_pc = .true. ! default is true, override with .petscrc
+    this%use_ims_cnvgopt = .true. ! use IMS convergence check, override with .petscrc
+    this%use_ims_pc = .true. ! use IMS preconditioning, override with .petscrc
     allocate (this%pc_context)
     call this%pc_context%create(this%matrix, linear_settings)
 
@@ -187,10 +190,19 @@ subroutine get_options_mf6(this)
     ! local
     PetscErrorCode :: ierr
     logical(LGP) :: found
+    logical(LGP) :: use_petsc_pc, use_petsc_cnvg
 
+    use_petsc_pc = .false.
     call PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, &
-                             '-ims_pc', this%use_ims_pc, found, ierr)
+                             '-use_petsc_pc', use_petsc_pc, found, ierr)
     CHKERRQ(ierr)
+    this%use_ims_pc = .not. use_petsc_pc
+
+    use_petsc_cnvg = .false.
+    call PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, &
+                             '-use_petsc_cnvg', use_petsc_cnvg, found, ierr)
+    CHKERRQ(ierr)
+    this%use_ims_cnvgopt = .not. use_petsc_cnvg
 
   end subroutine get_options_mf6
 
@@ -216,6 +228,8 @@ subroutine create_ksp(this)
     if (this%use_ims_pc) then
       call this%set_ims_pc()
     else
+      call dev_feature('PETSc preconditioning is under development, install the &
+            &nightly build or compile from source with IDEVELOPMODE = 1.')
       call this%set_petsc_pc()
     end if
 
@@ -251,10 +265,10 @@ subroutine set_petsc_pc(this)
     CHKERRQ(ierr)
     call PCSetType(sub_pc, this%sub_pc_type, ierr)
     CHKERRQ(ierr)
-    call PCSetFromOptions(sub_pc, ierr)
-    CHKERRQ(ierr)
     call PCFactorSetLevels(sub_pc, this%linear_settings%level, ierr)
     CHKERRQ(ierr)
+    call PCSetFromOptions(sub_pc, ierr)
+    CHKERRQ(ierr)
 
   end subroutine set_petsc_pc
 
@@ -302,8 +316,14 @@ subroutine create_convergence_check(this, convergence_summary)
 
     call this%petsc_ctx%create(this%mat_petsc, this%linear_settings, &
                                convergence_summary)
+    if (.not. this%use_ims_cnvgopt) then
+      ! use PETSc residual L2 norm for convergence
+      call dev_feature('Using PETSc convergence is under development, install &
+      &the nightly build or compile from source with IDEVELOPMODE = 1.')
+      this%petsc_ctx%icnvgopt = 100
+    end if
 
-    call KSPSetConvergenceTest(this%ksp_petsc, petsc_check_convergence, &
+    call KSPSetConvergenceTest(this%ksp_petsc, petsc_cnvg_check, &
                                this%petsc_ctx, PETSC_NULL_FUNCTION, ierr)
     CHKERRQ(ierr)
 
@@ -405,8 +425,13 @@ subroutine petsc_print_summary(this)
       "Dep. var. closure criterion:  ", trim(adjustl(dvclose_str))
     write (iout, '(1x,a,a)') &
       "Residual closure criterion:   ", trim(adjustl(rclose_str))
-    write (iout, '(1x,a,i0)') &
-      "Residual convergence option:  ", this%linear_settings%icnvgopt
+    if (this%use_ims_cnvgopt) then
+      write (iout, '(1x,a,i0)') &
+        "Residual convergence option:  ", this%linear_settings%icnvgopt
+    else
+      write (iout, '(1x,a)') &
+        "Residual convergence option:  PETSc L2 norm"
+    end if
     write (iout, '(1x,a,a)') &
       "Relaxation factor MILU(T):    ", trim(adjustl(relax_str))
     write (iout, '(1x,a,i0)') &