MODFLOW 6  version 6.5.0.dev2
MODFLOW 6 Code Documentation
Functions/Subroutines | Variables
imslinearbasemodule Module Reference

This module contains the IMS linear accelerator subroutines. More...

Functions/Subroutines

subroutine ims_base_cg (ICNVG, ITMAX, INNERIT, NEQ, NJA, NIAPC, NJAPC, IPC, ICNVGOPT, NORTH, DVCLOSE, RCLOSE, L2NORM0, EPFACT, IA0, JA0, A0, IAPC, JAPC, APC, X, B, D, P, Q, Z, NJLU, IW, JLU, NCONV, CONVNMOD, CONVMODSTART, CACCEL, summary)
 @ brief Preconditioned Conjugate Gradient linear accelerator More...
 
subroutine ims_base_bcgs (ICNVG, ITMAX, INNERIT, NEQ, NJA, NIAPC, NJAPC, IPC, ICNVGOPT, NORTH, ISCL, DSCALE, DVCLOSE, RCLOSE, L2NORM0, EPFACT, IA0, JA0, A0, IAPC, JAPC, APC, X, B, D, P, Q, T, V, DHAT, PHAT, QHAT, NJLU, IW, JLU, NCONV, CONVNMOD, CONVMODSTART, CACCEL, summary)
 @ brief Preconditioned BiConjugate Gradient Stabilized linear accelerator More...
 
subroutine ims_base_calc_order (IORD, NEQ, NJA, IA, JA, LORDER, IORDER)
 @ brief Calculate LORDER AND IORDER More...
 
subroutine ims_base_scale (IOPT, ISCL, NEQ, NJA, IA, JA, AMAT, X, B, DSCALE, DSCALE2)
 @ brief Scale the coefficient matrix More...
 
subroutine ims_base_pcu (IOUT, NJA, NEQ, NIAPC, NJAPC, IPC, RELAX, AMAT, IA, JA, APC, IAPC, JAPC, IW, W, LEVEL, DROPTOL, NJLU, NJW, NWLU, JLU, JW, WLU)
 @ brief Update the preconditioner More...
 
subroutine ims_base_pcjac (NJA, NEQ, AMAT, APC, IA, JA)
 @ brief Jacobi preconditioner More...
 
subroutine ims_base_jaca (NEQ, A, D1, D2)
 @ brief Apply the Jacobi preconditioner More...
 
subroutine ims_base_pcilu0 (NJA, NEQ, AMAT, IA, JA, APC, IAPC, JAPC, IW, W, RELAX, IPCFLAG, DELTA)
 @ brief Update the ILU0 preconditioner More...
 
subroutine ims_base_ilu0a (NJA, NEQ, APC, IAPC, JAPC, R, D)
 @ brief Apply the ILU0 and MILU0 preconditioners More...
 
subroutine ims_base_testcnvg (Icnvgopt, Icnvg, Iiter, Dvmax, Rmax, Rmax0, Epfact, Dvclose, Rclose)
 @ brief Test for solver convergence More...
 
subroutine ims_calc_pcdims (neq, nja, ia, level, ipc, niapc, njapc, njlu, njw, nwlu)
 
subroutine ims_base_pccrs (NEQ, NJA, IA, JA, IAPC, JAPC)
 @ brief Generate CRS pointers for the preconditioner More...
 
subroutine ims_base_isort (NVAL, IARRAY)
 In-place sorting for an integer array. More...
 
subroutine ims_base_residual (NEQ, NJA, X, B, D, A, IA, JA)
 Calculate residual. More...
 
real(dp) function ims_base_epfact (icnvgopt, kstp)
 Function returning EPFACT. More...
 

Variables

type(blockparsertype), private parser
 

Detailed Description

This module contains the IMS linear accelerator subroutines used by a MODFLOW 6 solution.

Function/Subroutine Documentation

◆ ims_base_bcgs()

subroutine imslinearbasemodule::ims_base_bcgs ( integer(i4b), intent(inout)  ICNVG,
integer(i4b), intent(in)  ITMAX,
integer(i4b), intent(inout)  INNERIT,
integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NIAPC,
integer(i4b), intent(in)  NJAPC,
integer(i4b), intent(in)  IPC,
integer(i4b), intent(in)  ICNVGOPT,
integer(i4b), intent(in)  NORTH,
integer(i4b), intent(in)  ISCL,
real(dp), dimension(neq), intent(in)  DSCALE,
real(dp), intent(in)  DVCLOSE,
real(dp), intent(in)  RCLOSE,
real(dp), intent(in)  L2NORM0,
real(dp), intent(in)  EPFACT,
integer(i4b), dimension(neq + 1), intent(in)  IA0,
integer(i4b), dimension(nja), intent(in)  JA0,
real(dp), dimension(nja), intent(in)  A0,
integer(i4b), dimension(niapc + 1), intent(in)  IAPC,
integer(i4b), dimension(njapc), intent(in)  JAPC,
real(dp), dimension(njapc), intent(in)  APC,
real(dp), dimension(neq), intent(inout)  X,
real(dp), dimension(neq), intent(in)  B,
real(dp), dimension(neq), intent(inout)  D,
real(dp), dimension(neq), intent(inout)  P,
real(dp), dimension(neq), intent(inout)  Q,
real(dp), dimension(neq), intent(inout)  T,
real(dp), dimension(neq), intent(inout)  V,
real(dp), dimension(neq), intent(inout)  DHAT,
real(dp), dimension(neq), intent(inout)  PHAT,
real(dp), dimension(neq), intent(inout)  QHAT,
integer(i4b), intent(in)  NJLU,
integer(i4b), dimension(niapc), intent(in)  IW,
integer(i4b), dimension(njlu), intent(in)  JLU,
integer(i4b), intent(in)  NCONV,
integer(i4b), intent(in)  CONVNMOD,
integer(i4b), dimension(convnmod + 1), intent(inout)  CONVMODSTART,
character(len=31), dimension(nconv), intent(inout)  CACCEL,
type(convergencesummarytype), intent(in), pointer  summary 
)

Apply the Preconditioned BiConjugate Gradient Stabilized linear accelerator to the current coefficient matrix, right-hand side, using the currentdependent-variable.

Parameters
[in,out]icnvgconvergence flag (1) non-convergence (0)
[in]itmaxmaximum number of inner iterations
[in,out]inneritinner iteration count
[in]neqnumber of equations
[in]njanumber of non-zero entries
[in]niapcpreconditioner number of rows
[in]njapcpreconditioner number of non-zero entries
[in]ipcpreconditioner option
[in]icnvgoptflow convergence criteria option
[in]northorthogonalization frequency
[in]isclscaling option
[in]dscalescaling vector
[in]dvclosedependent-variable closure criteria
[in]rcloseflow closure criteria
[in]l2norm0initial L-2 norm for system of equations
[in]epfactfactor for decreasing flow convergence criteria for subsequent Picard iterations
[in]ia0CRS row pointers
[in]ja0CRS column pointers
[in]a0coefficient matrix
[in]iapcpreconditioner CRS row pointers
[in]japcpreconditioner CRS column pointers
[in]apcpreconditioner matrix
[in,out]xdependent-variable vector
[in]bright-hand side vector
[in,out]dpreconditioner working vector
[in,out]ppreconditioner working vector
[in,out]qpreconditioner working vector
[in,out]tpreconditioner working vector
[in,out]vpreconditioner working vector
[in,out]dhatBCGS preconditioner working vector
[in,out]phatBCGS preconditioner working vector
[in,out]qhatBCGS preconditioner working vector
[in]njlupreconditioner length of JLU vector
[in]iwpreconditioner integer working vector
[in]jlupreconditioner JLU working vector
[in]nconvmaximum number of inner iterations in a time step (maxiter * maxinner)
[in]convnmodnumber of models in the solution
[in,out]convmodstartpointer to the start of each model in the convmod* arrays
[in,out]caccelconvergence string
[in]summaryConvergence summary report

Definition at line 252 of file ImsLinearBase.f90.

