Skip to content

Commit

Permalink
feat(par): add convergence check based on petsc L2 norm (MODFLOW-USGS…
Browse files Browse the repository at this point in the history 
…#1970)

* marked this, and petsc preconditioning, as development feature
* modify print summary to reflect new config setting
  • Loading branch information
@mjr-deltares
mjr-deltares authored Jul 30, 2024
1 parent a1af64c commit 5ba7681
Show file tree
Hide file tree
Showing 2 changed files with 148 additions and 90 deletions.
199 changes: 116 additions & 83 deletions src/Solution/PETSc/PetscConvergence.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
!<
Expand Down Expand Up @@ -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
Expand Down
39 changes: 32 additions & 7 deletions src/Solution/PETSc/PetscSolver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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

Expand All @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)') &
Expand Down

0 comments on commit 5ba7681

Please sign in to comment.