PetscConvergence.F90

Go to the documentation of this file.

module petscconvergencemodule
#include <petsc/finclude/petscksp.h>
  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_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
  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:
                             !! 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
    real(dp) :: rnorm_l2_init !< the initial L2 norm for (b - Ax)
    type(convergencesummarytype), pointer :: cnvg_summary => null() !< detailed convergence information
    real(dp) :: t_convergence_check !< the time spent convergence checking
  contains
    procedure :: create
    procedure :: destroy
  end type petsccnvgctxtype
 
  ! passing our context into PETSc requires an explicit interface
  ! on KSPSetConvergenceTest, it is defined here:
  interface
    subroutine cnvgcheckfunc(ksp, n, rnorm, flag, context, ierr)
      import tksp, petsccnvgctxtype
      type(tksp) :: ksp
      petscint :: n
      petscreal :: rnorm
      kspconvergedreason :: flag
      class(petsccnvgctxtype), pointer :: context
      petscerrorcode :: ierr
    end subroutine
 
    subroutine cnvgdestroyfunc(context, ierr)
      import petsccnvgctxtype
      class(petsccnvgctxtype), pointer :: context
      petscerrorcode :: ierr
    end subroutine
 
    subroutine kspsetconvergencetest(ksp, check_convergence, context, &
                                     destroy, ierr)
      import tksp, cnvgcheckfunc, petsccnvgctxtype, cnvgdestroyfunc
      type(tksp) :: ksp
      procedure(cnvgcheckfunc) :: check_convergence
      class(petsccnvgctxtype), pointer :: context
      procedure(cnvgdestroyfunc) :: destroy
      petscerrorcode :: ierr
    end subroutine
  end interface
 
contains
 
  subroutine create(this, mat, settings, summary)
    class(petsccnvgctxtype) :: this
    mat, pointer :: mat
    type(imslinearsettingstype), pointer :: settings
    type(convergencesummarytype), pointer :: summary
    ! local
    petscerrorcode :: ierr
 
    this%icnvg_ims = 0
    this%icnvgopt = settings%icnvgopt
    this%dvclose = settings%dvclose
    this%rclose = settings%rclose
    this%max_its = settings%iter1
    this%cnvg_summary => summary
    call matcreatevecs(mat, this%x_old, petsc_null_vec, ierr)
    chkerrq(ierr)
    call matcreatevecs(mat, this%delta_x, petsc_null_vec, ierr)
    chkerrq(ierr)
    call matcreatevecs(mat, this%residual, petsc_null_vec, ierr)
    chkerrq(ierr)
 
  end subroutine create
 
  !> @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
    kspconvergedreason :: flag !< Converged reason
    class(petsccnvgctxtype), pointer :: context !< context
    petscerrorcode :: ierr !< error
    ! local
    petscreal, parameter :: min_one = -1.0
    petscreal :: xnorm_inf, rnorm0, rnorm_inf_ims, rnorm_l2_ims
    vec :: x, res
    type(convergencesummarytype), pointer :: summary
 
    summary => context%cnvg_summary
 
    ! NB: KSPBuildResidual needs to have its vector destroyed
    ! to avoid a memory leak, KSPBuildSolution doesn't...
    call kspbuildsolution(ksp, petsc_null_vec, x, ierr)
    chkerrq(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)
 
    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 = rnorm0
      if (rnorm_l2 < rnorm_l2_tol) then
        ! exact solution found
        flag = ksp_converged_happy_breakdown
      else
        call veccopy(x, context%x_old, ierr)
        chkerrq(ierr)
        flag = ksp_converged_iterating
      end if
      ! early return
      call vecdestroy(res, ierr)
      chkerrq(ierr)
      return
    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, ierr)
    chkerrq(ierr)
 
    rnorm_inf_ims = huge(rnorm_inf_ims)
    if (context%icnvgopt == 0 .or. context%icnvgopt == 1) then
      call vecnorm(res, norm_infinity, rnorm_inf_ims, ierr)
      chkerrq(ierr)
    end if
 
    call veccopy(x, context%x_old, 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 if (context%icnvgopt < 100) then
      ! IMS check on convergence
      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
      ! not yet converged, max. iters reached? Then stop.
      if (n == context%max_its) then
        flag = ksp_diverged_its
      end if
    end if
 
    call vecdestroy(res, ierr)
    chkerrq(ierr)
 
  end subroutine petsc_cnvg_check
 
  !> @brief Apply the IMS convergence check
  !<
  function apply_check(ctx, nit, dvmax, rnorm_inf, rnorm_L2) result(flag)
    use tdismodule, only: kstp
    use imslinearbasemodule, only: ims_base_testcnvg, ims_base_epfact
    class(petsccnvgctxtype) :: ctx
    integer(I4B) :: nit !< iteration number
    real(dp) :: dvmax !< infinity norm of dep. var. change
    real(dp) :: rnorm_inf !< infinity norm of residual change
    real(dp) :: rnorm_l2 !< L2-norm of residual change
    kspconvergedreason :: flag !< the convergence status
    ! local
    real(dp) :: epfact
    real(dp) :: rcnvg
 
    ! Set to 'not converged'
    flag = ksp_converged_iterating
    ctx%icnvg_ims = 0
 
    epfact = ims_base_epfact(ctx%icnvgopt, kstp)
 
    if (ctx%icnvgopt == 2 .or. &
        ctx%icnvgopt == 3 .or. &
        ctx%icnvgopt == 4) then
      rcnvg = rnorm_l2
    else
      rcnvg = rnorm_inf
    end if
    call ims_base_testcnvg(ctx%icnvgopt, ctx%icnvg_ims, nit, &
                           dvmax, rcnvg, ctx%rnorm_L2_init, &
                           epfact, ctx%dvclose, ctx%rclose)
 
    if (ctx%icnvg_ims /= 0) then
      ! Set to 'converged'
      flag = ksp_converged_happy_breakdown
    end if
 
  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
    integer(I4B) :: ierr
 
    call vecdestroy(this%x_old, ierr)
    chkerrq(ierr)
    call vecdestroy(this%delta_x, ierr)
    chkerrq(ierr)
    call vecdestroy(this%residual, ierr)
    chkerrq(ierr)
 
  end subroutine destroy
 
end module petscconvergencemodule