262  ! -- dummy variables
263  integer(I4B), INTENT(INOUT) :: ICNVG !< convergence flag (1) non-convergence (0)
264  integer(I4B), INTENT(IN) :: ITMAX !< maximum number of inner iterations
265  integer(I4B), INTENT(INOUT) :: INNERIT !< inner iteration count
266  integer(I4B), INTENT(IN) :: NEQ !< number of equations
267  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
268  integer(I4B), INTENT(IN) :: NIAPC !< preconditioner number of rows
269  integer(I4B), INTENT(IN) :: NJAPC !< preconditioner number of non-zero entries
270  integer(I4B), INTENT(IN) :: IPC !< preconditioner option
271  integer(I4B), INTENT(IN) :: ICNVGOPT !< flow convergence criteria option
272  integer(I4B), INTENT(IN) :: NORTH !< orthogonalization frequency
273  integer(I4B), INTENT(IN) :: ISCL !< scaling option
274  real(DP), DIMENSION(NEQ), INTENT(IN) :: DSCALE !< scaling vector
275  real(DP), INTENT(IN) :: DVCLOSE !< dependent-variable closure criteria
276  real(DP), INTENT(IN) :: RCLOSE !< flow closure criteria
277  real(DP), INTENT(IN) :: L2NORM0 !< initial L-2 norm for system of equations
278  real(DP), INTENT(IN) :: EPFACT !< factor for decreasing flow convergence criteria for subsequent Picard iterations
279  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA0 !< CRS row pointers
280  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA0 !< CRS column pointers
281  real(DP), DIMENSION(NJA), INTENT(IN) :: A0 !< coefficient matrix
282  integer(I4B), DIMENSION(NIAPC + 1), INTENT(IN) :: IAPC !< preconditioner CRS row pointers
283  integer(I4B), DIMENSION(NJAPC), INTENT(IN) :: JAPC !< preconditioner CRS column pointers
284  real(DP), DIMENSION(NJAPC), INTENT(IN) :: APC !< preconditioner matrix
285  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: X !< dependent-variable vector
286  real(DP), DIMENSION(NEQ), INTENT(IN) :: B !< right-hand side vector
287  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: D !< preconditioner working vector
288  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: P !< preconditioner working vector
289  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: Q !< preconditioner working vector
290  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: T !< preconditioner working vector
291  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: V !< preconditioner working vector
292  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: DHAT !< BCGS preconditioner working vector
293  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: PHAT !< BCGS preconditioner working vector
294  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: QHAT !< BCGS preconditioner working vector
295  ! -- ILUT dummy variables
296  integer(I4B), INTENT(IN) :: NJLU !< preconditioner length of JLU vector
297  integer(I4B), DIMENSION(NIAPC), INTENT(IN) :: IW !< preconditioner integer working vector
298  integer(I4B), DIMENSION(NJLU), INTENT(IN) :: JLU !< preconditioner JLU working vector
299  ! -- convergence information dummy variables
300  integer(I4B), INTENT(IN) :: NCONV !< maximum number of inner iterations in a time step (maxiter * maxinner)
301  integer(I4B), INTENT(IN) :: CONVNMOD !< number of models in the solution
302  integer(I4B), DIMENSION(CONVNMOD + 1), INTENT(INOUT) :: CONVMODSTART !< pointer to the start of each model in the convmod* arrays
303  character(len=31), DIMENSION(NCONV), INTENT(INOUT) :: CACCEL !< convergence string
304  type(ConvergenceSummaryType), pointer, intent(in) :: summary !< Convergence summary report
305  ! -- local variables
306  LOGICAL :: LORTH
307  logical :: lsame
308  character(len=15) :: cval1, cval2
309  integer(I4B) :: n
310  integer(I4B) :: iiter
311  integer(I4B) :: xloc, rloc
312  integer(I4B) :: im, im0, im1
313  real(DP) :: ddot
314  real(DP) :: tv
315  real(DP) :: deltax
316  real(DP) :: rmax
317  real(DP) :: l2norm
318  real(DP) :: rcnvg
319  real(DP) :: alpha, alpha0
320  real(DP) :: beta
321  real(DP) :: rho, rho0
322  real(DP) :: omega, omega0
323  real(DP) :: numerator, denominator
324  !
325  ! -- initialize local variables
326  innerit = 0
327  alpha = dzero
328  alpha0 = dzero
329  beta = dzero
330  rho = dzero
331  rho0 = dzero
332  omega = dzero
333  omega0 = dzero
334  !
335  ! -- SAVE INITIAL RESIDUAL
336  DO n = 1, neq
337  dhat(n) = d(n)
338  END DO
339  !
340  ! -- INNER ITERATION
341  inner: DO iiter = 1, itmax
342  innerit = innerit + 1
343  summary%iter_cnt = summary%iter_cnt + 1
344  !
345  ! -- CALCULATE rho
346  rho = ddot(neq, dhat, 1, d, 1)
347  !
348  ! -- COMPUTE DIRECTIONAL VECTORS
349  IF (iiter == 1) THEN
350  DO n = 1, neq
351  p(n) = d(n)
352  END DO
353  ELSE
354  beta = (rho / rho0) * (alpha0 / omega0)
355  DO n = 1, neq
356  p(n) = d(n) + beta * (p(n) - omega0 * v(n))
357  END DO
358  END IF
359  !
360  ! -- APPLY PRECONDITIONER TO UPDATE PHAT
361  SELECT CASE (ipc)
362  !
363  ! -- ILU0 AND MILU0
364  CASE (1, 2)
365  CALL ims_base_ilu0a(nja, neq, apc, iapc, japc, p, phat)
366  !
367  ! -- ILUT AND MILUT
368  CASE (3, 4)
369  CALL lusol(neq, p, phat, apc, jlu, iw)
370  END SELECT
371  !
372  ! -- COMPUTE ITERATES
373  !
374  ! -- UPDATE V WITH A AND PHAT
375  call amux(neq, phat, v, a0, ja0, ia0)
376  !
377  ! -- UPDATE alpha WITH DHAT AND V
378  denominator = ddot(neq, dhat, 1, v, 1)
379  denominator = denominator + sign(dprec, denominator)
380  alpha = rho / denominator
381  !
382  ! -- UPDATE Q
383  DO n = 1, neq
384  q(n) = d(n) - alpha * v(n)
385  END DO
386  !
387  ! ! -- CALCULATE INFINITY NORM OF Q - TEST FOR TERMINATION
388  ! ! TERMINATE IF rmax IS LESS THAN MACHINE PRECISION (DPREC)
389  ! rmax = DZERO
390  ! DO n = 1, NEQ
391  ! tv = Q(n)
392  ! IF (ISCL.NE.0 ) tv = tv / DSCALE(n)
393  ! IF (ABS(tv) > ABS(rmax) ) rmax = tv
394  ! END DO
395  ! IF (ABS(rmax).LE.DPREC) THEN
396  ! deltax = DZERO
397  ! DO n = 1, NEQ
398  ! tv = alpha * PHAT(n)
399  ! IF (ISCL.NE.0) THEN
400  ! tv = tv * DSCALE(n)
401  ! END IF
402  ! X(n) = X(n) + tv
403  ! IF (ABS(tv) > ABS(deltax) ) deltax = tv
404  ! END DO
405  ! CALL IMSLINEARSUB_TESTCNVG(ICNVGOPT, ICNVG, INNERIT, &
406  ! deltax, rmax, &
407  ! rmax, EPFACT, DVCLOSE, RCLOSE )
408  ! IF (ICNVG.NE.0 ) EXIT INNER
409  ! END IF
410  !
411  ! -- APPLY PRECONDITIONER TO UPDATE QHAT
412  SELECT CASE (ipc)
413  !
414  ! -- ILU0 AND MILU0
415  CASE (1, 2)
416  CALL ims_base_ilu0a(nja, neq, apc, iapc, japc, q, qhat)
417  !
418  ! -- ILUT AND MILUT
419  CASE (3, 4)
420  CALL lusol(neq, q, qhat, apc, jlu, iw)
421  END SELECT
422  !
423  ! -- UPDATE T WITH A AND QHAT
424  call amux(neq, qhat, t, a0, ja0, ia0)
425  !
426  ! -- UPDATE omega
427  numerator = ddot(neq, t, 1, q, 1)
428  denominator = ddot(neq, t, 1, t, 1)
429  denominator = denominator + sign(dprec, denominator)
430  omega = numerator / denominator
431  !
432  ! -- UPDATE X AND RESIDUAL
433  deltax = dzero
434  rmax = dzero
435  l2norm = dzero
436  DO im = 1, convnmod
437  summary%dvmax(im) = dzero
438  summary%rmax(im) = dzero
439  END DO
440  im = 1
441  im0 = convmodstart(1)
442  im1 = convmodstart(2)
443  DO n = 1, neq
444  !
445  ! -- determine current model index
446  if (n == im1) then
447  im = im + 1
448  im0 = convmodstart(im)
449  im1 = convmodstart(im + 1)
450  end if
451  !
452  ! -- X AND DX
453  tv = alpha * phat(n) + omega * qhat(n)
454  x(n) = x(n) + tv
455  IF (iscl .NE. 0) THEN
456  tv = tv * dscale(n)
457  END IF
458  IF (abs(tv) > abs(deltax)) THEN
459  deltax = tv
460  xloc = n
461  END IF
462  IF (abs(tv) > abs(summary%dvmax(im))) THEN
463  summary%dvmax(im) = tv
464  summary%locdv(im) = n
465  END IF
466  !
467  ! -- RESIDUAL
468  tv = q(n) - omega * t(n)
469  d(n) = tv
470  IF (iscl .NE. 0) THEN
471  tv = tv / dscale(n)
472  END IF
473  IF (abs(tv) > abs(rmax)) THEN
474  rmax = tv
475  rloc = n
476  END IF
477  IF (abs(tv) > abs(summary%rmax(im))) THEN
478  summary%rmax(im) = tv
479  summary%locr(im) = n
480  END IF
481  l2norm = l2norm + tv * tv
482  END DO
483  l2norm = sqrt(l2norm)
484  !
485  ! -- SAVE SOLVER convergence information dummy variables
486  IF (nconv > 1) THEN !<
487  n = summary%iter_cnt
488  WRITE (cval1, '(g15.7)') alpha
489  WRITE (cval2, '(g15.7)') omega
490  caccel(n) = trim(adjustl(cval1))//','//trim(adjustl(cval2))
491  summary%itinner(n) = iiter
492  DO im = 1, convnmod
493  summary%convdvmax(im, n) = summary%dvmax(im)
494  summary%convlocdv(im, n) = summary%locdv(im)
495  summary%convrmax(im, n) = summary%rmax(im)
496  summary%convlocr(im, n) = summary%locr(im)
497  END DO
498  END IF
499  !
500  ! -- TEST FOR SOLVER CONVERGENCE
501  IF (icnvgopt == 2 .OR. icnvgopt == 3 .OR. icnvgopt == 4) THEN
502  rcnvg = l2norm
503  ELSE
504  rcnvg = rmax
505  END IF
506  CALL ims_base_testcnvg(icnvgopt, icnvg, innerit, &
507  deltax, rcnvg, &
508  l2norm0, epfact, dvclose, rclose)
509  !
510  ! -- CHECK FOR EXACT SOLUTION
511  IF (rcnvg == dzero) icnvg = 1
512  !
513  ! -- CHECK FOR STANDARD CONVERGENCE
514  IF (icnvg .NE. 0) EXIT inner
515  !
516  ! -- CHECK THAT CURRENT AND PREVIOUS rho, alpha, AND omega ARE
517  ! DIFFERENT
518  lsame = is_close(rho, rho0)
519  IF (lsame) THEN
520  EXIT inner
521  END IF
522  lsame = is_close(alpha, alpha0)
523  IF (lsame) THEN
524  EXIT inner
525  END IF
526  lsame = is_close(omega, omega0)
527  IF (lsame) THEN
528  EXIT inner
529  END IF
530  !
531  ! -- RECALCULATE THE RESIDUAL
532  IF (north > 0) THEN
533  lorth = mod(iiter + 1, north) == 0
534  IF (lorth) THEN
535  call ims_base_residual(neq, nja, x, b, d, a0, ia0, ja0)
536  END IF
537  END IF
538  !
539  ! -- exit inner if rho or omega are zero
540  if (rho * omega == dzero) then
541  exit inner
542  end if
543  !
544  ! -- SAVE CURRENT INNER ITERATES
545  rho0 = rho
546  alpha0 = alpha
547  omega0 = omega
548  END DO inner
549  !
550  ! -- RESET ICNVG
551  IF (icnvg < 0) icnvg = 0
552  !
553  ! -- RETURN
554  RETURN
ddot
real(kind=8) function ddot(n, dx, incx, dy, incy)
Definition: blas1_d.f90:296
lusol
subroutine lusol(n, y, x, alu, jlu, ju)
Definition: ilut.f90:466
amux
subroutine amux(n, x, y, a, ja, ia)
Definition: sparsekit.f90:2
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_calc_order()

