MODFLOW 6  version 6.7.0.dev2
USGS Modular Hydrologic Model
ImsLinear.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b
5  izero, dzero, dprec, dsame, &
6  dem8, dem6, dem5, dem4, dem3, dem2, dem1, &
7  dhalf, done, dtwo, &
8  vdebug
18 
19  IMPLICIT NONE
20  private
21 
22  TYPE, PUBLIC :: imslineardatatype
23  character(len=LENMEMPATH) :: memorypath !< the path for storing variables in the memory manager
24  integer(I4B), POINTER :: iout => null() !< simulation listing file unit
25  integer(I4B), POINTER :: iprims => null() !< print flag
26  ! input variables (pointing to fields in input structure)
27  real(dp), pointer :: dvclose => null() !< dependent variable closure criterion
28  real(dp), pointer :: rclose => null() !< residual closure criterion
29  integer(I4B), pointer :: icnvgopt => null() !< convergence option
30  integer(I4B), pointer :: iter1 => null() !< max. iterations
31  integer(I4B), pointer :: ilinmeth => null() !< linear solver method
32  integer(I4B), pointer :: iscl => null() !< scaling method
33  integer(I4B), pointer :: iord => null() !< reordering method
34  integer(I4B), pointer :: north => null() !< number of orthogonalizations
35 
36  real(dp), pointer :: relax => null() !< relaxation factor
37  integer(I4B), pointer :: level => null() !< nr. of preconditioner levels
38  real(dp), pointer :: droptol => null() !< drop tolerance for preconditioner
39  !
40  integer(I4B), POINTER :: ipc => null() !< preconditioner flag
41  integer(I4B), POINTER :: iacpc => null() !< preconditioner CRS row pointers
42  integer(I4B), POINTER :: niterc => null() !<
43  integer(I4B), POINTER :: niabcgs => null() !< size of working vectors for BCGS linear accelerator
44  integer(I4B), POINTER :: niapc => null() !< preconditioner number of rows
45  integer(I4B), POINTER :: njapc => null() !< preconditioner number of non-zero entries
46  real(dp), POINTER :: epfact => null() !< factor for decreasing convergence criteria in seubsequent Picard iterations
47  real(dp), POINTER :: l2norm0 => null() !< initial L2 norm
48  ! -- ilut variables
49  integer(I4B), POINTER :: njlu => null() !< length of jlu work vector
50  integer(I4B), POINTER :: njw => null() !< length of jw work vector
51  integer(I4B), POINTER :: nwlu => null() !< length of wlu work vector
52  ! -- pointers to solution variables
53  integer(I4B), POINTER :: neq => null() !< number of equations (rows in matrix)
54  integer(I4B), POINTER :: nja => null() !< number of non-zero values in amat
55  integer(I4B), dimension(:), pointer, contiguous :: ia => null() !< position of start of each row
56  integer(I4B), dimension(:), pointer, contiguous :: ja => null() !< column pointer
57  real(dp), dimension(:), pointer, contiguous :: amat => null() !< coefficient matrix
58  real(dp), dimension(:), pointer, contiguous :: rhs => null() !< right-hand side of equation
59  real(dp), dimension(:), pointer, contiguous :: x => null() !< dependent variable
60  ! VECTORS
61  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: dscale => null() !< scaling factor
62  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: dscale2 => null() !< unscaling factor
63  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: iapc => null() !< position of start of each row in preconditioner matrix
64  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: japc => null() !< preconditioner matrix column pointer
65  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: apc => null() !< preconditioner coefficient matrix
66  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: lorder => null() !< reordering mapping
67  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: iorder => null() !< mapping to restore reordered matrix
68  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: iaro => null() !< position of start of each row in reordered matrix
69  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: jaro => null() !< reordered matrix column pointer
70  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: aro => null() !< reordered coefficient matrix
71  ! WORKING ARRAYS
72  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: iw => null() !< integer working array
73  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: w => null() !< real working array
74  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: id => null() !< integer working array
75  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: d => null() !< real working array
76  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: p => null() !< real working array
77  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: q => null() !< real working array
78  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: z => null() !< real working array
79  ! BICGSTAB WORKING ARRAYS
80  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: t => null() !< BICGSTAB real working array
81  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: v => null() !< BICGSTAB real working array
82  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: dhat => null() !< BICGSTAB real working array
83  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: phat => null() !< BICGSTAB real working array
84  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: qhat => null() !< rBICGSTAB eal working array
85  ! POINTERS FOR USE WITH BOTH ORIGINAL AND RCM ORDERINGS
86  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: ia0 => null() !< pointer to current CRS row pointers
87  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: ja0 => null() !< pointer to current CRS column pointers
88  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: a0 => null() !< pointer to current coefficient matrix
89  ! ILUT WORKING ARRAYS
90  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: jlu => null() !< ilut integer working array
91  integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: jw => null() !< ilut integer working array
92  real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: wlu => null() !< ilut real working array
93 
94  ! PROCEDURES (METHODS)
95  CONTAINS
96  PROCEDURE :: imslinear_allocate => imslinear_ar
97  procedure :: imslinear_summary
98  PROCEDURE :: imslinear_apply => imslinear_ap
99  procedure :: imslinear_da => imslinear_da
100  procedure, private :: allocate_scalars
101  ! -- PRIVATE PROCEDURES
102  PROCEDURE, PRIVATE :: set_imslinear_input => imslinear_set_input
103  END TYPE imslineardatatype
104 
105 CONTAINS
106 
107  !> @ brief Allocate storage and read data
108  !!
109  !! Allocate storage for linear accelerators and read data
110  !!
111  !<
112  SUBROUTINE imslinear_ar(this, NAME, IOUT, IPRIMS, MXITER, &
113  NEQ, matrix, RHS, X, linear_settings)
114  ! -- modules
117  use simmodule, only: store_error, count_errors, &
119  ! -- dummy variables
120  CLASS(imslineardatatype), INTENT(INOUT) :: this !< ImsLinearDataType instance
121  CHARACTER(LEN=LENSOLUTIONNAME), INTENT(IN) :: NAME !< solution name
122  integer(I4B), INTENT(IN) :: IOUT !< simulation listing file unit
123  integer(I4B), TARGET, INTENT(IN) :: IPRIMS !< print option
124  integer(I4B), INTENT(IN) :: MXITER !< maximum outer iterations
125  integer(I4B), TARGET, INTENT(IN) :: NEQ !< number of equations
126  class(matrixbasetype), pointer :: matrix
127  real(DP), DIMENSION(NEQ), TARGET, INTENT(INOUT) :: RHS !< right-hand side
128  real(DP), DIMENSION(NEQ), TARGET, INTENT(INOUT) :: X !< dependent variables
129  type(imslinearsettingstype), pointer :: linear_settings !< the settings form the IMS file
130  ! -- local variables
131  character(len=LINELENGTH) :: errmsg
132  integer(I4B) :: n
133  integer(I4B) :: i0
134  integer(I4B) :: iscllen, iolen
135 
136  !
137  ! -- DEFINE NAME
138  this%memoryPath = create_mem_path(name, 'IMSLINEAR')
139  !
140  ! -- SET pointers to IMS settings
141  this%DVCLOSE => linear_settings%dvclose
142  this%RCLOSE => linear_settings%rclose
143  this%ICNVGOPT => linear_settings%icnvgopt
144  this%ITER1 => linear_settings%iter1
145  this%ILINMETH => linear_settings%ilinmeth
146  this%ISCL => linear_settings%iscl
147  this%IORD => linear_settings%iord
148  this%NORTH => linear_settings%north
149  this%RELAX => linear_settings%relax
150  this%LEVEL => linear_settings%level
151  this%DROPTOL => linear_settings%droptol
152  !
153  ! -- SET POINTERS TO SOLUTION STORAGE
154  this%IPRIMS => iprims
155  this%NEQ => neq
156  call matrix%get_aij(this%IA, this%JA, this%AMAT)
157  call mem_allocate(this%NJA, 'NJA', this%memoryPath)
158  this%NJA = size(this%AMAT)
159  this%RHS => rhs
160  this%X => x
161  !
162  ! -- ALLOCATE SCALAR VARIABLES
163  call this%allocate_scalars()
164  !
165  ! -- initialize iout
166  this%iout = iout
167  !
168  ! -- DEFAULT VALUES
169  this%IPC = 0
170  !
171  this%IACPC = 0
172  !
173  ! -- PRINT A MESSAGE IDENTIFYING IMSLINEAR SOLVER PACKAGE
174  write (iout, 2000)
175 02000 FORMAT(1x, /1x, 'IMSLINEAR -- UNSTRUCTURED LINEAR SOLUTION', &
176  ' PACKAGE, VERSION 8, 04/28/2017')
177  !
178  ! -- DETERMINE PRECONDITIONER
179  IF (this%LEVEL > 0 .OR. this%DROPTOL > dzero) THEN
180  this%IPC = 3
181  ELSE
182  this%IPC = 1
183  END IF
184  IF (this%RELAX > dzero) THEN
185  this%IPC = this%IPC + 1
186  END IF
187  !
188  ! -- ERROR CHECKING FOR OPTIONS
189  IF (this%ISCL < 0) this%ISCL = 0
190  IF (this%ISCL > 2) THEN
191  WRITE (errmsg, '(A)') 'IMSLINEAR7AR ISCL MUST BE <= 2'
192  call store_error(errmsg)
193  END IF
194  IF (this%IORD < 0) this%IORD = 0
195  IF (this%IORD > 2) THEN
196  WRITE (errmsg, '(A)') 'IMSLINEAR7AR IORD MUST BE <= 2'
197  call store_error(errmsg)
198  END IF
199  IF (this%NORTH < 0) THEN
200  WRITE (errmsg, '(A)') 'IMSLINEAR7AR NORTH MUST >= 0'
201  call store_error(errmsg)
202  END IF
203  IF (this%RCLOSE == dzero) THEN
204  IF (this%ICNVGOPT /= 3) THEN
205  WRITE (errmsg, '(A)') 'IMSLINEAR7AR RCLOSE MUST > 0.0'
206  call store_error(errmsg)
207  END IF
208  END IF
209  IF (this%RELAX < dzero) THEN
210  WRITE (errmsg, '(A)') 'IMSLINEAR7AR RELAX MUST BE >= 0.0'
211  call store_error(errmsg)
212  END IF
213  IF (this%RELAX > done) THEN
214  WRITE (errmsg, '(A)') 'IMSLINEAR7AR RELAX MUST BE <= 1.0'
215  call store_error(errmsg)
216  END IF
217  !
218  ! -- INITIALIZE IMSLINEAR VARIABLES
219  this%NITERC = 0
220  !
221  ! -- ALLOCATE AND INITIALIZE MEMORY FOR IMSLINEAR
222  iscllen = 1
223  IF (this%ISCL .NE. 0) iscllen = neq
224  CALL mem_allocate(this%DSCALE, iscllen, 'DSCALE', trim(this%memoryPath))
225  CALL mem_allocate(this%DSCALE2, iscllen, 'DSCALE2', trim(this%memoryPath))
226  !
227  ! -- determine dimensions for preconditing arrays
228  call ims_calc_pcdims(this%NEQ, this%NJA, this%IA, this%LEVEL, this%IPC, &
229  this%NIAPC, this%NJAPC, this%NJLU, this%NJW, this%NWLU)
230  !
231  ! -- ALLOCATE BASE PRECONDITIONER VECTORS
232  CALL mem_allocate(this%IAPC, this%NIAPC + 1, 'IAPC', trim(this%memoryPath))
233  CALL mem_allocate(this%JAPC, this%NJAPC, 'JAPC', trim(this%memoryPath))
234  CALL mem_allocate(this%APC, this%NJAPC, 'APC', trim(this%memoryPath))
235  !
236  ! -- ALLOCATE MEMORY FOR ILU0 AND MILU0 NON-ZERO ROW ENTRY VECTOR
237  CALL mem_allocate(this%IW, this%NIAPC, 'IW', trim(this%memoryPath))
238  CALL mem_allocate(this%W, this%NIAPC, 'W', trim(this%memoryPath))
239  !
240  ! -- ALLOCATE MEMORY FOR ILUT VECTORS
241  CALL mem_allocate(this%JLU, this%NJLU, 'JLU', trim(this%memoryPath))
242  CALL mem_allocate(this%JW, this%NJW, 'JW', trim(this%memoryPath))
243  CALL mem_allocate(this%WLU, this%NWLU, 'WLU', trim(this%memoryPath))
244  !
245  ! -- GENERATE IAPC AND JAPC FOR ILU0 AND MILU0
246  IF (this%IPC == 1 .OR. this%IPC == 2) THEN
247  CALL ims_base_pccrs(this%NEQ, this%NJA, this%IA, this%JA, &
248  this%IAPC, this%JAPC)
249  END IF
250  !
251  ! -- ALLOCATE SPACE FOR PERMUTATION VECTOR
252  i0 = 1
253  iolen = 1
254  IF (this%IORD .NE. 0) THEN
255  i0 = this%NEQ
256  iolen = this%NJA
257  END IF
258  CALL mem_allocate(this%LORDER, i0, 'LORDER', trim(this%memoryPath))
259  CALL mem_allocate(this%IORDER, i0, 'IORDER', trim(this%memoryPath))
260  CALL mem_allocate(this%IARO, i0 + 1, 'IARO', trim(this%memoryPath))
261  CALL mem_allocate(this%JARO, iolen, 'JARO', trim(this%memoryPath))
262  CALL mem_allocate(this%ARO, iolen, 'ARO', trim(this%memoryPath))
263  !
264  ! -- ALLOCATE WORKING VECTORS FOR IMSLINEAR SOLVER
265  CALL mem_allocate(this%ID, this%NEQ, 'ID', trim(this%memoryPath))
266  CALL mem_allocate(this%D, this%NEQ, 'D', trim(this%memoryPath))
267  CALL mem_allocate(this%P, this%NEQ, 'P', trim(this%memoryPath))
268  CALL mem_allocate(this%Q, this%NEQ, 'Q', trim(this%memoryPath))
269  CALL mem_allocate(this%Z, this%NEQ, 'Z', trim(this%memoryPath))
270  !
271  ! -- ALLOCATE MEMORY FOR BCGS WORKING ARRAYS
272  this%NIABCGS = 1
273  IF (this%ILINMETH == 2) THEN
274  this%NIABCGS = this%NEQ
275  END IF
276  CALL mem_allocate(this%T, this%NIABCGS, 'T', trim(this%memoryPath))
277  CALL mem_allocate(this%V, this%NIABCGS, 'V', trim(this%memoryPath))
278  CALL mem_allocate(this%DHAT, this%NIABCGS, 'DHAT', trim(this%memoryPath))
279  CALL mem_allocate(this%PHAT, this%NIABCGS, 'PHAT', trim(this%memoryPath))
280  CALL mem_allocate(this%QHAT, this%NIABCGS, 'QHAT', trim(this%memoryPath))
281  !
282  ! -- INITIALIZE IMSLINEAR VECTORS
283  DO n = 1, iscllen
284  this%DSCALE(n) = done
285  this%DSCALE2(n) = done
286  END DO
287  DO n = 1, this%NJAPC
288  this%APC(n) = dzero
289  END DO
290  !
291  ! -- WORKING VECTORS
292  DO n = 1, this%NEQ
293  this%ID(n) = izero
294  this%D(n) = dzero
295  this%P(n) = dzero
296  this%Q(n) = dzero
297  this%Z(n) = dzero
298  END DO
299  DO n = 1, this%NIAPC
300  this%IW(n) = izero
301  this%W(n) = dzero
302  END DO
303  !
304  ! -- BCGS WORKING VECTORS
305  DO n = 1, this%NIABCGS
306  this%T(n) = dzero
307  this%V(n) = dzero
308  this%DHAT(n) = dzero
309  this%PHAT(n) = dzero
310  this%QHAT(n) = dzero
311  END DO
312  !
313  ! -- ILUT AND MILUT WORKING VECTORS
314  DO n = 1, this%NJLU
315  this%JLU(n) = izero
316  END DO
317  DO n = 1, this%NJW
318  this%JW(n) = izero
319  END DO
320  DO n = 1, this%NWLU
321  this%WLU(n) = dzero
322  END DO
323  !
324  ! -- REORDERING VECTORS
325  DO n = 1, i0 + 1
326  this%IARO(n) = izero
327  END DO
328  DO n = 1, iolen
329  this%JARO(n) = izero
330  this%ARO(n) = dzero
331  END DO
332  !
333  ! -- REVERSE CUTHILL MCKEE AND MINIMUM DEGREE ORDERING
334  IF (this%IORD .NE. 0) THEN
335  CALL ims_base_calc_order(this%IORD, this%NEQ, this%NJA, this%IA, &
336  this%JA, this%LORDER, this%IORDER)
337  END IF
338  !
339  ! -- ALLOCATE MEMORY FOR STORING ITERATION CONVERGENCE DATA
340  end SUBROUTINE imslinear_ar
341 
342  !> @ brief Write summary of settings
343  !!
344  !! Write summary of linear accelerator settings.
345  !!
346  !<
347  subroutine imslinear_summary(this, mxiter)
348  ! -- dummy variables
349  class(imslineardatatype), intent(inout) :: this !< ImsLinearDataType instance
350  integer(I4B), intent(in) :: mxiter !< maximum number of outer iterations
351  ! -- local variables
352  CHARACTER(LEN=10) :: clin(0:2)
353  CHARACTER(LEN=31) :: clintit(0:2)
354  CHARACTER(LEN=20) :: cipc(0:4)
355  CHARACTER(LEN=20) :: cscale(0:2)
356  CHARACTER(LEN=25) :: corder(0:2)
357  CHARACTER(LEN=16), DIMENSION(0:4) :: ccnvgopt
358  CHARACTER(LEN=15) :: clevel
359  CHARACTER(LEN=15) :: cdroptol
360  integer(I4B) :: i
361  integer(I4B) :: j
362  ! -- data
363  DATA clin/'UNKNOWN ', &
364  &'CG ', &
365  &'BCGS '/
366  DATA clintit/' UNKNOWN ', &
367  &' CONJUGATE-GRADIENT ', &
368  &'BICONJUGATE-GRADIENT STABILIZED'/
369  DATA cipc/'UNKNOWN ', &
370  &'INCOMPLETE LU ', &
371  &'MOD. INCOMPLETE LU ', &
372  &'INCOMPLETE LUT ', &
373  &'MOD. INCOMPLETE LUT '/
374  DATA cscale/'NO SCALING ', &
375  &'SYMMETRIC SCALING ', &
376  &'L2 NORM SCALING '/
377  DATA corder/'ORIGINAL ORDERING ', &
378  &'RCM ORDERING ', &
379  &'MINIMUM DEGREE ORDERING '/
380  DATA ccnvgopt/'INFINITY NORM ', &
381  &'INFINITY NORM S ', &
382  &'L2 NORM ', &
383  &'RELATIVE L2NORM ', &
384  &'L2 NORM W. REL. '/
385  ! -- formats
386 02010 FORMAT(1x, /, 7x, 'SOLUTION BY THE', 1x, a31, 1x, 'METHOD', &
387  /, 1x, 66('-'), /, &
388  ' MAXIMUM OF ', i0, ' CALLS OF SOLUTION ROUTINE', /, &
389  ' MAXIMUM OF ', i0, &
390  ' INTERNAL ITERATIONS PER CALL TO SOLUTION ROUTINE', /, &
391  ' LINEAR ACCELERATION METHOD =', 1x, a, /, &
392  ' MATRIX PRECONDITIONING TYPE =', 1x, a, /, &
393  ' MATRIX SCALING APPROACH =', 1x, a, /, &
394  ' MATRIX REORDERING APPROACH =', 1x, a, /, &
395  ' NUMBER OF ORTHOGONALIZATIONS =', 1x, i0, /, &
396  ' HEAD CHANGE CRITERION FOR CLOSURE =', e15.5, /, &
397  ' RESIDUAL CHANGE CRITERION FOR CLOSURE =', e15.5, /, &
398  ' RESIDUAL CONVERGENCE OPTION =', 1x, i0, /, &
399  ' RESIDUAL CONVERGENCE NORM =', 1x, a, /, &
400  ' RELAXATION FACTOR =', e15.5)
401 02015 FORMAT(' NUMBER OF LEVELS =', a15, /, &
402  ' DROP TOLERANCE =', a15, //)
403 2030 FORMAT(1x, a20, 1x, 6(i6, 1x))
404 2040 FORMAT(1x, 20('-'), 1x, 6(6('-'), 1x))
405 2050 FORMAT(1x, 62('-'),/) !
406 ! -- -----------------------------------------------------------
407  !
408  ! -- initialize clevel and cdroptol
409  clevel = ''
410  cdroptol = ''
411  !
412  ! -- write common variables to all linear accelerators
413  write (this%iout, 2010) &
414  clintit(this%ILINMETH), mxiter, this%ITER1, &
415  clin(this%ILINMETH), cipc(this%IPC), &
416  cscale(this%ISCL), corder(this%IORD), &
417  this%NORTH, this%DVCLOSE, this%RCLOSE, &
418  this%ICNVGOPT, ccnvgopt(this%ICNVGOPT), &
419  this%RELAX
420  if (this%level > 0) then
421  write (clevel, '(i15)') this%level
422  end if
423  if (this%droptol > dzero) then
424  write (cdroptol, '(e15.5)') this%droptol
425  end if
426  IF (this%level > 0 .or. this%droptol > dzero) THEN
427  write (this%iout, 2015) trim(adjustl(clevel)), &
428  trim(adjustl(cdroptol))
429  ELSE
430  write (this%iout, '(//)')
431  END IF
432 
433  if (this%iord /= 0) then
434  !
435  ! -- WRITE SUMMARY OF REORDERING INFORMATION TO LIST FILE
436  if (this%iprims == 2) then
437  DO i = 1, this%neq, 6
438  write (this%iout, 2030) 'ORIGINAL NODE :', &
439  (j, j=i, min(i + 5, this%neq))
440  write (this%iout, 2040)
441  write (this%iout, 2030) 'REORDERED INDEX :', &
442  (this%lorder(j), j=i, min(i + 5, this%neq))
443  write (this%iout, 2030) 'REORDERED NODE :', &
444  (this%iorder(j), j=i, min(i + 5, this%neq))
445  write (this%iout, 2050)
446  END DO
447  END IF
448  end if
449  end subroutine imslinear_summary
450 
451  !> @ brief Allocate and initialize scalars
452  !!
453  !! Allocate and initialize linear accelerator scalars
454  !!
455  !<
456  subroutine allocate_scalars(this)
457  ! -- modules
459  ! -- dummy variables
460  class(imslineardatatype), intent(inout) :: this !< ImsLinearDataType instance
461  !
462  ! -- allocate scalars
463  call mem_allocate(this%iout, 'IOUT', this%memoryPath)
464  call mem_allocate(this%ipc, 'IPC', this%memoryPath)
465  call mem_allocate(this%iacpc, 'IACPC', this%memoryPath)
466  call mem_allocate(this%niterc, 'NITERC', this%memoryPath)
467  call mem_allocate(this%niabcgs, 'NIABCGS', this%memoryPath)
468  call mem_allocate(this%niapc, 'NIAPC', this%memoryPath)
469  call mem_allocate(this%njapc, 'NJAPC', this%memoryPath)
470  call mem_allocate(this%epfact, 'EPFACT', this%memoryPath)
471  call mem_allocate(this%l2norm0, 'L2NORM0', this%memoryPath)
472  call mem_allocate(this%njlu, 'NJLU', this%memoryPath)
473  call mem_allocate(this%njw, 'NJW', this%memoryPath)
474  call mem_allocate(this%nwlu, 'NWLU', this%memoryPath)
475  !
476  ! -- initialize scalars
477  this%iout = 0
478  this%ipc = 0
479  this%iacpc = 0
480  this%niterc = 0
481  this%niabcgs = 0
482  this%niapc = 0
483  this%njapc = 0
484  this%epfact = dzero
485  this%l2norm0 = 0
486  this%njlu = 0
487  this%njw = 0
488  this%nwlu = 0
489  end subroutine allocate_scalars
490 
491  !> @ brief Deallocate memory
492  !!
493  !! Deallocate linear accelerator memory.
494  !!
495  !<
496  subroutine imslinear_da(this)
497  ! -- modules
499  ! -- dummy variables
500  class(imslineardatatype), intent(inout) :: this !< linear datatype instance
501  !
502  ! -- arrays
503  call mem_deallocate(this%dscale)
504  call mem_deallocate(this%dscale2)
505  call mem_deallocate(this%iapc)
506  call mem_deallocate(this%japc)
507  call mem_deallocate(this%apc)
508  call mem_deallocate(this%iw)
509  call mem_deallocate(this%w)
510  call mem_deallocate(this%jlu)
511  call mem_deallocate(this%jw)
512  call mem_deallocate(this%wlu)
513  call mem_deallocate(this%lorder)
514  call mem_deallocate(this%iorder)
515  call mem_deallocate(this%iaro)
516  call mem_deallocate(this%jaro)
517  call mem_deallocate(this%aro)
518  call mem_deallocate(this%id)
519  call mem_deallocate(this%d)
520  call mem_deallocate(this%p)
521  call mem_deallocate(this%q)
522  call mem_deallocate(this%z)
523  call mem_deallocate(this%t)
524  call mem_deallocate(this%v)
525  call mem_deallocate(this%dhat)
526  call mem_deallocate(this%phat)
527  call mem_deallocate(this%qhat)
528  !
529  ! -- scalars
530  call mem_deallocate(this%iout)
531  call mem_deallocate(this%ipc)
532  call mem_deallocate(this%iacpc)
533  call mem_deallocate(this%niterc)
534  call mem_deallocate(this%niabcgs)
535  call mem_deallocate(this%niapc)
536  call mem_deallocate(this%njapc)
537  call mem_deallocate(this%epfact)
538  call mem_deallocate(this%l2norm0)
539  call mem_deallocate(this%njlu)
540  call mem_deallocate(this%njw)
541  call mem_deallocate(this%nwlu)
542  call mem_deallocate(this%NJA)
543  !
544  ! -- nullify pointers
545  nullify (this%iprims)
546  nullify (this%neq)
547  nullify (this%nja)
548  nullify (this%ia)
549  nullify (this%ja)
550  nullify (this%amat)
551  nullify (this%rhs)
552  nullify (this%x)
553  end subroutine imslinear_da
554 
555  !> @ brief Set default settings
556  !!
557  !! Set default linear accelerator settings.
558  !!
559  !<
560  SUBROUTINE imslinear_set_input(this, IFDPARAM)
561  ! -- dummy variables
562  CLASS(imslineardatatype), INTENT(INOUT) :: this !< ImsLinearDataType instance
563  integer(I4B), INTENT(IN) :: IFDPARAM !< complexity option
564  ! -- code
565  SELECT CASE (ifdparam)
566  !
567  ! -- Simple option
568  CASE (1)
569  this%ITER1 = 50
570  this%ILINMETH = 1
571  this%IPC = 1
572  this%ISCL = 0
573  this%IORD = 0
574  this%DVCLOSE = dem3
575  this%RCLOSE = dem1
576  this%RELAX = dzero
577  this%LEVEL = 0
578  this%DROPTOL = dzero
579  this%NORTH = 0
580  !
581  ! -- Moderate
582  CASE (2)
583  this%ITER1 = 100
584  this%ILINMETH = 2
585  this%IPC = 2
586  this%ISCL = 0
587  this%IORD = 0
588  this%DVCLOSE = dem2
589  this%RCLOSE = dem1
590  this%RELAX = 0.97d0
591  this%LEVEL = 0
592  this%DROPTOL = dzero
593  this%NORTH = 0
594  !
595  ! -- Complex
596  CASE (3)
597  this%ITER1 = 500
598  this%ILINMETH = 2
599  this%IPC = 3
600  this%ISCL = 0
601  this%IORD = 0
602  this%DVCLOSE = dem1
603  this%RCLOSE = dem1
604  this%RELAX = dzero
605  this%LEVEL = 5
606  this%DROPTOL = dem4
607  this%NORTH = 2
608  END SELECT
609  end SUBROUTINE imslinear_set_input
610 
611  !> @ brief Base linear accelerator subroutine
612  !!
613  !! Base linear accelerator subroutine that scales and reorders
614  !! the system of equations, if necessary, updates the preconditioner,
615  !! and calls the appropriate linear accelerator.
616  !!
617  !<
618  SUBROUTINE imslinear_ap(this, ICNVG, KSTP, KITER, IN_ITER, &
619  NCONV, CONVNMOD, CONVMODSTART, &
620  CACCEL, summary)
621  ! -- modules
622  USE simmodule
623  ! -- dummy variables
624  CLASS(imslineardatatype), INTENT(INOUT) :: this !< ImsLinearDataType instance
625  integer(I4B), INTENT(INOUT) :: ICNVG !< convergence flag (1) non-convergence (0)
626  integer(I4B), INTENT(IN) :: KSTP !< time step number
627  integer(I4B), INTENT(IN) :: KITER !< outer iteration number
628  integer(I4B), INTENT(INOUT) :: IN_ITER !< inner iteration number
629  ! -- convergence information dummy variables
630  integer(I4B), INTENT(IN) :: NCONV !<
631  integer(I4B), INTENT(IN) :: CONVNMOD !<
632  integer(I4B), DIMENSION(CONVNMOD + 1), INTENT(INOUT) :: CONVMODSTART !<
633  character(len=31), DIMENSION(NCONV), INTENT(INOUT) :: CACCEL !<
634  type(convergencesummarytype), pointer, intent(in) :: summary !< Convergence summary report
635  ! -- local variables
636  integer(I4B) :: n
637  integer(I4B) :: innerit
638  integer(I4B) :: irc
639  integer(I4B) :: itmax
640  real(DP) :: dnrm2
641  !
642  ! -- set epfact based on timestep
643  this%EPFACT = ims_base_epfact(this%ICNVGOPT, kstp)
644  !
645  ! -- SCALE PROBLEM
646  IF (this%ISCL .NE. 0) THEN
647  CALL ims_base_scale(0, this%ISCL, &
648  this%NEQ, this%NJA, this%IA, this%JA, &
649  this%AMAT, this%X, this%RHS, &
650  this%DSCALE, this%DSCALE2)
651  END IF
652  !
653  ! -- PERMUTE ROWS, COLUMNS, AND RHS
654  IF (this%IORD /= 0) THEN
655  CALL dperm(this%NEQ, this%AMAT, this%JA, this%IA, &
656  this%ARO, this%JARO, this%IARO, &
657  this%LORDER, this%ID, 1)
658  CALL dvperm(this%NEQ, this%X, this%LORDER)
659  CALL dvperm(this%NEQ, this%RHS, this%LORDER)
660  this%IA0 => this%IARO
661  this%JA0 => this%JARO
662  this%A0 => this%ARO
663  ELSE
664  this%IA0 => this%IA
665  this%JA0 => this%JA
666  this%A0 => this%AMAT
667  END IF
668  !
669  ! -- UPDATE PRECONDITIONER
670  CALL ims_base_pcu(this%iout, this%NJA, this%NEQ, this%NIAPC, this%NJAPC, &
671  this%IPC, this%RELAX, this%A0, this%IA0, this%JA0, &
672  this%APC, this%IAPC, this%JAPC, this%IW, this%W, &
673  this%LEVEL, this%DROPTOL, this%NJLU, this%NJW, &
674  this%NWLU, this%JLU, this%JW, this%WLU)
675  !
676  ! -- INITIALIZE SOLUTION VARIABLE AND ARRAYS
677  IF (kiter == 1) then
678  this%NITERC = 0
679  summary%iter_cnt = 0
680  end if
681  irc = 1
682  icnvg = 0
683  DO n = 1, this%NEQ
684  this%D(n) = dzero
685  this%P(n) = dzero
686  this%Q(n) = dzero
687  this%Z(n) = dzero
688  END DO
689  !
690  ! -- CALCULATE INITIAL RESIDUAL
691  call ims_base_residual(this%NEQ, this%NJA, this%X, this%RHS, this%D, &
692  this%A0, this%IA0, this%JA0)
693  this%L2NORM0 = dnrm2(this%NEQ, this%D, 1)
694  !
695  ! -- CHECK FOR EXACT SOLUTION
696  itmax = this%ITER1
697  IF (this%L2NORM0 == dzero) THEN
698  itmax = 0
699  icnvg = 1
700  END IF
701  !
702  ! -- SOLUTION BY THE CONJUGATE GRADIENT METHOD
703  IF (this%ILINMETH == 1) THEN
704  CALL ims_base_cg(icnvg, itmax, innerit, &
705  this%NEQ, this%NJA, this%NIAPC, this%NJAPC, &
706  this%IPC, this%ICNVGOPT, this%NORTH, &
707  this%DVCLOSE, this%RCLOSE, this%L2NORM0, &
708  this%EPFACT, this%IA0, this%JA0, this%A0, &
709  this%IAPC, this%JAPC, this%APC, &
710  this%X, this%RHS, this%D, this%P, this%Q, this%Z, &
711  this%NJLU, this%IW, this%JLU, &
712  nconv, convnmod, convmodstart, &
713  caccel, summary)
714  !
715  ! -- SOLUTION BY THE BICONJUGATE GRADIENT STABILIZED METHOD
716  ELSE IF (this%ILINMETH == 2) THEN
717  CALL ims_base_bcgs(icnvg, itmax, innerit, &
718  this%NEQ, this%NJA, this%NIAPC, this%NJAPC, &
719  this%IPC, this%ICNVGOPT, this%NORTH, &
720  this%ISCL, this%DSCALE, &
721  this%DVCLOSE, this%RCLOSE, this%L2NORM0, &
722  this%EPFACT, this%IA0, this%JA0, this%A0, &
723  this%IAPC, this%JAPC, this%APC, &
724  this%X, this%RHS, this%D, this%P, this%Q, &
725  this%T, this%V, this%DHAT, this%PHAT, this%QHAT, &
726  this%NJLU, this%IW, this%JLU, &
727  nconv, convnmod, convmodstart, &
728  caccel, summary)
729  END IF
730  !
731  ! -- BACK PERMUTE AMAT, SOLUTION, AND RHS
732  IF (this%IORD /= 0) THEN
733  CALL dperm(this%NEQ, this%A0, this%JA0, this%IA0, &
734  this%AMAT, this%JA, this%IA, &
735  this%IORDER, this%ID, 1)
736  CALL dvperm(this%NEQ, this%X, this%IORDER)
737  CALL dvperm(this%NEQ, this%RHS, this%IORDER)
738  END IF
739  !
740  ! -- UNSCALE PROBLEM
741  IF (this%ISCL .NE. 0) THEN
742  CALL ims_base_scale(1, this%ISCL, &
743  this%NEQ, this%NJA, this%IA, this%JA, &
744  this%AMAT, this%X, this%RHS, &
745  this%DSCALE, this%DSCALE2)
746  END IF
747  !
748  ! -- SET IMS INNER ITERATION NUMBER (IN_ITER) TO NUMBER OF
749  ! IMSLINEAR INNER ITERATIONS (innerit)
750  in_iter = innerit
751  end SUBROUTINE imslinear_ap
752 
753 END MODULE imslinearmodule
This module contains block parser methods.
Definition: BlockParser.f90:7
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter dsame
real constant for values that are considered the same based on machine precision
Definition: Constants.f90:122
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
integer(i4b), parameter lensolutionname
maximum length of the solution name
Definition: Constants.f90:21
real(dp), parameter dem8
real constant 1e-8
Definition: Constants.f90:111
real(dp), parameter dem1
real constant 1e-1
Definition: Constants.f90:103
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
real(dp), parameter dem3
real constant 1e-3
Definition: Constants.f90:106
integer(i4b), parameter izero
integer constant zero
Definition: Constants.f90:51
real(dp), parameter dem4
real constant 1e-4
Definition: Constants.f90:107
real(dp), parameter dem6
real constant 1e-6
Definition: Constants.f90:109
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
real(dp), parameter dem5
real constant 1e-5
Definition: Constants.f90:108
real(dp), parameter dprec
real constant machine precision
Definition: Constants.f90:120
@ vdebug
write debug output
Definition: Constants.f90:190
real(dp), parameter dem2
real constant 1e-2
Definition: Constants.f90:105
real(dp), parameter dtwo
real constant 2
Definition: Constants.f90:79
integer(i4b), parameter lenmempath
maximum length of the memory path
Definition: Constants.f90:27
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
This module contains the IMS linear accelerator subroutines.
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
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
subroutine ims_base_pccrs(NEQ, NJA, IA, JA, IAPC, JAPC)
@ brief Generate CRS pointers for the preconditioner
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
subroutine ims_calc_pcdims(neq, nja, ia, level, ipc, niapc, njapc, njlu, njw, nwlu)
subroutine ims_base_scale(IOPT, ISCL, NEQ, NJA, IA, JA, AMAT, X, B, DSCALE, DSCALE2)
@ brief Scale the coefficient matrix
real(dp) function ims_base_epfact(icnvgopt, kstp)
Function returning EPFACT.
subroutine ims_base_calc_order(IORD, NEQ, NJA, IA, JA, LORDER, IORDER)
@ brief Calculate LORDER AND IORDER
subroutine ims_base_residual(NEQ, NJA, X, B, D, A, IA, JA)
Calculate residual.
subroutine allocate_scalars(this)
@ brief Allocate and initialize scalars
Definition: ImsLinear.f90:457
subroutine imslinear_summary(this, mxiter)
@ brief Write summary of settings
Definition: ImsLinear.f90:348
subroutine imslinear_set_input(this, IFDPARAM)
@ brief Set default settings
Definition: ImsLinear.f90:561
subroutine imslinear_da(this)
@ brief Deallocate memory
Definition: ImsLinear.f90:497
subroutine imslinear_ar(this, NAME, IOUT, IPRIMS, MXITER, NEQ, matrix, RHS, X, linear_settings)
@ brief Allocate storage and read data
Definition: ImsLinear.f90:114
subroutine imslinear_ap(this, ICNVG, KSTP, KITER, IN_ITER, NCONV, CONVNMOD, CONVMODSTART, CACCEL, summary)
@ brief Base linear accelerator subroutine
Definition: ImsLinear.f90:621
This module defines variable data types.
Definition: kind.f90:8
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
subroutine, public deprecation_warning(cblock, cvar, cver, endmsg, iunit)
Store deprecation warning message.
Definition: Sim.f90:256
subroutine dvperm(n, x, perm)
Definition: sparsekit.f90:62
subroutine dperm(nrow, a, ja, ia, ao, jao, iao, perm, qperm, job)
Definition: sparsekit.f90:354
This structure stores the generic convergence info for a solution.