MODFLOW 6  version 6.7.0.dev2
USGS Modular Hydrologic Model
petscconvergencemodule Module Reference

Data Types

type  petsccnvgctxtype
 x vector from the previous iteration More...
 
interface  CnvgCheckFunc
 
interface  CnvgDestroyFunc
 
interface  KSPSetConvergenceTest
 

Functions/Subroutines

subroutine create (this, mat, settings, summary)
 
subroutine, public petsc_cnvg_check (ksp, n, rnorm_L2, flag, context, ierr)
 Routine to check the convergence following the configuration. More...
 
function apply_check (ctx, nit, dvmax, rnorm_inf, rnorm_L2)
 Apply the IMS convergence check. More...
 
subroutine, public petsc_cnvg_check_internal (ksp, n, rnorm_L2, flag, context, ierr)
 Routine to check the convergence following the configuration. More...
 
subroutine fill_cnvg_summary (summary, dx, res, n)
 Fill the convergence summary from the context. More...
 
subroutine fill_cnvg_summary_internal (summary, dx, n)
 Fill the convergence summary from the context. More...
 
subroutine destroy (this)
 

Variables

real(dp), parameter, private rnorm_l2_tol = DPREC
 

Function/Subroutine Documentation

◆ apply_check()

function petscconvergencemodule::apply_check ( class(petsccnvgctxtype ctx,
integer(i4b)  nit,
real(dp)  dvmax,
real(dp)  rnorm_inf,
real(dp)  rnorm_L2 
)
private
Parameters
nititeration number
dvmaxinfinity norm of dep. var. change
rnorm_infinfinity norm of residual change
rnorm_l2L2-norm of residual change
rnorm_l2the convergence status

Definition at line 192 of file PetscConvergence.F90.

193  use tdismodule, only: kstp
195  class(PetscCnvgCtxType) :: ctx
196  integer(I4B) :: nit !< iteration number
197  real(DP) :: dvmax !< infinity norm of dep. var. change
198  real(DP) :: rnorm_inf !< infinity norm of residual change
199  real(DP) :: rnorm_L2 !< L2-norm of residual change
200  kspconvergedreason :: flag !< the convergence status
201  ! local
202  real(DP) :: epfact
203  real(DP) :: rcnvg
204 
205  ! Set to 'not converged'
206  flag = ksp_converged_iterating
207  ctx%icnvg_ims = 0
208 
209  epfact = ims_base_epfact(ctx%icnvgopt, kstp)
210 
211  if (ctx%icnvgopt == 2 .or. &
212  ctx%icnvgopt == 3 .or. &
213  ctx%icnvgopt == 4) then
214  rcnvg = rnorm_l2
215  else
216  rcnvg = rnorm_inf
217  end if
218  call ims_base_testcnvg(ctx%icnvgopt, ctx%icnvg_ims, nit, &
219  dvmax, rcnvg, ctx%rnorm_L2_init, &
220  epfact, ctx%dvclose, ctx%rclose)
221 
222  if (ctx%icnvg_ims /= 0) then
223  ! Set to 'converged'
224  flag = ksp_converged_happy_breakdown
225  end if
226 
This module contains the IMS linear accelerator subroutines.
subroutine ims_base_testcnvg(Icnvgopt, Icnvg, Iiter, Dvmax, Rmax, Rmax0, Epfact, Dvclose, Rclose)
@ brief Test for solver convergence
real(dp) function ims_base_epfact(icnvgopt, kstp)
Function returning EPFACT.
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
Here is the call graph for this function:
Here is the caller graph for this function:

◆ create()

subroutine petscconvergencemodule::create ( class(petsccnvgctxtype this,
pointer  mat,
type(imslinearsettingstype), pointer  settings,
type(convergencesummarytype), pointer  summary 
)
private

Definition at line 71 of file PetscConvergence.F90.

72  class(PetscCnvgCtxType) :: this
73  mat, pointer :: mat
74  type(ImsLinearSettingsType), pointer :: settings
75  type(ConvergenceSummaryType), pointer :: summary
76  ! local
77  petscerrorcode :: ierr
78 
79  this%icnvg_ims = 0
80  this%icnvgopt = settings%icnvgopt
81  this%dvclose = settings%dvclose
82  this%rclose = settings%rclose
83  this%max_its = settings%iter1
84  this%cnvg_summary => summary
85  call matcreatevecs(mat, this%x_old, petsc_null_vec, ierr)
86  chkerrq(ierr)
87  call matcreatevecs(mat, this%delta_x, petsc_null_vec, ierr)
88  chkerrq(ierr)
89  call matcreatevecs(mat, this%residual, petsc_null_vec, ierr)
90  chkerrq(ierr)
91 

◆ destroy()

subroutine petscconvergencemodule::destroy ( class(petsccnvgctxtype this)
private

Definition at line 431 of file PetscConvergence.F90.

432  class(PetscCnvgCtxType) :: this
433  ! local
434  integer(I4B) :: ierr
435 
436  call vecdestroy(this%x_old, ierr)
437  chkerrq(ierr)
438  call vecdestroy(this%delta_x, ierr)
439  chkerrq(ierr)
440  call vecdestroy(this%residual, ierr)
441  chkerrq(ierr)
442 

◆ fill_cnvg_summary()

subroutine petscconvergencemodule::fill_cnvg_summary ( type(convergencesummarytype), pointer  summary,
  dx,
  res,
  n 
)
private
Parameters
summarythe convergence summary
summarythe vector with changes in x
summarythe residual vector
summarythe PETSc iteration number

Definition at line 312 of file PetscConvergence.F90.

313  type(ConvergenceSummaryType), pointer :: summary !< the convergence summary
314  vec :: dx !< the vector with changes in x
315  vec :: res !< the residual vector
316  petscint :: n !< the PETSc iteration number
317  ! local
318  petscreal, dimension(:), pointer :: local_dx, local_res
319  petscreal :: dvmax_model, rmax_model
320  petscerrorcode :: ierr
321  petscint :: idx_dv, idx_r
322  petscint :: i, j, istart, iend
323  petscint :: iter_cnt
324 
325  ! increment iteration counter
326  summary%iter_cnt = summary%iter_cnt + 1
327  iter_cnt = summary%iter_cnt
328 
329  if (summary%nitermax > 1) then
330  summary%itinner(iter_cnt) = n
331  do i = 1, summary%convnmod
332  summary%convdvmax(i, iter_cnt) = dzero
333  summary%convlocdv(i, iter_cnt) = 0
334  summary%convrmax(i, iter_cnt) = dzero
335  summary%convlocr(i, iter_cnt) = 0
336  end do
337  end if
338 
339  ! get dv and dr per local model (readonly!)
340  call vecgetarrayreadf90(dx, local_dx, ierr)
341  chkerrq(ierr)
342  call vecgetarrayreadf90(res, local_res, ierr)
343  chkerrq(ierr)
344  do i = 1, summary%convnmod
345  ! reset
346  dvmax_model = dzero
347  idx_dv = 0
348  rmax_model = dzero
349  idx_r = 0
350  ! get first and last model index
351  istart = summary%model_bounds(i)
352  iend = summary%model_bounds(i + 1) - 1
353  do j = istart, iend
354  if (abs(local_dx(j)) > abs(dvmax_model)) then
355  dvmax_model = local_dx(j)
356  idx_dv = j
357  end if
358  if (abs(local_res(j)) > abs(rmax_model)) then
359  rmax_model = local_res(j)
360  idx_r = j
361  end if
362  end do
363  if (summary%nitermax > 1) then
364  summary%convdvmax(i, iter_cnt) = dvmax_model
365  summary%convlocdv(i, iter_cnt) = idx_dv
366  summary%convrmax(i, iter_cnt) = rmax_model
367  summary%convlocr(i, iter_cnt) = idx_r
368  end if
369  end do
370  call vecrestorearrayf90(dx, local_dx, ierr)
371  chkerrq(ierr)
372  call vecrestorearrayf90(res, local_res, ierr)
373  chkerrq(ierr)
374 
Here is the caller graph for this function:

◆ fill_cnvg_summary_internal()

subroutine petscconvergencemodule::fill_cnvg_summary_internal ( type(convergencesummarytype), pointer  summary,
  dx,
  n 
)
private
Parameters
summarythe convergence summary
summarythe vector with changes in x
summarythe PETSc iteration number

Definition at line 379 of file PetscConvergence.F90.

380  type(ConvergenceSummaryType), pointer :: summary !< the convergence summary
381  vec :: dx !< the vector with changes in x
382  petscint :: n !< the PETSc iteration number
383  ! local
384  petscreal, dimension(:), pointer :: local_dx
385  petscreal :: dvmax_model
386  petscerrorcode :: ierr
387  petscint :: idx_dv
388  petscint :: i, j, istart, iend
389  petscint :: iter_cnt
390 
391  ! increment iteration counter
392  summary%iter_cnt = summary%iter_cnt + 1
393  iter_cnt = summary%iter_cnt
394 
395  if (summary%nitermax > 1) then
396  summary%itinner(iter_cnt) = n
397  do i = 1, summary%convnmod
398  summary%convdvmax(i, iter_cnt) = dzero
399  summary%convlocdv(i, iter_cnt) = 0
400  summary%convrmax(i, iter_cnt) = dzero
401  summary%convlocr(i, iter_cnt) = 0
402  end do
403  end if
404 
405  ! get dv and dr per local model (readonly!)
406  call vecgetarrayreadf90(dx, local_dx, ierr)
407  chkerrq(ierr)
408  do i = 1, summary%convnmod
409  ! reset
410  dvmax_model = dzero
411  idx_dv = 0
412  ! get first and last model index
413  istart = summary%model_bounds(i)
414  iend = summary%model_bounds(i + 1) - 1
415  do j = istart, iend
416  if (abs(local_dx(j)) > abs(dvmax_model)) then
417  dvmax_model = local_dx(j)
418  idx_dv = j
419  end if
420  end do
421  if (summary%nitermax > 1) then
422  summary%convdvmax(i, iter_cnt) = dvmax_model
423  summary%convlocdv(i, iter_cnt) = idx_dv
424  end if
425  end do
426  call vecrestorearrayf90(dx, local_dx, ierr)
427  chkerrq(ierr)
428 
Here is the caller graph for this function:

◆ petsc_cnvg_check()

subroutine, public petscconvergencemodule::petsc_cnvg_check (   ksp,
  n,
  rnorm_L2,
  flag,
class(petsccnvgctxtype), pointer  context,
  ierr 
)
Parameters
contexterror

Definition at line 96 of file PetscConvergence.F90.

97  ksp :: ksp !< Iterative context
98  petscint :: n !< Iteration number
99  petscreal :: rnorm_l2 !< 2-norm (preconditioned) residual value
100  kspconvergedreason :: flag !< Converged reason
101  class(PetscCnvgCtxType), pointer :: context !< context
102  petscerrorcode :: ierr !< error
103  ! local
104  petscreal, parameter :: min_one = -1.0
105  petscreal :: xnorm_inf, rnorm0, rnorm_inf_ims, rnorm_l2_ims
106  vec :: x, res
107  type(ConvergenceSummaryType), pointer :: summary
108 
109  summary => context%cnvg_summary
110 
111  ! NB: KSPBuildResidual needs to have its vector destroyed
112  ! to avoid a memory leak, KSPBuildSolution doesn't...
113  call kspbuildsolution(ksp, petsc_null_vec, x, ierr)
114  chkerrq(ierr)
115 
116  ! for CG the KSPBuildResidual returns the work vector directly,
117  ! but BCGS (and possibly others) will do the matrix multiplication
118  call kspbuildresidual(ksp, petsc_null_vec, petsc_null_vec, res, ierr)
119  chkerrq(ierr)
120 
121  rnorm0 = huge(rnorm0)
122  if (context%icnvgopt == 2 .or. &
123  context%icnvgopt == 3 .or. &
124  context%icnvgopt == 4) then
125  call vecnorm(res, norm_2, rnorm_l2_ims, ierr)
126  rnorm0 = rnorm_l2_ims
127  chkerrq(ierr)
128  end if
129 
130  ! n == 0 is before the iteration starts
131  if (n == 0) then
132  context%rnorm_L2_init = rnorm0
133  if (rnorm_l2 < rnorm_l2_tol) then
134  ! exact solution found
135  flag = ksp_converged_happy_breakdown
136  else
137  call veccopy(x, context%x_old, ierr)
138  chkerrq(ierr)
139  flag = ksp_converged_iterating
140  end if
141  ! early return
142  call vecdestroy(res, ierr)
143  chkerrq(ierr)
144  return
145  end if
146 
147  call vecwaxpy(context%delta_x, min_one, context%x_old, x, ierr)
148  chkerrq(ierr)
149 
150  call vecnorm(context%delta_x, norm_infinity, xnorm_inf, ierr)
151  chkerrq(ierr)
152 
153  rnorm_inf_ims = huge(rnorm_inf_ims)
154  if (context%icnvgopt == 0 .or. context%icnvgopt == 1) then
155  call vecnorm(res, norm_infinity, rnorm_inf_ims, ierr)
156  chkerrq(ierr)
157  end if
158 
159  call veccopy(x, context%x_old, ierr)
160  chkerrq(ierr)
161 
162  ! fill the summary for reporting
163  call fill_cnvg_summary(summary, context%delta_x, res, n)
164 
165  if (rnorm_l2 < rnorm_l2_tol) then
166  ! exact solution, set to 'converged'
167  flag = ksp_converged_happy_breakdown
168  context%icnvg_ims = 1
169  ! check for strict option
170  if (n > 1 .and. context%icnvgopt == 1) then
171  context%icnvg_ims = -1
172  end if
173  else
174  ! IMS check on convergence
175  flag = apply_check(context, n, xnorm_inf, rnorm_inf_ims, rnorm_l2_ims)
176  end if
177 
178  if (flag == ksp_converged_iterating) then
179  ! not yet converged, max. iters reached? Then stop.
180  if (n == context%max_its) then
181  flag = ksp_diverged_its
182  end if
183  end if
184 
185  call vecdestroy(res, ierr)
186  chkerrq(ierr)
187 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ petsc_cnvg_check_internal()

subroutine, public petscconvergencemodule::petsc_cnvg_check_internal (   ksp,
  n,
  rnorm_L2,
  flag,
class(petsccnvgctxtype), pointer  context,
  ierr 
)
Parameters
contexterror

Definition at line 231 of file PetscConvergence.F90.

232  ksp :: ksp !< Iterative context
233  petscint :: n !< Iteration number
234  petscreal :: rnorm_l2 !< 2-norm (preconditioned) residual value
235  kspconvergedreason :: flag !< Converged reason
236  class(PetscCnvgCtxType), pointer :: context !< context
237  petscerrorcode :: ierr !< error
238  ! local
239  petscreal, parameter :: min_one = -1.0
240  petscreal :: xnorm_inf, rnorm0
241  vec :: x, res
242  type(ConvergenceSummaryType), pointer :: summary
243 
244  summary => context%cnvg_summary
245 
246  ! NB: KSPBuildResidual needs to have its vector destroyed
247  ! to avoid a memory leak, KSPBuildSolution doesn't...
248  call kspbuildsolution(ksp, petsc_null_vec, x, ierr)
249  chkerrq(ierr)
250 
251  rnorm0 = rnorm_l2
252 
253  ! n == 0 is before the iteration starts
254  if (n == 0) then
255  context%rnorm_L2_init = rnorm0
256  if (rnorm_l2 < rnorm_l2_tol) then
257  ! exact solution found
258  flag = ksp_converged_happy_breakdown
259  context%icnvg_ims = 1
260  else
261  call veccopy(x, context%x_old, ierr)
262  chkerrq(ierr)
263  flag = ksp_converged_iterating
264  end if
265  ! early return
266  return
267  end if
268 
269  call vecwaxpy(context%delta_x, min_one, context%x_old, x, ierr)
270  chkerrq(ierr)
271 
272  call vecnorm(context%delta_x, norm_infinity, xnorm_inf, ierr)
273  chkerrq(ierr)
274 
275  call veccopy(x, context%x_old, ierr)
276  chkerrq(ierr)
277 
278  ! fill the summary for reporting
279  call fill_cnvg_summary_internal(summary, context%delta_x, n)
280 
281  context%icnvg_ims = 0
282  if (rnorm_l2 < rnorm_l2_tol) then
283  ! exact solution, set to 'converged'
284  flag = ksp_converged_happy_breakdown
285  context%icnvg_ims = 1
286  ! apply 'strict'
287  if (n > 1 .and. context%icnvgopt == 1) context%icnvg_ims = -1
288  else
289  ! L2norm
290  if (xnorm_inf <= context%dvclose .and. rnorm_l2 <= context%rclose) then
291  flag = ksp_converged_happy_breakdown
292  context%icnvg_ims = 1
293  ! apply 'strict'
294  if (n > 1 .and. context%icnvgopt == 1) context%icnvg_ims = -1
295  end if
296  end if
297 
298  if (flag == ksp_converged_iterating) then
299  ! not yet converged, max. iters reached? Then stop.
300  if (n == context%max_its) then
301  flag = ksp_diverged_its
302  end if
303  end if
304 
305  call vecdestroy(res, ierr)
306  chkerrq(ierr)
307 
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ rnorm_l2_tol

real(dp), parameter, private petscconvergencemodule::rnorm_l2_tol = DPREC
private

Definition at line 18 of file PetscConvergence.F90.

18  real(DP), private, parameter :: RNORM_L2_TOL = dprec