subroutine imslinearbasemodule::ims_base_calc_order ( integer(i4b), intent(in)  IORD,
integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA,
integer(i4b), dimension(neq), intent(inout)  LORDER,
integer(i4b), dimension(neq), intent(inout)  IORDER 
)

Calculate LORDER and IORDER for reordering.

Parameters
[in]iordreordering option
[in]neqnumber of rows
[in]njanumber of non-zero entries
[in]iarow pointer
[in]jacolumn pointer
[in,out]lorderreorder vector
[in,out]iorderinverse of reorder vector

Definition at line 562 of file ImsLinearBase.f90.

563  ! -- modules
564  use simmodule, only: store_error, count_errors
565  ! -- dummy variables
566  integer(I4B), INTENT(IN) :: IORD !< reordering option
567  integer(I4B), INTENT(IN) :: NEQ !< number of rows
568  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
569  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< row pointer
570  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< column pointer
571  integer(I4B), DIMENSION(NEQ), INTENT(INOUT) :: LORDER !< reorder vector
572  integer(I4B), DIMENSION(NEQ), INTENT(INOUT) :: IORDER !< inverse of reorder vector
573  ! -- local variables
574  character(len=LINELENGTH) :: errmsg
575  integer(I4B) :: n
576  integer(I4B) :: nsp
577  integer(I4B), DIMENSION(:), ALLOCATABLE :: iwork0
578  integer(I4B), DIMENSION(:), ALLOCATABLE :: iwork1
579  integer(I4B) :: iflag
580  !
581  ! -- initialize lorder and iorder
582  DO n = 1, neq
583  lorder(n) = izero
584  iorder(n) = izero
585  END DO
586  ! ALLOCATE (iwork0(NEQ))
587  SELECT CASE (iord)
588  CASE (1)
589  CALL genrcm(neq, nja, ia, ja, lorder)
590  CASE (2)
591  nsp = 3 * neq + 4 * nja
592  allocate (iwork0(neq))
593  allocate (iwork1(nsp))
594  CALL ims_odrv(neq, nja, nsp, ia, ja, lorder, iwork0, &
595  iwork1, iflag)
596  IF (iflag .NE. 0) THEN
597  write (errmsg, '(A,1X,A)') &
598  'IMSLINEARSUB_CALC_ORDER error creating minimum degree ', &
599  'order permutation '
600  call store_error(errmsg)
601  END IF
602  !
603  ! -- DEALLOCATE TEMPORARY STORAGE
604  deallocate (iwork0, iwork1)
605  END SELECT
606  !
607  ! -- GENERATE INVERSE OF LORDER
608  DO n = 1, neq
609  iorder(lorder(n)) = n
610  END DO
611  !
612  ! -- terminate if errors occurred
613  if (count_errors() > 0) then
614  call parser%StoreErrorUnit()
615  end if
616  !
617  ! -- RETURN
618  RETURN
simmodule
This module contains simulation methods.
Definition: Sim.f90:10
simmodule::store_error
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
simmodule::count_errors
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
genrcm
subroutine genrcm(node_num, adj_num, adj_row, adj, perm)
Definition: rcm.f90:1004
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_cg()

subroutine imslinearbasemodule::ims_base_cg ( integer(i4b), intent(inout)  ICNVG,
integer(i4b), intent(in)  ITMAX,
integer(i4b), intent(inout)  INNERIT,
integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NIAPC,
integer(i4b), intent(in)  NJAPC,
integer(i4b), intent(in)  IPC,
integer(i4b), intent(in)  ICNVGOPT,
integer(i4b), intent(in)  NORTH,
real(dp), intent(in)  DVCLOSE,
real(dp), intent(in)  RCLOSE,
real(dp), intent(in)  L2NORM0,
real(dp), intent(in)  EPFACT,
integer(i4b), dimension(neq + 1), intent(in)  IA0,
integer(i4b), dimension(nja), intent(in)  JA0,
real(dp), dimension(nja), intent(in)  A0,
integer(i4b), dimension(niapc + 1), intent(in)  IAPC,
integer(i4b), dimension(njapc), intent(in)  JAPC,
real(dp), dimension(njapc), intent(in)  APC,
real(dp), dimension(neq), intent(inout)  X,
real(dp), dimension(neq), intent(inout)  B,
real(dp), dimension(neq), intent(inout)  D,
real(dp), dimension(neq), intent(inout)  P,
real(dp), dimension(neq), intent(inout)  Q,
real(dp), dimension(neq), intent(inout)  Z,
integer(i4b), intent(in)  NJLU,
integer(i4b), dimension(niapc), intent(in)  IW,
integer(i4b), dimension(njlu), intent(in)  JLU,
integer(i4b), intent(in)  NCONV,
integer(i4b), intent(in)  CONVNMOD,
integer(i4b), dimension(convnmod + 1), intent(inout)  CONVMODSTART,
character(len=31), dimension(nconv), intent(inout)  CACCEL,
type(convergencesummarytype), intent(in), pointer  summary 
)

Apply the Preconditioned Conjugate Gradient linear accelerator to the current coefficient matrix, right-hand side, using the current dependent-variable.

Parameters
[in,out]icnvgconvergence flag (1) non-convergence (0)
[in]itmaxmaximum number of inner iterations
[in,out]inneritinner iteration count
[in]neqnumber of equations
[in]njanumber of non-zero entries
[in]niapcpreconditioner number of rows
[in]njapcpreconditioner number of non-zero entries
[in]ipcpreconditioner option
[in]icnvgoptflow convergence criteria option
[in]northorthogonalization frequency
[in]dvclosedependent-variable closure criteria
[in]rcloseflow closure criteria
[in]l2norm0initial L-2 norm for system of equations
[in]epfactfactor for decreasing flow convergence criteria for subsequent Picard iterations
[in]ia0CRS row pointers
[in]ja0CRS column pointers
[in]a0coefficient matrix
[in]iapcpreconditioner CRS row pointers
[in]japcpreconditioner CRS column pointers
[in]apcpreconditioner matrix
[in,out]xdependent-variable vector
[in,out]bright-hand side vector
[in,out]dworking vector
[in,out]pworking vector
[in,out]qworking vector
[in,out]zworking vector
[in]njlupreconditioner length of JLU vector
[in]iwpreconditioner integer working vector
[in]jlupreconditioner JLU working vector
[in]nconvmaximum number of inner iterations in a time step (maxiter * maxinner)
[in]convnmodnumber of models in the solution
[in,out]convmodstartpointer to the start of each model in the convmod* arrays
[in,out]caccelconvergence string
[in]summaryConvergence summary report

Definition at line 30 of file ImsLinearBase.f90.

39  ! -- dummy variables
40  integer(I4B), INTENT(INOUT) :: ICNVG !< convergence flag (1) non-convergence (0)
41  integer(I4B), INTENT(IN) :: ITMAX !< maximum number of inner iterations
42  integer(I4B), INTENT(INOUT) :: INNERIT !< inner iteration count
43  integer(I4B), INTENT(IN) :: NEQ !< number of equations
44  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
45  integer(I4B), INTENT(IN) :: NIAPC !< preconditioner number of rows
46  integer(I4B), INTENT(IN) :: NJAPC !< preconditioner number of non-zero entries
47  integer(I4B), INTENT(IN) :: IPC !< preconditioner option
48  integer(I4B), INTENT(IN) :: ICNVGOPT !< flow convergence criteria option
49  integer(I4B), INTENT(IN) :: NORTH !< orthogonalization frequency
50  real(DP), INTENT(IN) :: DVCLOSE !< dependent-variable closure criteria
51  real(DP), INTENT(IN) :: RCLOSE !< flow closure criteria
52  real(DP), INTENT(IN) :: L2NORM0 !< initial L-2 norm for system of equations
53  real(DP), INTENT(IN) :: EPFACT !< factor for decreasing flow convergence criteria for subsequent Picard iterations
54  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA0 !< CRS row pointers
55  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA0 !< CRS column pointers
56  real(DP), DIMENSION(NJA), INTENT(IN) :: A0 !< coefficient matrix
57  integer(I4B), DIMENSION(NIAPC + 1), INTENT(IN) :: IAPC !< preconditioner CRS row pointers
58  integer(I4B), DIMENSION(NJAPC), INTENT(IN) :: JAPC !< preconditioner CRS column pointers
59  real(DP), DIMENSION(NJAPC), INTENT(IN) :: APC !< preconditioner matrix
60  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: X !< dependent-variable vector
61  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: B !< right-hand side vector
62  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: D !< working vector
63  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: P !< working vector
64  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: Q !< working vector
65  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: Z !< working vector
66  ! -- ILUT dummy variables
67  integer(I4B), INTENT(IN) :: NJLU !< preconditioner length of JLU vector
68  integer(I4B), DIMENSION(NIAPC), INTENT(IN) :: IW !< preconditioner integer working vector
69  integer(I4B), DIMENSION(NJLU), INTENT(IN) :: JLU !< preconditioner JLU working vector
70  ! -- convergence information dummy variables dummy variables
71  integer(I4B), INTENT(IN) :: NCONV !< maximum number of inner iterations in a time step (maxiter * maxinner)
72  integer(I4B), INTENT(IN) :: CONVNMOD !< number of models in the solution
73  integer(I4B), DIMENSION(CONVNMOD + 1), INTENT(INOUT) :: CONVMODSTART !< pointer to the start of each model in the convmod* arrays
74  character(len=31), DIMENSION(NCONV), INTENT(INOUT) :: CACCEL !< convergence string
75  type(ConvergenceSummaryType), pointer, intent(in) :: summary !< Convergence summary report
76  ! -- local variables
77  LOGICAL :: lorth
78  logical :: lsame
79  character(len=31) :: cval
80  integer(I4B) :: n
81  integer(I4B) :: iiter
82  integer(I4B) :: xloc, rloc
83  integer(I4B) :: im, im0, im1
84  real(DP) :: ddot
85  real(DP) :: tv
86  real(DP) :: deltax
87  real(DP) :: rmax
88  real(DP) :: l2norm
89  real(DP) :: rcnvg
90  real(DP) :: denominator
91  real(DP) :: alpha, beta
92  real(DP) :: rho, rho0
93  !
94  ! -- initialize local variables
95  rho0 = dzero
96  rho = dzero
97  innerit = 0
98  !
99  ! -- INNER ITERATION
100  inner: DO iiter = 1, itmax
101  innerit = innerit + 1
102  summary%iter_cnt = summary%iter_cnt + 1
103  !
104  ! -- APPLY PRECONDITIONER
105  SELECT CASE (ipc)
106  !
107  ! -- ILU0 AND MILU0
108  CASE (1, 2)
109  CALL ims_base_ilu0a(nja, neq, apc, iapc, japc, d, z)
110  !
111  ! -- ILUT AND MILUT
112  CASE (3, 4)
113  CALL lusol(neq, d, z, apc, jlu, iw)
114  END SELECT
115  rho = ddot(neq, d, 1, z, 1)
116  !
117  ! -- COMPUTE DIRECTIONAL VECTORS
118  IF (iiter == 1) THEN
119  DO n = 1, neq
120  p(n) = z(n)
121  END DO
122  ELSE
123  beta = rho / rho0
124  DO n = 1, neq
125  p(n) = z(n) + beta * p(n)
126  END DO
127  END IF
128  !
129  ! -- COMPUTE ITERATES
130  !
131  ! -- UPDATE Q
132  call amux(neq, p, q, a0, ja0, ia0)
133  denominator = ddot(neq, p, 1, q, 1)
134  denominator = denominator + sign(dprec, denominator)
135  alpha = rho / denominator
136  !
137  ! -- UPDATE X AND RESIDUAL
138  deltax = dzero
139  rmax = dzero
140  l2norm = dzero
141  DO im = 1, convnmod
142  summary%locdv(im) = 0
143  summary%dvmax(im) = dzero
144  summary%locr(im) = 0
145  summary%rmax(im) = dzero
146  END DO
147  im = 1
148  im0 = convmodstart(1)
149  im1 = convmodstart(2)
150  DO n = 1, neq
151  !
152  ! -- determine current model index
153  if (n == im1) then
154  im = im + 1
155  im0 = convmodstart(im)
156  im1 = convmodstart(im + 1)
157  end if
158  !
159  ! -- identify deltax and rmax
160  tv = alpha * p(n)
161  x(n) = x(n) + tv
162  IF (abs(tv) > abs(deltax)) THEN
163  deltax = tv
164  xloc = n
165  END IF
166  IF (abs(tv) > abs(summary%dvmax(im))) THEN
167  summary%dvmax(im) = tv
168  summary%locdv(im) = n
169  END IF
170  tv = d(n)
171  tv = tv - alpha * q(n)
172  d(n) = tv
173  IF (abs(tv) > abs(rmax)) THEN
174  rmax = tv
175  rloc = n
176  END IF
177  IF (abs(tv) > abs(summary%rmax(im))) THEN
178  summary%rmax(im) = tv
179  summary%locr(im) = n
180  END IF
181  l2norm = l2norm + tv * tv
182  END DO
183  l2norm = sqrt(l2norm)
184  !
185  ! -- SAVE SOLVER convergence information dummy variables
186  IF (nconv > 1) THEN !<
187  n = summary%iter_cnt
188  WRITE (cval, '(g15.7)') alpha
189  caccel(n) = cval
190  summary%itinner(n) = iiter
191  DO im = 1, convnmod
192  summary%convlocdv(im, n) = summary%locdv(im)
193  summary%convlocr(im, n) = summary%locr(im)
194  summary%convdvmax(im, n) = summary%dvmax(im)
195  summary%convrmax(im, n) = summary%rmax(im)
196  END DO
197  END IF
198  !
199  ! -- TEST FOR SOLVER CONVERGENCE
200  IF (icnvgopt == 2 .OR. icnvgopt == 3 .OR. icnvgopt == 4) THEN
201  rcnvg = l2norm
202  ELSE
203  rcnvg = rmax
204  END IF
205  CALL ims_base_testcnvg(icnvgopt, icnvg, innerit, &
206  deltax, rcnvg, &
207  l2norm0, epfact, dvclose, rclose)
208  !
209  ! -- CHECK FOR EXACT SOLUTION
210  IF (rcnvg == dzero) icnvg = 1
211  !
212  ! -- CHECK FOR STANDARD CONVERGENCE
213  IF (icnvg .NE. 0) EXIT inner
214  !
215  ! -- CHECK THAT CURRENT AND PREVIOUS rho ARE DIFFERENT
216  lsame = is_close(rho, rho0)
217  IF (lsame) THEN
218  EXIT inner
219  END IF
220  !
221  ! -- RECALCULATE THE RESIDUAL
222  IF (north > 0) THEN
223  lorth = mod(iiter + 1, north) == 0
224  IF (lorth) THEN
225  call ims_base_residual(neq, nja, x, b, d, a0, ia0, ja0)
226  END IF
227  END IF
228  !
229  ! -- exit inner if rho is zero
230  if (rho == dzero) then
231  exit inner
232  end if
233  !
234  ! -- SAVE CURRENT INNER ITERATES
235  rho0 = rho
236  END DO inner
237  !
238  ! -- RESET ICNVG
239  IF (icnvg < 0) icnvg = 0
240  !
241  ! -- RETURN
242  RETURN
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_epfact()

real(dp) function imslinearbasemodule::ims_base_epfact ( integer(i4b)  icnvgopt,
integer(i4b)  kstp 
)
Parameters
icnvgoptIMS convergence option
kstptime step number
Returns
factor for decreasing convergence criteria in subsequent Picard iterations

Definition at line 1355 of file ImsLinearBase.f90.

1356  integer(I4B) :: icnvgopt !< IMS convergence option
1357  integer(I4B) :: kstp !< time step number
1358  real(DP) :: epfact !< factor for decreasing convergence criteria in subsequent Picard iterations
1359 
1360  if (icnvgopt == 2) then
1361  if (kstp == 1) then
1362  epfact = 0.01
1363  else
1364  epfact = 0.10
1365  end if
1366  else if (icnvgopt == 4) then
1367  epfact = dem4
1368  else
1369  epfact = done
1370  end if
1371 
Here is the caller graph for this function:

◆ ims_base_ilu0a()

subroutine imslinearbasemodule::ims_base_ilu0a ( integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NEQ,
real(dp), dimension(nja), intent(in)  APC,
integer(i4b), dimension(neq + 1), intent(in)  IAPC,
integer(i4b), dimension(nja), intent(in)  JAPC,
real(dp), dimension(neq), intent(in)  R,
real(dp), dimension(neq), intent(inout)  D 
)

Apply the ILU0 and MILU0 preconditioners to the passed vector (R).

Parameters
[in]njanumber of non-zero entries
[in]neqnumber of equations
[in]apcILU0/MILU0 preconditioner matrix
[in]iapcILU0/MILU0 preconditioner CRS row pointers
[in]japcILU0/MILU0 preconditioner CRS column pointers
[in]rinput vector
[in,out]doutput vector after applying APC to R

Definition at line 1073 of file ImsLinearBase.f90.

1074  ! -- dummy variables
1075  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
1076  integer(I4B), INTENT(IN) :: NEQ !< number of equations
1077  real(DP), DIMENSION(NJA), INTENT(IN) :: APC !< ILU0/MILU0 preconditioner matrix
1078  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IAPC !< ILU0/MILU0 preconditioner CRS row pointers
1079  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JAPC !< ILU0/MILU0 preconditioner CRS column pointers
1080  real(DP), DIMENSION(NEQ), INTENT(IN) :: R !< input vector
1081  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: D !< output vector after applying APC to R
1082  ! -- local variables
1083  integer(I4B) :: ic0, ic1
1084  integer(I4B) :: iu
1085  integer(I4B) :: jcol
1086  integer(I4B) :: j, n
1087  real(DP) :: tv
1088  !
1089  ! -- FORWARD SOLVE - APC * D = R
1090  forward: DO n = 1, neq
1091  tv = r(n)
1092  ic0 = iapc(n)
1093  ic1 = iapc(n + 1) - 1
1094  iu = japc(n) - 1
1095  lower: DO j = ic0, iu
1096  jcol = japc(j)
1097  tv = tv - apc(j) * d(jcol)
1098  END DO lower
1099  d(n) = tv
1100  END DO forward
1101  !
1102  ! -- BACKWARD SOLVE - D = D / U
1103  backward: DO n = neq, 1, -1
1104  ic0 = iapc(n)
1105  ic1 = iapc(n + 1) - 1
1106  iu = japc(n)
1107  tv = d(n)
1108  upper: DO j = iu, ic1
1109  jcol = japc(j)
1110  tv = tv - apc(j) * d(jcol)
1111  END DO upper
1112  !
1113  ! -- COMPUTE D FOR DIAGONAL - D = D / U
1114  d(n) = tv * apc(n)
1115  END DO backward
1116  !
1117  ! -- RETURN
1118  RETURN
Here is the caller graph for this function:

◆ ims_base_isort()

subroutine imslinearbasemodule::ims_base_isort ( integer(i4b), intent(in)  NVAL,
integer(i4b), dimension(nval), intent(inout)  IARRAY 
)

Subroutine sort an integer array in-place.

Parameters
[in]nvallength of the integer array
[in,out]iarrayinteger array to be sorted

Definition at line 1301 of file ImsLinearBase.f90.

1302  ! -- dummy variables
1303  integer(I4B), INTENT(IN) :: NVAL !< length of the integer array
1304  integer(I4B), DIMENSION(NVAL), INTENT(INOUT) :: IARRAY !< integer array to be sorted
1305  ! -- local variables
1306  integer(I4B) :: i, j, itemp
1307  ! -- code
1308  DO i = 1, nval - 1
1309  DO j = i + 1, nval
1310  if (iarray(i) > iarray(j)) then
1311  itemp = iarray(j)
1312  iarray(j) = iarray(i)
1313  iarray(i) = itemp
1314  END IF
1315  END DO
1316  END DO
1317  !
1318  ! -- RETURN
1319  RETURN
Here is the caller graph for this function:

◆ ims_base_jaca()

subroutine imslinearbasemodule::ims_base_jaca ( integer(i4b), intent(in)  NEQ,
real(dp), dimension(neq), intent(in)  A,
real(dp), dimension(neq), intent(in)  D1,
real(dp), dimension(neq), intent(inout)  D2 
)

Apply the Jacobi preconditioner and return the resultant vector.

Parameters
[in]neqnumber of equations
[in]aJacobi preconditioner
[in]d1input vector
[in,out]d2resultant vector

Definition at line 925 of file ImsLinearBase.f90.

926  ! -- dummy variables
927  integer(I4B), INTENT(IN) :: NEQ !< number of equations
928  real(DP), DIMENSION(NEQ), INTENT(IN) :: A !< Jacobi preconditioner
929  real(DP), DIMENSION(NEQ), INTENT(IN) :: D1 !< input vector
930  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: D2 !< resultant vector
931  ! -- local variables
932  integer(I4B) :: n
933  real(DP) :: tv
934  ! -- code
935  DO n = 1, neq
936  tv = a(n) * d1(n)
937  d2(n) = tv
938  END DO
939  !
940  ! -- RETURN
941  RETURN

◆ ims_base_pccrs()

subroutine imslinearbasemodule::ims_base_pccrs ( integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA,
integer(i4b), dimension(neq + 1), intent(inout)  IAPC,
integer(i4b), dimension(nja), intent(inout)  JAPC 
)

Generate the CRS row and column pointers for the preconditioner. JAPC(1:NEQ) hHas the position of the upper entry for a row, JAPC(NEQ+1:NJA) is the column position for entry, APC(1:NEQ) is the preconditioned inverse of the diagonal, and APC(NEQ+1:NJA) are the preconditioned entries for off diagonals.

Definition at line 1237 of file ImsLinearBase.f90.

1239  ! -- dummy variables
1240  integer(I4B), INTENT(IN) :: NEQ !<
1241  integer(I4B), INTENT(IN) :: NJA !<
1242  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !<
1243  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !<
1244  integer(I4B), DIMENSION(NEQ + 1), INTENT(INOUT) :: IAPC !<
1245  integer(I4B), DIMENSION(NJA), INTENT(INOUT) :: JAPC !<
1246  ! -- local variables
1247  integer(I4B) :: n, j
1248  integer(I4B) :: i0, i1
1249  integer(I4B) :: nlen
1250  integer(I4B) :: ic, ip
1251  integer(I4B) :: jcol
1252  integer(I4B), DIMENSION(:), ALLOCATABLE :: iarr
1253  ! -- code
1254  ip = neq + 1
1255  DO n = 1, neq
1256  i0 = ia(n)
1257  i1 = ia(n + 1) - 1
1258  nlen = i1 - i0
1259  ALLOCATE (iarr(nlen))
1260  ic = 0
1261  DO j = i0, i1
1262  jcol = ja(j)
1263  IF (jcol == n) cycle
1264  ic = ic + 1
1265  iarr(ic) = jcol
1266  END DO
1267  CALL ims_base_isort(nlen, iarr)
1268  iapc(n) = ip
1269  DO j = 1, nlen
1270  jcol = iarr(j)
1271  japc(ip) = jcol
1272  ip = ip + 1
1273  END DO
1274  DEALLOCATE (iarr)
1275  END DO
1276  iapc(neq + 1) = nja + 1
1277  !
1278  ! -- POSITION OF THE FIRST UPPER ENTRY FOR ROW
1279  DO n = 1, neq
1280  i0 = iapc(n)
1281  i1 = iapc(n + 1) - 1
1282  japc(n) = iapc(n + 1)
1283  DO j = i0, i1
1284  jcol = japc(j)
1285  IF (jcol > n) THEN
1286  japc(n) = j
1287  EXIT
1288  END IF
1289  END DO
1290  END DO
1291  !
1292  ! -- RETURN
1293  RETURN
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_pcilu0()

subroutine imslinearbasemodule::ims_base_pcilu0 ( integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NEQ,
real(dp), dimension(nja), intent(in)  AMAT,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA,
real(dp), dimension(nja), intent(inout)  APC,
integer(i4b), dimension(neq + 1), intent(inout)  IAPC,
integer(i4b), dimension(nja), intent(inout)  JAPC,
integer(i4b), dimension(neq), intent(inout)  IW,
real(dp), dimension(neq), intent(inout)  W,
real(dp), intent(in)  RELAX,
integer(i4b), intent(inout)  IPCFLAG,
real(dp), intent(in)  DELTA 
)

Update the ILU0 preconditioner using the current coefficient matrix.

Parameters
[in]njanumber of non-zero entries
[in]neqnumber of equations
[in]amatcoefficient matrix
[in]iaCRS row pointers
[in]jaCRS column pointers
[in,out]apcpreconditioned matrix
[in,out]iapcpreconditioner CRS row pointers
[in,out]japcpreconditioner CRS column pointers
[in,out]iwpreconditioner integer work vector
[in,out]wpreconditioner work vector
[in]relaxMILU0 preconditioner relaxation factor
[in,out]ipcflagpreconditioner error flag
[in]deltafactor used to correct non-diagonally dominant matrices

Definition at line 949 of file ImsLinearBase.f90.

952  ! -- dummy variables
953  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
954  integer(I4B), INTENT(IN) :: NEQ !< number of equations
955  real(DP), DIMENSION(NJA), INTENT(IN) :: AMAT !< coefficient matrix
956  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< CRS row pointers
957  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< CRS column pointers
958  real(DP), DIMENSION(NJA), INTENT(INOUT) :: APC !< preconditioned matrix
959  integer(I4B), DIMENSION(NEQ + 1), INTENT(INOUT) :: IAPC !< preconditioner CRS row pointers
960  integer(I4B), DIMENSION(NJA), INTENT(INOUT) :: JAPC !< preconditioner CRS column pointers
961  integer(I4B), DIMENSION(NEQ), INTENT(INOUT) :: IW !< preconditioner integer work vector
962  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: W !< preconditioner work vector
963  real(DP), INTENT(IN) :: RELAX !< MILU0 preconditioner relaxation factor
964  integer(I4B), INTENT(INOUT) :: IPCFLAG !< preconditioner error flag
965  real(DP), INTENT(IN) :: DELTA !< factor used to correct non-diagonally dominant matrices
966  ! -- local variables
967  integer(I4B) :: ic0, ic1
968  integer(I4B) :: iic0, iic1
969  integer(I4B) :: iu, iiu
970  integer(I4B) :: j, n
971  integer(I4B) :: jj
972  integer(I4B) :: jcol, jw
973  integer(I4B) :: jjcol
974  real(DP) :: drelax
975  real(DP) :: sd1
976  real(DP) :: tl
977  real(DP) :: rs
978  real(DP) :: d
979  !
980  ! -- initialize local variables
981  drelax = relax
982  DO n = 1, neq
983  iw(n) = 0
984  w(n) = dzero
985  END DO
986  main: DO n = 1, neq
987  ic0 = ia(n)
988  ic1 = ia(n + 1) - 1
989  DO j = ic0, ic1
990  jcol = ja(j)
991  iw(jcol) = 1
992  w(jcol) = w(jcol) + amat(j)
993  END DO
994  ic0 = iapc(n)
995  ic1 = iapc(n + 1) - 1
996  iu = japc(n)
997  rs = dzero
998  lower: DO j = ic0, iu - 1
999  jcol = japc(j)
1000  iic0 = iapc(jcol)
1001  iic1 = iapc(jcol + 1) - 1
1002  iiu = japc(jcol)
1003  tl = w(jcol) * apc(jcol)
1004  w(jcol) = tl
1005  DO jj = iiu, iic1
1006  jjcol = japc(jj)
1007  jw = iw(jjcol)
1008  IF (jw .NE. 0) THEN
1009  w(jjcol) = w(jjcol) - tl * apc(jj)
1010  ELSE
1011  rs = rs + tl * apc(jj)
1012  END IF
1013  END DO
1014  END DO lower
1015  !
1016  ! -- DIAGONAL - CALCULATE INVERSE OF DIAGONAL FOR SOLUTION
1017  d = w(n)
1018  tl = (done + delta) * d - (drelax * rs)
1019  !
1020  ! -- ENSURE THAT THE SIGN OF THE DIAGONAL HAS NOT CHANGED AND IS
1021  sd1 = sign(d, tl)
1022  IF (sd1 .NE. d) THEN
1023  !
1024  ! -- USE SMALL VALUE IF DIAGONAL SCALING IS NOT EFFECTIVE FOR
1025  ! PIVOTS THAT CHANGE THE SIGN OF THE DIAGONAL
1026  IF (ipcflag > 1) THEN
1027  tl = sign(dem6, d)
1028  !
1029  ! -- DIAGONAL SCALING CONTINUES TO BE EFFECTIVE
1030  ELSE
1031  ipcflag = 1
1032  EXIT main
1033  END IF
1034  END IF
1035  IF (abs(tl) == dzero) THEN
1036  !
1037  ! -- USE SMALL VALUE IF DIAGONAL SCALING IS NOT EFFECTIVE FOR
1038  ! ZERO PIVOTS
1039  IF (ipcflag > 1) THEN
1040  tl = sign(dem6, d)
1041  !
1042  ! -- DIAGONAL SCALING CONTINUES TO BE EFFECTIVE FOR ELIMINATING
1043  ELSE
1044  ipcflag = 1
1045  EXIT main
1046  END IF
1047  END IF
1048  apc(n) = done / tl
1049  !
1050  ! -- RESET POINTER FOR IW TO ZERO
1051  iw(n) = 0
1052  w(n) = dzero
1053  DO j = ic0, ic1
1054  jcol = japc(j)
1055  apc(j) = w(jcol)
1056  iw(jcol) = 0
1057  w(jcol) = dzero
1058  END DO
1059  END DO main
1060  !
1061  ! -- RESET IPCFLAG IF SUCCESSFUL COMPLETION OF MAIN
1062  ipcflag = 0
1063  !
1064  ! -- RETURN
1065  RETURN
Here is the caller graph for this function:

◆ ims_base_pcjac()

subroutine imslinearbasemodule::ims_base_pcjac ( integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NEQ,
real(dp), dimension(nja), intent(in)  AMAT,
real(dp), dimension(neq), intent(inout)  APC,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA 
)

Calculate the Jacobi preconditioner (inverse of the diagonal) using the current coefficient matrix.

Parameters
[in]njanumber of non-zero entries
[in]neqnumber of equations
[in]amatcoefficient matrix
[in,out]apcpreconditioner matrix
[in]iaCRS row pointers
[in]jaCRS column pointers

Definition at line 887 of file ImsLinearBase.f90.

888  ! -- dummy variables
889  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
890  integer(I4B), INTENT(IN) :: NEQ !< number of equations
891  real(DP), DIMENSION(NJA), INTENT(IN) :: AMAT !< coefficient matrix
892  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: APC !< preconditioner matrix
893  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< CRS row pointers
894  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< CRS column pointers
895  ! -- local variables
896  integer(I4B) :: i, n
897  integer(I4B) :: ic0, ic1
898  integer(I4B) :: id
899  real(DP) :: tv
900  ! -- code
901  DO n = 1, neq
902  ic0 = ia(n)
903  ic1 = ia(n + 1) - 1
904  id = ia(n)
905  DO i = ic0, ic1
906  IF (ja(i) == n) THEN
907  id = i
908  EXIT
909  END IF
910  END DO
911  tv = amat(id)
912  IF (abs(tv) > dzero) tv = done / tv
913  apc(n) = tv
914  END DO
915  !
916  ! -- RETURN
917  RETURN

◆ ims_base_pcu()

subroutine imslinearbasemodule::ims_base_pcu ( integer(i4b), intent(in)  IOUT,
integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NIAPC,
integer(i4b), intent(in)  NJAPC,
integer(i4b), intent(in)  IPC,
real(dp), intent(in)  RELAX,
real(dp), dimension(nja), intent(in)  AMAT,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA,
real(dp), dimension(njapc), intent(inout)  APC,
integer(i4b), dimension(niapc + 1), intent(inout)  IAPC,
integer(i4b), dimension(njapc), intent(inout)  JAPC,
integer(i4b), dimension(niapc), intent(inout)  IW,
real(dp), dimension(niapc), intent(inout)  W,
integer(i4b), intent(in)  LEVEL,
real(dp), intent(in)  DROPTOL,
integer(i4b), intent(in)  NJLU,
integer(i4b), intent(in)  NJW,
integer(i4b), intent(in)  NWLU,
integer(i4b), dimension(njlu), intent(inout)  JLU,
integer(i4b), dimension(njw), intent(inout)  JW,
real(dp), dimension(nwlu), intent(inout)  WLU 
)

Update the preconditioner using the current coefficient matrix.

Parameters
[in]ioutsimulation listing file unit
[in]njanumber of non-zero entries
[in]neqnumber of equations
[in]niapcpreconditioner number of rows
[in]njapcpreconditioner number of non-zero entries
[in]ipcprecoditioner (1) ILU0 (2) MILU0 (3) ILUT (4) MILUT
[in]relaxpreconditioner relaxation factor for MILU0 and MILUT
[in]amatcoefficient matrix
[in]iaCRS row pointers
[in]jaCRS column pointers
[in,out]apcpreconditioner matrix
[in,out]iapcpreconditioner CRS row pointers
[in,out]japcpreconditioner CRS column pointers
[in,out]iwpreconditioner integed work vector
[in,out]wpreconditioner work vector
[in]levelnumber of levels of fill for ILUT and MILUT
[in]droptoldrop tolerance
[in]njlulength of JLU working vector
[in]njwlength of JW working vector
[in]nwlulength of WLU working vector
[in,out]jluILUT/MILUT JLU working vector
[in,out]jwILUT/MILUT JW working vector
[in,out]wluILUT/MILUT WLU working vector

Definition at line 773 of file ImsLinearBase.f90.

776  ! -- modules
777  use simmodule, only: store_error, count_errors
778  ! -- dummy variables
779  integer(I4B), INTENT(IN) :: IOUT !< simulation listing file unit
780  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
781  integer(I4B), INTENT(IN) :: NEQ !< number of equations
782  integer(I4B), INTENT(IN) :: NIAPC !< preconditioner number of rows
783  integer(I4B), INTENT(IN) :: NJAPC !< preconditioner number of non-zero entries
784  integer(I4B), INTENT(IN) :: IPC !< precoditioner (1) ILU0 (2) MILU0 (3) ILUT (4) MILUT
785  real(DP), INTENT(IN) :: RELAX !< preconditioner relaxation factor for MILU0 and MILUT
786  real(DP), DIMENSION(NJA), INTENT(IN) :: AMAT !< coefficient matrix
787  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< CRS row pointers
788  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< CRS column pointers
789  real(DP), DIMENSION(NJAPC), INTENT(INOUT) :: APC !< preconditioner matrix
790  integer(I4B), DIMENSION(NIAPC + 1), INTENT(INOUT) :: IAPC !< preconditioner CRS row pointers
791  integer(I4B), DIMENSION(NJAPC), INTENT(INOUT) :: JAPC !< preconditioner CRS column pointers
792  integer(I4B), DIMENSION(NIAPC), INTENT(INOUT) :: IW !< preconditioner integed work vector
793  real(DP), DIMENSION(NIAPC), INTENT(INOUT) :: W !< preconditioner work vector
794  ! -- ILUT dummy variables
795  integer(I4B), INTENT(IN) :: LEVEL !< number of levels of fill for ILUT and MILUT
796  real(DP), INTENT(IN) :: DROPTOL !< drop tolerance
797  integer(I4B), INTENT(IN) :: NJLU !< length of JLU working vector
798  integer(I4B), INTENT(IN) :: NJW !< length of JW working vector
799  integer(I4B), INTENT(IN) :: NWLU !< length of WLU working vector
800  integer(I4B), DIMENSION(NJLU), INTENT(INOUT) :: JLU !< ILUT/MILUT JLU working vector
801  integer(I4B), DIMENSION(NJW), INTENT(INOUT) :: JW !< ILUT/MILUT JW working vector
802  real(DP), DIMENSION(NWLU), INTENT(INOUT) :: WLU !< ILUT/MILUT WLU working vector
803  ! -- local variables
804  character(len=LINELENGTH) :: errmsg
805  character(len=100), dimension(5), parameter :: cerr = &
806  ["Elimination process has generated a row in L or U whose length is > n.", &
807  &"The matrix L overflows the array al. ", &
808  &"The matrix U overflows the array alu. ", &
809  &"Illegal value for lfil. ", &
810  &"Zero row encountered. "]
811  integer(i4b) :: ipcflag
812  integer(I4B) :: icount
813  integer(I4B) :: ierr
814  real(DP) :: delta
815  ! -- formats
816 2000 FORMAT(/, ' MATRIX IS SEVERELY NON-DIAGONALLY DOMINANT.', &
817  /, ' ADDED SMALL VALUE TO PIVOT ', i0, ' TIMES IN', &
818  ' IMSLINEARSUB_PCU.')
819  !
820  ! -- initialize local variables
821  ipcflag = 0
822  icount = 0
823  delta = dzero
824  pcscale: DO
825  SELECT CASE (ipc)
826  !
827  ! -- ILU0 AND MILU0
828  CASE (1, 2)
829  CALL ims_base_pcilu0(nja, neq, amat, ia, ja, &
830  apc, iapc, japc, iw, w, &
831  relax, ipcflag, delta)
832  !
833  ! -- ILUT AND MILUT
834  CASE (3, 4)
835  ierr = 0
836  CALL ilut(neq, amat, ja, ia, level, droptol, &
837  apc, jlu, iw, njapc, wlu, jw, ierr, &
838  relax, ipcflag, delta)
839  if (ierr /= 0) then
840  if (ierr > 0) then
841  write (errmsg, '(a,1x,i0,1x,a)') &
842  'ILUT: zero pivot encountered at step number', ierr, '.'
843  else
844  write (errmsg, '(a,1x,a)') 'ILUT:', cerr(-ierr)
845  end if
846  call store_error(errmsg)
847  call parser%StoreErrorUnit()
848  end if
849  !
850  ! -- ADDITIONAL PRECONDITIONERS
851  CASE DEFAULT
852  ipcflag = 0
853  END SELECT
854  IF (ipcflag < 1) THEN
855  EXIT pcscale
856  END IF
857  delta = 1.5d0 * delta + dem3
858  ipcflag = 0
859  IF (delta > dhalf) THEN
860  delta = dhalf
861  ipcflag = 2
862  END IF
863  icount = icount + 1
864  !
865  ! -- terminate pcscale loop if not making progress
866  if (icount > 10) then
867  exit pcscale
868  end if
869 
870  END DO pcscale
871  !
872  ! -- write error message if small value added to pivot
873  if (icount > 0) then
874  write (iout, 2000) icount
875  end if
876  !
877  ! -- RETURN
878  RETURN
ilut
subroutine ilut(n, a, ja, ia, lfil, droptol, alu, jlu, ju, iwk, w, jw, ierr, relax, izero, delta)
Definition: ilut.f90:51
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_residual()

subroutine imslinearbasemodule::ims_base_residual ( integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
real(dp), dimension(neq), intent(in)  X,
real(dp), dimension(neq), intent(in)  B,
real(dp), dimension(neq), intent(inout)  D,
real(dp), dimension(nja), intent(in)  A,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA 
)

Subroutine to calculate the residual.

Parameters
[in]neqlength of vectors
[in]njalength of coefficient matrix
[in]xdependent variable
[in]bright-hand side
[in,out]dresidual
[in]acoefficient matrix
[in]iaCRS row pointers
[in]jaCRS column pointers

Definition at line 1327 of file ImsLinearBase.f90.

1328  ! -- dummy variables
1329  integer(I4B), INTENT(IN) :: NEQ !< length of vectors
1330  integer(I4B), INTENT(IN) :: NJA !< length of coefficient matrix
1331  real(DP), DIMENSION(NEQ), INTENT(IN) :: X !< dependent variable
1332  real(DP), DIMENSION(NEQ), INTENT(IN) :: B !< right-hand side
1333  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: D !< residual
1334  real(DP), DIMENSION(NJA), INTENT(IN) :: A !< coefficient matrix
1335  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< CRS row pointers
1336  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< CRS column pointers
1337  ! -- local variables
1338  integer(I4B) :: n
1339  ! -- code
1340  !
1341  ! -- calculate matrix-vector product
1342  call amux(neq, x, d, a, ja, ia)
1343  !
1344  ! -- subtract matrix-vector product from right-hand side
1345  DO n = 1, neq
1346  d(n) = b(n) - d(n)
1347  END DO
1348  !
1349  ! -- return
1350  RETURN
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_scale()

subroutine imslinearbasemodule::ims_base_scale ( integer(i4b), intent(in)  IOPT,
integer(i4b), intent(in)  ISCL,
integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA,
real(dp), dimension(nja), intent(inout)  AMAT,
real(dp), dimension(neq), intent(inout)  X,
real(dp), dimension(neq), intent(inout)  B,
real(dp), dimension(neq), intent(inout)  DSCALE,
real(dp), dimension(neq), intent(inout)  DSCALE2 
)

Scale the coefficient matrix (AMAT), the right-hand side (B), and the estimate of the dependent variable (X).

Parameters
[in]ioptflag to scale (0) or unscale the system of equations
[in]isclscaling option (1) symmetric (2) L-2 norm
[in]neqnumber of equations
[in]njanumber of non-zero entries
[in]iaCRS row pointer
[in]jaCRS column pointer
[in,out]amatcoefficient matrix
[in,out]xdependent variable
[in,out]bright-hand side
[in,out]dscalefirst scaling vector
[in,out]dscale2second scaling vector

Definition at line 628 of file ImsLinearBase.f90.

630  ! -- dummy variables
631  integer(I4B), INTENT(IN) :: IOPT !< flag to scale (0) or unscale the system of equations
632  integer(I4B), INTENT(IN) :: ISCL !< scaling option (1) symmetric (2) L-2 norm
633  integer(I4B), INTENT(IN) :: NEQ !< number of equations
634  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
635  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< CRS row pointer
636  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< CRS column pointer
637  real(DP), DIMENSION(NJA), INTENT(INOUT) :: AMAT !< coefficient matrix
638  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: X !< dependent variable
639  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: B !< right-hand side
640  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: DSCALE !< first scaling vector
641  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: DSCALE2 !< second scaling vector
642  ! -- local variables
643  integer(I4B) :: i, n
644  integer(I4B) :: id, jc
645  integer(I4B) :: i0, i1
646  real(DP) :: v, c1, c2
647  !
648  ! -- SCALE SCALE AMAT, X, AND B
649  IF (iopt == 0) THEN
650  !
651  ! -- SYMMETRIC SCALING
652  SELECT CASE (iscl)
653  CASE (1)
654  DO n = 1, neq
655  id = ia(n)
656  v = amat(id)
657  c1 = done / sqrt(abs(v))
658  dscale(n) = c1
659  dscale2(n) = c1
660  END DO
661  !
662  ! -- SCALE AMAT -- AMAT = DSCALE(row) * AMAT(i) * DSCALE2(col)
663  DO n = 1, neq
664  c1 = dscale(n)
665  i0 = ia(n)
666  i1 = ia(n + 1) - 1
667  DO i = i0, i1
668  jc = ja(i)
669  c2 = dscale2(jc)
670  amat(i) = c1 * amat(i) * c2
671  END DO
672  END DO
673  !
674  ! -- L-2 NORM SCALING
675  CASE (2)
676  !
677  ! -- SCALE EACH ROW SO THAT THE L-2 NORM IS 1
678  DO n = 1, neq
679  c1 = dzero
680  i0 = ia(n)
681  i1 = ia(n + 1) - 1
682  DO i = i0, i1
683  c1 = c1 + amat(i) * amat(i)
684  END DO
685  c1 = sqrt(c1)
686  IF (c1 == dzero) THEN
687  c1 = done
688  ELSE
689  c1 = done / c1
690  END IF
691  dscale(n) = c1
692  !
693  ! -- INITIAL SCALING OF AMAT -- AMAT = DSCALE(row) * AMAT(i)
694  DO i = i0, i1
695  amat(i) = c1 * amat(i)
696  END DO
697  END DO
698  !
699  ! -- SCALE EACH COLUMN SO THAT THE L-2 NORM IS 1
700  DO n = 1, neq
701  dscale2(n) = dzero
702  END DO
703  c2 = dzero
704  DO n = 1, neq
705  i0 = ia(n)
706  i1 = ia(n + 1) - 1
707  DO i = i0, i1
708  jc = ja(i)
709  c2 = amat(i)
710  dscale2(jc) = dscale2(jc) + c2 * c2
711  END DO
712  END DO
713  DO n = 1, neq
714  c2 = dscale2(n)
715  IF (c2 == dzero) THEN
716  c2 = done
717  ELSE
718  c2 = done / sqrt(c2)
719  END IF
720  dscale2(n) = c2
721  END DO
722  !
723  ! -- FINAL SCALING OF AMAT -- AMAT = DSCALE2(col) * AMAT(i)
724  DO n = 1, neq
725  i0 = ia(n)
726  i1 = ia(n + 1) - 1
727  DO i = i0, i1
728  jc = ja(i)
729  c2 = dscale2(jc)
730  amat(i) = c2 * amat(i)
731  END DO
732  END DO
733  END SELECT
734  !
735  ! -- SCALE X AND B
736  DO n = 1, neq
737  c1 = dscale(n)
738  c2 = dscale2(n)
739  x(n) = x(n) / c2
740  b(n) = b(n) * c1
741  END DO
742  !
743  ! -- UNSCALE SCALE AMAT, X, AND B
744  ELSE
745  DO n = 1, neq
746  c1 = dscale(n)
747  i0 = ia(n)
748  i1 = ia(n + 1) - 1
749  !
750  ! -- UNSCALE AMAT
751  DO i = i0, i1
752  jc = ja(i)
753  c2 = dscale2(jc)
754  amat(i) = (done / c1) * amat(i) * (done / c2)
755  END DO
756  !
757  ! -- UNSCALE X AND B
758  c2 = dscale2(n)
759  x(n) = x(n) * c2
760  b(n) = b(n) / c1
761  END DO
762  END IF
763  !
764  ! -- RETURN
765  RETURN
Here is the caller graph for this function:

◆ ims_base_testcnvg()

subroutine imslinearbasemodule::ims_base_testcnvg ( integer(i4b), intent(in)  Icnvgopt,
integer(i4b), intent(inout)  Icnvg,
integer(i4b), intent(in)  Iiter,
real(dp), intent(in)  Dvmax,
real(dp), intent(in)  Rmax,
real(dp), intent(in)  Rmax0,
real(dp), intent(in)  Epfact,
real(dp), intent(in)  Dvclose,
real(dp), intent(in)  Rclose 
)

General routine for testing for solver convergence based on the user-specified convergence option (Icnvgopt).

Parameters
[in]icnvgoptconvergence option - see documentation for option
[in,out]icnvgflag indicating if convergence achieved (1) or not (0)
[in]iiterinner iteration number (used for strict convergence option)
[in]dvmaxmaximum dependent-variable change
[in]rmaxmaximum flow change
[in]rmax0initial flow change (initial L2-norm)
[in]epfactfactor for reducing convergence criteria in subsequent Picard iterations
[in]dvcloseMaximum depenendent-variable change allowed
[in]rcloseMaximum flow change allowed

Definition at line 1128 of file ImsLinearBase.f90.

1131  ! -- dummy variables
1132  integer(I4B), INTENT(IN) :: Icnvgopt !< convergence option - see documentation for option
1133  integer(I4B), INTENT(INOUT) :: Icnvg !< flag indicating if convergence achieved (1) or not (0)
1134  integer(I4B), INTENT(IN) :: Iiter !< inner iteration number (used for strict convergence option)
1135  real(DP), INTENT(IN) :: Dvmax !< maximum dependent-variable change
1136  real(DP), INTENT(IN) :: Rmax !< maximum flow change
1137  real(DP), INTENT(IN) :: Rmax0 !< initial flow change (initial L2-norm)
1138  real(DP), INTENT(IN) :: Epfact !< factor for reducing convergence criteria in subsequent Picard iterations
1139  real(DP), INTENT(IN) :: Dvclose !< Maximum depenendent-variable change allowed
1140  real(DP), INTENT(IN) :: Rclose !< Maximum flow change allowed
1141  ! -- code
1142  IF (icnvgopt == 0) THEN
1143  IF (abs(dvmax) <= dvclose .AND. abs(rmax) <= rclose) THEN
1144  icnvg = 1
1145  END IF
1146  ELSE IF (icnvgopt == 1) THEN
1147  IF (abs(dvmax) <= dvclose .AND. abs(rmax) <= rclose) THEN
1148  IF (iiter == 1) THEN
1149  icnvg = 1
1150  ELSE
1151  icnvg = -1
1152  END IF
1153  END IF
1154  ELSE IF (icnvgopt == 2) THEN
1155  IF (abs(dvmax) <= dvclose .OR. rmax <= rclose) THEN
1156  icnvg = 1
1157  ELSE IF (rmax <= rmax0 * epfact) THEN
1158  icnvg = -1
1159  END IF
1160  ELSE IF (icnvgopt == 3) THEN
1161  IF (abs(dvmax) <= dvclose) THEN
1162  icnvg = 1
1163  ELSE IF (rmax <= rmax0 * rclose) THEN
1164  icnvg = -1
1165  END IF
1166  ELSE IF (icnvgopt == 4) THEN
1167  IF (abs(dvmax) <= dvclose .AND. rmax <= rclose) THEN
1168  icnvg = 1
1169  ELSE IF (rmax <= rmax0 * epfact) THEN
1170  icnvg = -1
1171  END IF
1172  END IF
1173  !
1174  ! -- return
1175  RETURN
Here is the caller graph for this function:

◆ ims_calc_pcdims()

subroutine imslinearbasemodule::ims_calc_pcdims ( integer(i4b), intent(in)  neq,
integer(i4b), intent(in)  nja,
integer(i4b), dimension(:), intent(in)  ia,
integer(i4b), intent(in)  level,
integer(i4b), intent(in)  ipc,
integer(i4b), intent(inout)  niapc,
integer(i4b), intent(inout)  njapc,
integer(i4b), intent(inout)  njlu,
integer(i4b), intent(inout)  njw,
integer(i4b), intent(inout)  nwlu 
)
Parameters
[in]neqnr. of rows A
[in]njanr. of nonzeros A
[in]iaCSR row pointers A
[in]levelfill level ILU
[in]ipcIMS preconditioner type
[in,out]niapcwork array size
[in,out]njapcwork array size
[in,out]njluwork array size
[in,out]njwwork array size
[in,out]nwluwork array size

Definition at line 1178 of file ImsLinearBase.f90.

1180  integer(I4B), intent(in) :: neq !< nr. of rows A
1181  integer(I4B), intent(in) :: nja !< nr. of nonzeros A
1182  integer(I4B), dimension(:), intent(in) :: ia !< CSR row pointers A
1183  integer(I4B), intent(in) :: level !< fill level ILU
1184  integer(I4B), intent(in) :: ipc !< IMS preconditioner type
1185  integer(I4B), intent(inout) :: niapc !< work array size
1186  integer(I4B), intent(inout) :: njapc !< work array size
1187  integer(I4B), intent(inout) :: njlu !< work array size
1188  integer(I4B), intent(inout) :: njw !< work array size
1189  integer(I4B), intent(inout) :: nwlu !< work array size
1190  ! local
1191  integer(I4B) :: n, i
1192  integer(I4B) :: ijlu, ijw, iwlu, iwk
1193 
1194  ijlu = 1
1195  ijw = 1
1196  iwlu = 1
1197 
1198  ! ILU0 and MILU0
1199  niapc = neq
1200  njapc = nja
1201 
1202  ! ILUT and MILUT
1203  if (ipc == 3 .or. ipc == 4) then
1204  niapc = neq
1205  if (level > 0) then
1206  iwk = neq * (level * 2 + 1)
1207  else
1208  iwk = 0
1209  do n = 1, neq
1210  i = ia(n + 1) - ia(n)
1211  if (i > iwk) then
1212  iwk = i
1213  end if
1214  end do
1215  iwk = neq * iwk
1216  end if
1217  njapc = iwk
1218  ijlu = iwk
1219  ijw = 2 * neq
1220  iwlu = neq + 1
1221  end if
1222 
1223  njlu = ijlu
1224  njw = ijw
1225  nwlu = iwlu
1226 
Here is the caller graph for this function:

Variable Documentation

◆ parser

type(blockparsertype), private imslinearbasemodule::parser
private

Definition at line 19 of file ImsLinearBase.f90.

19  type(BlockParserType), private :: parser