MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
NumericalSolution.f90
Go to the documentation of this file.
1 ! This is the numerical solution module.
2 
4  use kindmodule, only: dp, i4b, lgp
5  use errorutilmodule, only: pstop
6  use timermodule, only: code_timer
9  dprec, dzero, dem20, dem15, dem6, &
12  tableft, tabright, &
13  mnormal, mvalidate, &
16  use tablemodule, only: tabletype, table_cr
17  Use messagemodule, only: write_message
18  use mathutilmodule, only: is_close
19  use versionmodule, only: idevelopmode
23  use listmodule, only: listtype
24  use listsmodule, only: basesolutionlist
32  use sparsemodule, only: sparsematrix
33  use simvariablesmodule, only: iout, isim_mode, errmsg, &
35  use simstagesmodule
37  use imslinearmodule
45 
46  implicit none
47  private
48 
49  public :: numericalsolutiontype
53 
54  integer(I4B), parameter :: ims_solver = 1
55  integer(I4B), parameter :: petsc_solver = 2
56 
58  character(len=LENMEMPATH) :: memory_path !< the path for storing solution variables in the memory manager
59  character(len=LINELENGTH) :: fname !< input file name
60  character(len=16) :: solver_mode !< the type of solve: sequential, parallel, mayve block, etc.
61  type(listtype), pointer :: modellist !< list of models in solution
62  type(listtype), pointer :: exchangelist !< list of exchanges in solution
63  integer(I4B), pointer :: id !< solution number
64  integer(I4B), pointer :: iu !< input file unit
65  real(dp), pointer :: ttform !< timer - total formulation time
66  real(dp), pointer :: ttsoln !< timer - total solution time
67  integer(I4B), pointer :: isymmetric => null() !< flag indicating if matrix symmetry is required
68  integer(I4B), pointer :: neq => null() !< number of equations
69  integer(I4B), pointer :: matrix_offset => null() !< offset of linear system when part of distributed solution
70  class(linearsolverbasetype), pointer :: linear_solver => null() !< the linear solver for this solution
71  class(matrixbasetype), pointer :: system_matrix => null() !< sparse A-matrix for the system of equations
72  class(vectorbasetype), pointer :: vec_rhs => null() !< the right-hand side vector
73  class(vectorbasetype), pointer :: vec_x => null() !< the dependent-variable vector
74  real(dp), dimension(:), pointer, contiguous :: rhs => null() !< right-hand side vector values
75  real(dp), dimension(:), pointer, contiguous :: x => null() !< dependent-variable vector values
76  integer(I4B), dimension(:), pointer, contiguous :: active => null() !< active cell array
77  real(dp), dimension(:), pointer, contiguous :: xtemp => null() !< temporary vector for previous dependent-variable iterate
78  type(blockparsertype) :: parser !< block parser object
79  !
80  ! -- sparse matrix data
81  real(dp), pointer :: theta => null() !< under-relaxation theta
82  real(dp), pointer :: akappa => null() !< under-relaxation kappa
83  real(dp), pointer :: gamma => null() !< under-relaxation gamma
84  real(dp), pointer :: amomentum => null() !< under-relaxation momentum term
85  real(dp), pointer :: breduc => null() !< backtracking reduction factor
86  real(dp), pointer :: btol => null() !< backtracking tolerance
87  real(dp), pointer :: res_lim => null() !< backtracking residual threshold
88  real(dp), pointer :: dvclose => null() !< dependent-variable closure criteria
89  real(dp), pointer :: bigchold => null() !< cooley under-relaxation weight
90  real(dp), pointer :: bigch => null() !< under-relaxation maximum dependent-variable change
91  real(dp), pointer :: relaxold => null() !< under-relaxation previous relaxation factor
92  real(dp), pointer :: res_prev => null() !< previous L-2 norm
93  real(dp), pointer :: res_new => null() !< current L-2 norm
94  integer(I4B), pointer :: icnvg => null() !< convergence flag (1) non-convergence (0)
95  integer(I4B), pointer :: itertot_timestep => null() !< total nr. of linear solves per call to sln_ca
96  integer(I4B), pointer :: iouttot_timestep => null() !< total nr. of outer iterations per call to sln_ca
97  integer(I4B), pointer :: itertot_sim => null() !< total nr. of inner iterations for simulation
98  integer(I4B), pointer :: mxiter => null() !< maximum number of Picard iterations
99  integer(I4B), pointer :: linsolver => null() !< linear solver used (IMS, PETSC, ...)
100  integer(I4B), pointer :: nonmeth => null() !< under-relaxation method used
101  integer(I4B), pointer :: numtrack => null() !< maximum number of backtracks
102  integer(I4B), pointer :: iprims => null() !< solver print option
103  integer(I4B), pointer :: ibflag => null() !< backtracking flag (1) on (0) off
104  integer(I4B), dimension(:, :), pointer, contiguous :: lrch => null() !< location of the largest dependent-variable change at the end of a Picard iteration
105  real(dp), dimension(:), pointer, contiguous :: hncg => null() !< largest dependent-variable change at the end of a Picard iteration
106  real(dp), dimension(:), pointer, contiguous :: dxold => null() !< DBD under-relaxation previous dependent-variable change
107  real(dp), dimension(:), pointer, contiguous :: deold => null() !< DBD under-relaxation dependent-variable change variable
108  real(dp), dimension(:), pointer, contiguous :: wsave => null() !< DBD under-relaxation sign-change factor
109  real(dp), dimension(:), pointer, contiguous :: hchold => null() !< DBD under-relaxation weighted dependent-variable change
110  !
111  ! -- convergence summary information
112  character(len=31), dimension(:), pointer, contiguous :: caccel => null() !< convergence string
113  integer(I4B), pointer :: icsvouterout => null() !< Picard iteration CSV output flag and file unit
114  integer(I4B), pointer :: icsvinnerout => null() !< Inner iteration CSV output flag and file unit
115  integer(I4B), pointer :: nitermax => null() !< maximum number of iterations in a time step (maxiter * maxinner)
116  integer(I4B), pointer :: convnmod => null() !< number of models in the solution
117  integer(I4B), dimension(:), pointer, contiguous :: convmodstart => null() !< pointer to the start of each model in the convmod* arrays
118  !
119  ! -- refactoring
120  type(convergencesummarytype), pointer :: cnvg_summary => null() !< details on the convergence behavior within a timestep
121  type(imslinearsettingstype), pointer :: linear_settings => null() !< IMS settings for linear solver
122  !
123  ! -- pseudo-transient continuation
124  integer(I4B), pointer :: iallowptc => null() !< flag indicating if ptc applied this time step
125  integer(I4B), pointer :: iptcopt => null() !< option for how to calculate the initial PTC value (ptcdel0)
126  integer(I4B), pointer :: iptcout => null() !< PTC output flag and file unit
127  real(dp), pointer :: l2norm0 => null() !< L-2 norm at the start of the first Picard iteration
128  real(dp), pointer :: ptcdel => null() !< PTC delta value
129  real(dp), pointer :: ptcdel0 => null() !< initial PTC delta value
130  real(dp), pointer :: ptcexp => null() !< PTC exponent
131  !
132  ! -- timer handles
133  integer(I4B) :: tmr_prep_solve !< timer - prepare solve
134  integer(I4B) :: tmr_solve !< timer - solve
135  integer(I4B) :: tmr_final_solve !< timer - finalize solve
136  integer(I4B) :: tmr_formulate !< timer - formulate
137  integer(I4B) :: tmr_linsolve !< timer - linear solve
138  integer(I4B) :: tmr_flows !< timer - calculate flows
139  integer(I4B) :: tmr_budgets !< timer - calculate budgets
140  character(len=24) :: id_postfix !< solution id based postfix for timer titles
141  !
142  ! -- adaptive time step
143  real(dp), pointer :: atsfrac => null() !< adaptive time step faction
144  !
145  ! -- linear accelerator storage
146  type(imslineardatatype), pointer :: imslinear => null() !< IMS linear acceleration object
147  !
148  ! -- sparse object
149  type(sparsematrix) :: sparse !< sparse object
150  !
151  ! -- table objects
152  type(tabletype), pointer :: innertab => null() !< inner iteration table object
153  type(tabletype), pointer :: outertab => null() !< Picard iteration table object
154  !
155  ! -- for synchronization of exchanges
156  class(*), pointer :: synchronize_ctx => null()
157  procedure(synchronize_iface), pointer :: synchronize => null()
158 
159  contains
160  procedure :: sln_df
161  procedure :: sln_ar
162  procedure :: sln_dt
163  procedure :: sln_ad
164  procedure :: sln_ot
165  procedure :: sln_ca
166  procedure :: sln_fp
167  procedure :: sln_da
168  procedure :: add_model
169  procedure :: add_exchange
170  procedure :: get_models
171  procedure :: get_exchanges
172  procedure :: save
173 
174  ! 'protected' (this can be overridden)
175  procedure :: sln_has_converged
179  procedure :: sln_calc_ptc
180  procedure :: sln_underrelax
183  procedure :: apply_backtracking
184 
185  ! private
186  procedure, private :: sln_connect
187  procedure, private :: sln_reset
188  procedure, private :: sln_ls
189  procedure, private :: sln_setouter
190  procedure, private :: sln_backtracking
191  procedure, private :: sln_maxval
192  procedure, private :: sln_calcdx
193  procedure, private :: sln_calc_residual
194  procedure, private :: sln_l2norm
195  procedure, private :: sln_get_dxmax
196  procedure, private :: sln_get_loc
197  procedure, private :: sln_get_nodeu
198  procedure, private :: allocate_scalars
199  procedure, private :: allocate_arrays
200  procedure, private :: convergence_summary
201  procedure, private :: csv_convergence_summary
202  procedure, private :: sln_buildsystem
203  procedure, private :: writecsvheader
204  procedure, private :: writeptcinfotofile
205 
206  ! Expose these for use through the BMI/XMI:
207  procedure, public :: preparesolve
208  procedure, public :: solve
209  procedure, public :: finalizesolve
210 
211  end type numericalsolutiontype
212 
213  abstract interface
214  subroutine synchronize_iface(solution, stage, ctx)
215  import numericalsolutiontype
216  import i4b
217  class(numericalsolutiontype) :: solution
218  integer(I4B) :: stage
219  class(*), pointer :: ctx
220  end subroutine synchronize_iface
221  end interface
222 
223 contains
224 
225  !> @ brief Create a new solution
226  !!
227  !! Create a new solution using the data in filename, assign this new
228  !! solution an id number and store the solution in the basesolutionlist.
229  !! Also open the filename for later reading.
230  !!
231  !<
232  subroutine create_numerical_solution(num_sol, filename, id)
233  ! -- modules
234  use simvariablesmodule, only: iout
236  ! -- dummy variables
237  class(numericalsolutiontype), pointer :: num_sol !< the create solution
238  character(len=*), intent(in) :: filename !< solution input file name
239  integer(I4B), intent(in) :: id !< solution id
240  ! -- local variables
241  integer(I4B) :: inunit
242  class(basesolutiontype), pointer :: solbase => null()
243  character(len=LENSOLUTIONNAME) :: solutionname
244  !
245  ! -- Create a new solution and add it to the basesolutionlist container
246  solbase => num_sol
247  write (solutionname, '(a, i0)') 'SLN_', id
248  !
249  num_sol%name = solutionname
250  num_sol%memory_path = create_mem_path(solutionname)
251  allocate (num_sol%modellist)
252  allocate (num_sol%exchangelist)
253  !
254  call num_sol%allocate_scalars()
255  !
257  !
258  num_sol%id = id
259  !
260  ! -- Open solution input file for reading later after problem size is known
261  ! Check to see if the file is already opened, which can happen when
262  ! running in single model mode
263  inquire (file=filename, number=inunit)
264 
265  if (inunit < 0) inunit = getunit()
266  num_sol%iu = inunit
267  write (iout, '(/a,a)') ' Creating solution: ', num_sol%name
268  call openfile(num_sol%iu, iout, filename, 'IMS')
269  !
270  ! -- Initialize block parser
271  call num_sol%parser%Initialize(num_sol%iu, iout)
272  end subroutine create_numerical_solution
273 
274  !> @ brief Allocate scalars
275  !!
276  !! Allocate scalars for a new solution.
277  !!
278  !<
279  subroutine allocate_scalars(this)
280  ! -- modules
282  ! -- dummy variables
283  class(numericalsolutiontype) :: this
284  !
285  ! -- allocate scalars
286  call mem_allocate(this%id, 'ID', this%memory_path)
287  call mem_allocate(this%iu, 'IU', this%memory_path)
288  call mem_allocate(this%ttform, 'TTFORM', this%memory_path)
289  call mem_allocate(this%ttsoln, 'TTSOLN', this%memory_path)
290  call mem_allocate(this%isymmetric, 'ISYMMETRIC', this%memory_path)
291  call mem_allocate(this%neq, 'NEQ', this%memory_path)
292  call mem_allocate(this%matrix_offset, 'MATRIX_OFFSET', this%memory_path)
293  call mem_allocate(this%dvclose, 'DVCLOSE', this%memory_path)
294  call mem_allocate(this%bigchold, 'BIGCHOLD', this%memory_path)
295  call mem_allocate(this%bigch, 'BIGCH', this%memory_path)
296  call mem_allocate(this%relaxold, 'RELAXOLD', this%memory_path)
297  call mem_allocate(this%res_prev, 'RES_PREV', this%memory_path)
298  call mem_allocate(this%res_new, 'RES_NEW', this%memory_path)
299  call mem_allocate(this%icnvg, 'ICNVG', this%memory_path)
300  call mem_allocate(this%itertot_timestep, 'ITERTOT_TIMESTEP', this%memory_path)
301  call mem_allocate(this%iouttot_timestep, 'IOUTTOT_TIMESTEP', this%memory_path)
302  call mem_allocate(this%itertot_sim, 'INNERTOT_SIM', this%memory_path)
303  call mem_allocate(this%mxiter, 'MXITER', this%memory_path)
304  call mem_allocate(this%linsolver, 'LINSOLVER', this%memory_path)
305  call mem_allocate(this%nonmeth, 'NONMETH', this%memory_path)
306  call mem_allocate(this%iprims, 'IPRIMS', this%memory_path)
307  call mem_allocate(this%theta, 'THETA', this%memory_path)
308  call mem_allocate(this%akappa, 'AKAPPA', this%memory_path)
309  call mem_allocate(this%gamma, 'GAMMA', this%memory_path)
310  call mem_allocate(this%amomentum, 'AMOMENTUM', this%memory_path)
311  call mem_allocate(this%breduc, 'BREDUC', this%memory_path)
312  call mem_allocate(this%btol, 'BTOL', this%memory_path)
313  call mem_allocate(this%res_lim, 'RES_LIM', this%memory_path)
314  call mem_allocate(this%numtrack, 'NUMTRACK', this%memory_path)
315  call mem_allocate(this%ibflag, 'IBFLAG', this%memory_path)
316  call mem_allocate(this%icsvouterout, 'ICSVOUTEROUT', this%memory_path)
317  call mem_allocate(this%icsvinnerout, 'ICSVINNEROUT', this%memory_path)
318  call mem_allocate(this%nitermax, 'NITERMAX', this%memory_path)
319  call mem_allocate(this%convnmod, 'CONVNMOD', this%memory_path)
320  call mem_allocate(this%iallowptc, 'IALLOWPTC', this%memory_path)
321  call mem_allocate(this%iptcopt, 'IPTCOPT', this%memory_path)
322  call mem_allocate(this%iptcout, 'IPTCOUT', this%memory_path)
323  call mem_allocate(this%l2norm0, 'L2NORM0', this%memory_path)
324  call mem_allocate(this%ptcdel, 'PTCDEL', this%memory_path)
325  call mem_allocate(this%ptcdel0, 'PTCDEL0', this%memory_path)
326  call mem_allocate(this%ptcexp, 'PTCEXP', this%memory_path)
327  call mem_allocate(this%atsfrac, 'ATSFRAC', this%memory_path)
328  !
329  ! -- initialize scalars
330  this%isymmetric = 0
331  this%id = 0
332  this%iu = 0
333  this%ttform = dzero
334  this%ttsoln = dzero
335  this%neq = 0
336  this%dvclose = dzero
337  this%bigchold = dzero
338  this%bigch = dzero
339  this%relaxold = dzero
340  this%res_prev = dzero
341  this%icnvg = 0
342  this%itertot_timestep = 0
343  this%iouttot_timestep = 0
344  this%itertot_sim = 0
345  this%mxiter = 0
346  this%linsolver = ims_solver
347  this%nonmeth = 0
348  this%iprims = 0
349  this%theta = done
350  this%akappa = dzero
351  this%gamma = done
352  this%amomentum = dzero
353  this%breduc = dzero
354  this%btol = 0
355  this%res_lim = dzero
356  this%numtrack = 0
357  this%ibflag = 0
358  this%icsvouterout = 0
359  this%icsvinnerout = 0
360  this%nitermax = 0
361  this%convnmod = 0
362  this%iallowptc = 1
363  this%iptcopt = 0
364  this%iptcout = 0
365  this%l2norm0 = dzero
366  this%ptcdel = dzero
367  this%ptcdel0 = dzero
368  this%ptcexp = done
369  this%atsfrac = donethird
370  end subroutine allocate_scalars
371 
372  !> @ brief Allocate arrays
373  !!
374  !! Allocate arrays for a new solution.
375  !!
376  !<
377  subroutine allocate_arrays(this)
378  ! -- modules
380  ! -- dummy variables
381  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
382  ! -- local variables
383  class(numericalmodeltype), pointer :: mp => null()
384  integer(I4B) :: i
385  integer(I4B) :: ieq
386  !
387  ! -- initialize the number of models in the solution
388  this%convnmod = this%modellist%Count()
389  !
390  ! -- allocate arrays
391  call mem_allocate(this%active, this%neq, 'IACTIVE', this%memory_path)
392  call mem_allocate(this%xtemp, this%neq, 'XTEMP', this%memory_path)
393  call mem_allocate(this%dxold, this%neq, 'DXOLD', this%memory_path)
394  call mem_allocate(this%hncg, 0, 'HNCG', this%memory_path)
395  call mem_allocate(this%lrch, 3, 0, 'LRCH', this%memory_path)
396  call mem_allocate(this%wsave, 0, 'WSAVE', this%memory_path)
397  call mem_allocate(this%hchold, 0, 'HCHOLD', this%memory_path)
398  call mem_allocate(this%deold, 0, 'DEOLD', this%memory_path)
399  call mem_allocate(this%convmodstart, this%convnmod + 1, 'CONVMODSTART', &
400  this%memory_path)
401  !
402  ! -- initialize allocated arrays
403  do i = 1, this%neq
404  this%xtemp(i) = dzero
405  this%dxold(i) = dzero
406  this%active(i) = 1 !default is active
407  end do
408  !
409  ! -- initialize convmodstart
410  ieq = 1
411  this%convmodstart(1) = ieq
412  do i = 1, this%modellist%Count()
413  mp => getnumericalmodelfromlist(this%modellist, i)
414  ieq = ieq + mp%neq
415  this%convmodstart(i + 1) = ieq
416  end do
417  end subroutine allocate_arrays
418 
419  !> @ brief Define the solution
420  !!
421  !! Define a new solution. Must be called after the models and exchanges have
422  !! been added to solution. The order of the steps is (1) Allocate neq and nja,
423  !! (2) Assign model offsets and solution ids, (3) Allocate and initialize
424  !! the solution arrays, (4) Point each model's x and rhs arrays, and
425  !! (5) Initialize the sparsematrix instance
426  !!
427  !<
428  subroutine sln_df(this)
429  ! modules
432  ! -- dummy variables
433  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
434  ! -- local variables
435  class(numericalmodeltype), pointer :: mp => null()
436  integer(I4B) :: i
437  integer(I4B), allocatable, dimension(:) :: rowmaxnnz
438  integer(I4B) :: ncol, irow_start, irow_end
439  integer(I4B) :: mod_offset
440  !
441  ! -- set sol id and determine nr. of equation in this solution
442  do i = 1, this%modellist%Count()
443  mp => getnumericalmodelfromlist(this%modellist, i)
444  call mp%set_idsoln(this%id)
445  this%neq = this%neq + mp%neq
446  end do
447  !
448  ! -- set up the (possibly parallel) linear system
449  if (simulation_mode == 'PARALLEL') then
450  this%solver_mode = 'PETSC'
451  else
452  this%solver_mode = 'IMS'
453  end if
454  !
455  ! -- allocate settings structure
456  allocate (this%linear_settings)
457  !
458  ! -- create linear system matrix and compatible vectors
459  this%linear_solver => create_linear_solver(this%solver_mode, this%name)
460  this%system_matrix => this%linear_solver%create_matrix()
461  this%vec_x => this%system_matrix%create_vec_mm(this%neq, 'X', &
462  this%memory_path)
463  this%x => this%vec_x%get_array()
464  this%vec_rhs => this%system_matrix%create_vec_mm(this%neq, 'RHS', &
465  this%memory_path)
466  this%rhs => this%vec_rhs%get_array()
467  !
468  call this%vec_rhs%get_ownership_range(irow_start, irow_end)
469  ncol = this%vec_rhs%get_size()
470  !
471  ! -- calculate and set offsets
472  mod_offset = irow_start - 1
473  this%matrix_offset = irow_start - 1
474  do i = 1, this%modellist%Count()
475  mp => getnumericalmodelfromlist(this%modellist, i)
476  call mp%set_moffset(mod_offset)
477  mod_offset = mod_offset + mp%neq
478  end do
479  !
480  ! -- Allocate and initialize solution arrays
481  call this%allocate_arrays()
482  !
483  ! -- Create convergence summary report
484  allocate (this%cnvg_summary)
485  call this%cnvg_summary%init(this%modellist%Count(), this%convmodstart, &
486  this%memory_path)
487  !
488  ! -- Go through each model and point x, ibound, and rhs to solution
489  do i = 1, this%modellist%Count()
490  mp => getnumericalmodelfromlist(this%modellist, i)
491  call mp%set_xptr(this%x, this%matrix_offset, 'X', this%name)
492  call mp%set_rhsptr(this%rhs, this%matrix_offset, 'RHS', this%name)
493  call mp%set_iboundptr(this%active, this%matrix_offset, 'IBOUND', this%name)
494  end do
495  !
496  ! -- Create the sparsematrix instance
497  allocate (rowmaxnnz(this%neq))
498  do i = 1, this%neq
499  rowmaxnnz(i) = 4
500  end do
501  call this%sparse%init(this%neq, ncol, rowmaxnnz)
502  this%sparse%offset = this%matrix_offset
503  deallocate (rowmaxnnz)
504  !
505  ! -- Assign connections, fill ia/ja, map connections
506  call this%sln_connect()
507 
508  ! add timers
509  write (this%id_postfix, '(a,i0,a)') " (", this%id, ")"
510  this%tmr_prep_solve = -1
511  this%tmr_solve = -1
512  this%tmr_final_solve = -1
513  this%tmr_formulate = -1
514  this%tmr_linsolve = -1
515  this%tmr_flows = -1
516  this%tmr_budgets = -1
517 
518  end subroutine sln_df
519 
520  !> @ brief Allocate and read data
521  !!
522  !! Allocate and read data for a solution.
523  !!
524  !<
525  subroutine sln_ar(this)
526  ! -- modules
528  use simvariablesmodule, only: iout
531  ! -- dummy variables
532  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
533  ! -- local variables
534  class(numericalmodeltype), pointer :: mp => null()
535  class(numericalexchangetype), pointer :: cp => null()
536  character(len=linelength) :: warnmsg
537  character(len=linelength) :: keyword
538  character(len=linelength) :: fname
539  character(len=linelength) :: msg
540  integer(I4B) :: i
541  integer(I4B) :: ifdparam, mxvl, npp
542  integer(I4B) :: ierr
543  logical(LGP) :: isfound, endOfBlock
544  integer(I4B) :: ival
545  real(DP) :: rval
546  character(len=*), parameter :: fmtcsvout = &
547  "(4x, 'CSV OUTPUT WILL BE SAVED TO FILE: ', a, &
548  &/4x, 'OPENED ON UNIT: ', I7)"
549  character(len=*), parameter :: fmtptcout = &
550  "(4x, 'PTC OUTPUT WILL BE SAVED TO FILE: ', a, &
551  &/4x, 'OPENED ON UNIT: ', I7)"
552  character(len=*), parameter :: fmterrasym = &
553  "(a,' **',a,'** PRODUCES AN ASYMMETRIC COEFFICIENT MATRIX, BUT THE &
554  &CONJUGATE GRADIENT METHOD WAS SELECTED. USE BICGSTAB INSTEAD. ')"
555  !
556  ! identify package and initialize.
557  WRITE (iout, 1) this%iu
558 00001 FORMAT(1x, /1x, 'IMS -- ITERATIVE MODEL SOLUTION PACKAGE, VERSION 6', &
559  ', 4/28/2017', /, 9x, 'INPUT READ FROM UNIT', i5)
560  !
561  ! -- initialize
562  i = 1
563  ifdparam = 1
564  npp = 0
565  mxvl = 0
566  !
567  ! -- get options block
568  call this%parser%GetBlock('OPTIONS', isfound, ierr, &
569  supportopenclose=.true., blockrequired=.false.)
570  !
571  ! -- parse options block if detected
572  if (isfound) then
573  write (iout, '(/1x,a)') 'PROCESSING IMS OPTIONS'
574  do
575  call this%parser%GetNextLine(endofblock)
576  if (endofblock) exit
577  call this%parser%GetStringCaps(keyword)
578  select case (keyword)
579  case ('PRINT_OPTION')
580  call this%parser%GetStringCaps(keyword)
581  if (keyword .eq. 'NONE') then
582  this%iprims = 0
583  else if (keyword .eq. 'SUMMARY') then
584  this%iprims = 1
585  else if (keyword .eq. 'ALL') then
586  this%iprims = 2
587  else
588  write (errmsg, '(3a)') &
589  'Unknown IMS print option (', trim(keyword), ').'
590  call store_error(errmsg)
591  end if
592  case ('COMPLEXITY')
593  call this%parser%GetStringCaps(keyword)
594  if (keyword .eq. 'SIMPLE') then
595  ifdparam = 1
596  WRITE (iout, 21)
597  else if (keyword .eq. 'MODERATE') then
598  ifdparam = 2
599  WRITE (iout, 23)
600  else if (keyword .eq. 'COMPLEX') then
601  ifdparam = 3
602  WRITE (iout, 25)
603  else
604  write (errmsg, '(3a)') &
605  'Unknown IMS COMPLEXITY option (', trim(keyword), ').'
606  call store_error(errmsg)
607  end if
608  case ('CSV_OUTER_OUTPUT')
609  call this%parser%GetStringCaps(keyword)
610  if (keyword == 'FILEOUT') then
611  call this%parser%GetString(fname)
612  if (nr_procs > 1) then
613  call append_processor_id(fname, proc_id)
614  end if
615  this%icsvouterout = getunit()
616  call openfile(this%icsvouterout, iout, fname, 'CSV_OUTER_OUTPUT', &
617  filstat_opt='REPLACE')
618  write (iout, fmtcsvout) trim(fname), this%icsvouterout
619  else
620  write (errmsg, '(a)') 'Optional CSV_OUTER_OUTPUT '// &
621  'keyword must be followed by FILEOUT'
622  call store_error(errmsg)
623  end if
624  case ('CSV_INNER_OUTPUT')
625  call this%parser%GetStringCaps(keyword)
626  if (keyword == 'FILEOUT') then
627  call this%parser%GetString(fname)
628  if (nr_procs > 1) then
629  call append_processor_id(fname, proc_id)
630  end if
631  this%icsvinnerout = getunit()
632  call openfile(this%icsvinnerout, iout, fname, 'CSV_INNER_OUTPUT', &
633  filstat_opt='REPLACE')
634  write (iout, fmtcsvout) trim(fname), this%icsvinnerout
635  else
636  write (errmsg, '(a)') 'Optional CSV_INNER_OUTPUT '// &
637  'keyword must be followed by FILEOUT'
638  call store_error(errmsg)
639  end if
640  case ('NO_PTC')
641  call this%parser%GetStringCaps(keyword)
642  select case (keyword)
643  case ('ALL')
644  ival = 0
645  msg = 'ALL'
646  case ('FIRST')
647  ival = -1
648  msg = 'THE FIRST'
649  case default
650  ival = 0
651  msg = 'ALL'
652  end select
653  this%iallowptc = ival
654  write (iout, '(1x,A)') 'PSEUDO-TRANSIENT CONTINUATION DISABLED FOR'// &
655  ' '//trim(adjustl(msg))//' STRESS-PERIOD(S)'
656  case ('ATS_OUTER_MAXIMUM_FRACTION')
657  rval = this%parser%GetDouble()
658  if (rval < dzero .or. rval > dhalf) then
659  write (errmsg, '(a,G0)') 'Value for ATS_OUTER_MAXIMUM_FRAC must be &
660  &between 0 and 0.5. Found ', rval
661  call store_error(errmsg)
662  end if
663  this%atsfrac = rval
664  write (iout, '(1x,A,G0)') 'ADAPTIVE TIME STEP SETTING FOUND. FRACTION &
665  &OF OUTER MAXIMUM USED TO INCREASE OR DECREASE TIME STEP SIZE IS ',&
666  &this%atsfrac
667  !
668  ! -- DEPRECATED OPTIONS
669  case ('CSV_OUTPUT')
670  call this%parser%GetStringCaps(keyword)
671  if (keyword == 'FILEOUT') then
672  call this%parser%GetString(fname)
673  this%icsvouterout = getunit()
674  call openfile(this%icsvouterout, iout, fname, 'CSV_OUTPUT', &
675  filstat_opt='REPLACE')
676  write (iout, fmtcsvout) trim(fname), this%icsvouterout
677  !
678  ! -- create warning message
679  write (warnmsg, '(a)') &
680  'OUTER ITERATION INFORMATION WILL BE SAVED TO '//trim(fname)
681  !
682  ! -- create deprecation warning
683  call deprecation_warning('OPTIONS', 'CSV_OUTPUT', '6.1.1', &
684  warnmsg, this%parser%GetUnit())
685  else
686  write (errmsg, '(a)') 'Optional CSV_OUTPUT '// &
687  'keyword must be followed by FILEOUT'
688  call store_error(errmsg)
689  end if
690  !
691  ! -- right now these are options that are only available in the
692  ! development version and are not included in the documentation.
693  ! These options are only available when IDEVELOPMODE in
694  ! constants module is set to 1
695  case ('DEV_PTC')
696  call this%parser%DevOpt()
697  this%iallowptc = 1
698  write (iout, '(1x,A)') 'PSEUDO-TRANSIENT CONTINUATION ENABLED'
699  case ('DEV_PTC_OUTPUT')
700  call this%parser%DevOpt()
701  this%iallowptc = 1
702  call this%parser%GetStringCaps(keyword)
703  if (keyword == 'FILEOUT') then
704  call this%parser%GetString(fname)
705  if (nr_procs > 1) then
706  call append_processor_id(fname, proc_id)
707  end if
708  this%iptcout = getunit()
709  call openfile(this%iptcout, iout, fname, 'PTC-OUT', &
710  filstat_opt='REPLACE')
711  write (iout, fmtptcout) trim(fname), this%iptcout
712  else
713  write (errmsg, '(a)') &
714  'Optional PTC_OUTPUT keyword must be followed by FILEOUT'
715  call store_error(errmsg)
716  end if
717  case ('DEV_PTC_OPTION')
718  call this%parser%DevOpt()
719  this%iallowptc = 1
720  this%iptcopt = 1
721  write (iout, '(1x,A)') &
722  'PSEUDO-TRANSIENT CONTINUATION USES BNORM AND L2NORM TO '// &
723  'SET INITIAL VALUE'
724  case ('DEV_PTC_EXPONENT')
725  call this%parser%DevOpt()
726  rval = this%parser%GetDouble()
727  if (rval < dzero) then
728  write (errmsg, '(a)') 'PTC_EXPONENT must be > 0.'
729  call store_error(errmsg)
730  else
731  this%iallowptc = 1
732  this%ptcexp = rval
733  write (iout, '(1x,A,1x,g15.7)') &
734  'PSEUDO-TRANSIENT CONTINUATION EXPONENT', this%ptcexp
735  end if
736  case ('DEV_PTC_DEL0')
737  call this%parser%DevOpt()
738  rval = this%parser%GetDouble()
739  if (rval < dzero) then
740  write (errmsg, '(a)') 'IMS sln_ar: PTC_DEL0 must be > 0.'
741  call store_error(errmsg)
742  else
743  this%iallowptc = 1
744  this%ptcdel0 = rval
745  write (iout, '(1x,A,1x,g15.7)') &
746  'PSEUDO-TRANSIENT CONTINUATION INITIAL TIMESTEP', this%ptcdel0
747  end if
748  case default
749  write (errmsg, '(a,2(1x,a))') &
750  'Unknown IMS option (', trim(keyword), ').'
751  call store_error(errmsg)
752  end select
753  end do
754  write (iout, '(1x,a)') 'END OF IMS OPTIONS'
755  else
756  write (iout, '(1x,a)') 'NO IMS OPTION BLOCK DETECTED.'
757  end if
758 
759 00021 FORMAT(1x, 'SIMPLE OPTION:', /, &
760  1x, 'DEFAULT SOLVER INPUT VALUES FOR FAST SOLUTIONS')
761 00023 FORMAT(1x, 'MODERATE OPTION:', /, 1x, 'DEFAULT SOLVER', &
762  ' INPUT VALUES REFLECT MODERATELY NONLINEAR MODEL')
763 00025 FORMAT(1x, 'COMPLEX OPTION:', /, 1x, 'DEFAULT SOLVER', &
764  ' INPUT VALUES REFLECT STRONGLY NONLINEAR MODEL')
765 
766  !-------READ NONLINEAR ITERATION PARAMETERS AND LINEAR SOLVER SELECTION INDEX
767  ! -- set default nonlinear parameters
768  call this%sln_setouter(ifdparam)
769  !
770  ! -- get NONLINEAR block
771  call this%parser%GetBlock('NONLINEAR', isfound, ierr, &
772  supportopenclose=.true., blockrequired=.false.)
773  !
774  ! -- parse NONLINEAR block if detected
775  if (isfound) then
776  write (iout, '(/1x,a)') 'PROCESSING IMS NONLINEAR'
777  do
778  call this%parser%GetNextLine(endofblock)
779  if (endofblock) exit
780  call this%parser%GetStringCaps(keyword)
781  ! -- parse keyword
782  select case (keyword)
783  case ('OUTER_DVCLOSE')
784  this%dvclose = this%parser%GetDouble()
785  case ('OUTER_MAXIMUM')
786  this%mxiter = this%parser%GetInteger()
787  case ('UNDER_RELAXATION')
788  call this%parser%GetStringCaps(keyword)
789  ival = 0
790  if (keyword == 'NONE') then
791  ival = 0
792  else if (keyword == 'SIMPLE') then
793  ival = 1
794  else if (keyword == 'COOLEY') then
795  ival = 2
796  else if (keyword == 'DBD') then
797  ival = 3
798  else
799  write (errmsg, '(3a)') &
800  'Unknown UNDER_RELAXATION specified (', trim(keyword), ').'
801  call store_error(errmsg)
802  end if
803  this%nonmeth = ival
804  case ('LINEAR_SOLVER')
805  call this%parser%GetStringCaps(keyword)
806  ival = ims_solver
807  if (keyword .eq. 'DEFAULT' .or. &
808  keyword .eq. 'LINEAR') then
809  ival = ims_solver
810  else
811  write (errmsg, '(3a)') &
812  'Unknown LINEAR_SOLVER specified (', trim(keyword), ').'
813  call store_error(errmsg)
814  end if
815  this%linsolver = ival
816  case ('UNDER_RELAXATION_THETA')
817  this%theta = this%parser%GetDouble()
818  case ('UNDER_RELAXATION_KAPPA')
819  this%akappa = this%parser%GetDouble()
820  case ('UNDER_RELAXATION_GAMMA')
821  this%gamma = this%parser%GetDouble()
822  case ('UNDER_RELAXATION_MOMENTUM')
823  this%amomentum = this%parser%GetDouble()
824  case ('BACKTRACKING_NUMBER')
825  this%numtrack = this%parser%GetInteger()
826  IF (this%numtrack > 0) this%ibflag = 1
827  case ('BACKTRACKING_TOLERANCE')
828  this%btol = this%parser%GetDouble()
829  case ('BACKTRACKING_REDUCTION_FACTOR')
830  this%breduc = this%parser%GetDouble()
831  case ('BACKTRACKING_RESIDUAL_LIMIT')
832  this%res_lim = this%parser%GetDouble()
833  !
834  ! -- deprecated variables
835  case ('OUTER_HCLOSE')
836  this%dvclose = this%parser%GetDouble()
837  !
838  ! -- create warning message
839  write (warnmsg, '(a)') &
840  'SETTING OUTER_DVCLOSE TO OUTER_HCLOSE VALUE'
841  !
842  ! -- create deprecation warning
843  call deprecation_warning('NONLINEAR', 'OUTER_HCLOSE', '6.1.1', &
844  warnmsg, this%parser%GetUnit())
845  case ('OUTER_RCLOSEBND')
846  !
847  ! -- create warning message
848  write (warnmsg, '(a)') &
849  'OUTER_DVCLOSE IS USED TO EVALUATE PACKAGE CONVERGENCE'
850  !
851  ! -- create deprecation warning
852  call deprecation_warning('NONLINEAR', 'OUTER_RCLOSEBND', '6.1.1', &
853  warnmsg, this%parser%GetUnit())
854  case default
855  write (errmsg, '(3a)') &
856  'Unknown IMS NONLINEAR keyword (', trim(keyword), ').'
857  call store_error(errmsg)
858  end select
859  end do
860  write (iout, '(1x,a)') 'END OF IMS NONLINEAR DATA'
861  else
862  if (ifdparam .EQ. 0) then
863  write (errmsg, '(a)') 'NO IMS NONLINEAR block detected.'
864  call store_error(errmsg)
865  end if
866  end if
867  !
868  if (this%theta < dem3) then
869  this%theta = dem3
870  end if
871  !
872  ! -- backtracking should only be used if this%nonmeth > 0
873  if (this%nonmeth < 1) then
874  this%ibflag = 0
875  end if
876  !
877  ! -- check that MXITER is greater than zero
878  if (this%mxiter <= 0) then
879  write (errmsg, '(a)') 'Outer iteration number must be > 0.'
880  call store_error(errmsg)
881  END IF
882  !
883  ! -- write under-relaxation option
884  if (this%nonmeth > 0) then
885  WRITE (iout, *) '**UNDER-RELAXATION WILL BE USED***'
886  WRITE (iout, *)
887  elseif (this%nonmeth == 0) then
888  WRITE (iout, *) '***UNDER-RELAXATION WILL NOT BE USED***'
889  WRITE (iout, *)
890  else
891  WRITE (errmsg, '(a)') &
892  'Incorrect value for variable NONMETH was specified.'
893  call store_error(errmsg)
894  end if
895  !
896  ! -- ensure gamma is > 0 for simple
897  if (this%nonmeth == 1) then
898  if (this%gamma == 0) then
899  WRITE (errmsg, '(a)') &
900  'GAMMA must be greater than zero if SIMPLE under-relaxation is used.'
901  call store_error(errmsg)
902  end if
903  end if
904 
905  if (this%solver_mode == 'PETSC') then
906  this%linsolver = petsc_solver
907  end if
908 
909  ! configure linear settings
910  call this%linear_settings%init(this%memory_path)
911  call this%linear_settings%preset_config(ifdparam)
912  call this%linear_settings%read_from_file(this%parser, iout)
913  !
914  if (this%linear_settings%ilinmeth == cg_method) then
915  this%isymmetric = 1
916  end if
917  !
918  ! -- call secondary subroutine to initialize and read linear
919  ! solver parameters IMSLINEAR solver
920  if (this%solver_mode == "IMS") then
921  allocate (this%imslinear)
922  WRITE (iout, *) '***IMS LINEAR SOLVER WILL BE USED***'
923  call this%imslinear%imslinear_allocate(this%name, iout, this%iprims, &
924  this%mxiter, this%neq, &
925  this%system_matrix, this%rhs, &
926  this%x, this%linear_settings)
927  !
928  ! -- petsc linear solver flag
929  else if (this%solver_mode == "PETSC") then
930  call this%linear_solver%initialize(this%system_matrix, &
931  this%linear_settings, &
932  this%cnvg_summary)
933  !
934  ! -- incorrect linear solver flag
935  else
936  write (errmsg, '(a)') &
937  'Incorrect value for linear solution method specified.'
938  call store_error(errmsg)
939  end if
940  !
941  ! -- write message about matrix symmetry
942  if (this%isymmetric == 1) then
943  write (iout, '(1x,a,/)') 'A symmetric matrix will be solved'
944  else
945  write (iout, '(1x,a,/)') 'An asymmetric matrix will be solved'
946  end if
947  !
948  ! -- If CG, then go through each model and each exchange and check
949  ! for asymmetry
950  if (this%isymmetric == 1) then
951  !
952  ! -- Models
953  do i = 1, this%modellist%Count()
954  mp => getnumericalmodelfromlist(this%modellist, i)
955  if (mp%get_iasym() /= 0) then
956  write (errmsg, fmterrasym) 'MODEL', trim(adjustl(mp%name))
957  call store_error(errmsg)
958  end if
959  end do
960  !
961  ! -- Exchanges
962  do i = 1, this%exchangelist%Count()
963  cp => getnumericalexchangefromlist(this%exchangelist, i)
964  if (cp%get_iasym() /= 0) then
965  write (errmsg, fmterrasym) 'EXCHANGE', trim(adjustl(cp%name))
966  call store_error(errmsg)
967  end if
968  end do
969  !
970  end if
971  !
972  ! -- write solver data to output file
973  !
974  ! -- non-linear solver data
975  WRITE (iout, 9002) this%dvclose, this%mxiter, &
976  this%iprims, this%nonmeth, this%linsolver
977  !
978  ! -- standard outer iteration formats
979 9002 FORMAT(1x, 'OUTER ITERATION CONVERGENCE CRITERION (DVCLOSE) = ', e15.6, &
980  /1x, 'MAXIMUM NUMBER OF OUTER ITERATIONS (MXITER) = ', i0, &
981  /1x, 'SOLVER PRINTOUT INDEX (IPRIMS) = ', i0, &
982  /1x, 'NONLINEAR ITERATION METHOD (NONLINMETH) = ', i0, &
983  /1x, 'LINEAR SOLUTION METHOD (LINMETH) = ', i0)
984  !
985  if (this%nonmeth == 1) then ! simple
986  write (iout, 9003) this%gamma
987  else if (this%nonmeth == 2) then ! cooley
988  write (iout, 9004) this%gamma
989  else if (this%nonmeth == 3) then ! delta bar delta
990  write (iout, 9005) this%theta, this%akappa, this%gamma, this%amomentum
991  end if
992  !
993  ! -- write backtracking information
994  if (this%numtrack /= 0) write (iout, 9006) this%numtrack, this%btol, &
995  this%breduc, this%res_lim
996  !
997  ! -- under-relaxation formats (simple, cooley, dbd)
998 9003 FORMAT(1x, 'UNDER-RELAXATION FACTOR (GAMMA) = ', e15.6)
999 9004 FORMAT(1x, 'UNDER-RELAXATION PREVIOUS HISTORY FACTOR (GAMMA) = ', e15.6)
1000 9005 FORMAT(1x, 'UNDER-RELAXATION WEIGHT REDUCTION FACTOR (THETA) = ', e15.6, &
1001  /1x, 'UNDER-RELAXATION WEIGHT INCREASE INCREMENT (KAPPA) = ', e15.6, &
1002  /1x, 'UNDER-RELAXATION PREVIOUS HISTORY FACTOR (GAMMA) = ', e15.6, &
1003  /1x, 'UNDER-RELAXATION MOMENTUM TERM (AMOMENTUM) = ', e15.6)
1004  !
1005  ! -- backtracking formats
1006 9006 FORMAT(1x, 'MAXIMUM NUMBER OF BACKTRACKS (NUMTRACK) = ', i0, &
1007  /1x, 'BACKTRACKING TOLERANCE FACTOR (BTOL) = ', e15.6, &
1008  /1x, 'BACKTRACKING REDUCTION FACTOR (BREDUC) = ', e15.6, &
1009  /1x, 'BACKTRACKING RESIDUAL LIMIT (RES_LIM) = ', e15.6)
1010  !
1011  ! -- linear solver data
1012  if (this%linsolver == ims_solver) then
1013  call this%imslinear%imslinear_summary(this%mxiter)
1014  else
1015  call this%linear_solver%print_summary()
1016  end if
1017 
1018  ! -- write summary of solver error messages
1019  ierr = count_errors()
1020  if (ierr > 0) then
1021  call this%parser%StoreErrorUnit()
1022  end if
1023  !
1024  ! reallocate space for nonlinear arrays and initialize
1025  call mem_reallocate(this%hncg, this%mxiter, 'HNCG', this%name)
1026  call mem_reallocate(this%lrch, 3, this%mxiter, 'LRCH', this%name)
1027 
1028  ! delta-bar-delta under-relaxation
1029  if (this%nonmeth == 3) then
1030  call mem_reallocate(this%wsave, this%neq, 'WSAVE', this%name)
1031  call mem_reallocate(this%hchold, this%neq, 'HCHOLD', this%name)
1032  call mem_reallocate(this%deold, this%neq, 'DEOLD', this%name)
1033  do i = 1, this%neq
1034  this%wsave(i) = dzero
1035  this%hchold(i) = dzero
1036  this%deold(i) = dzero
1037  end do
1038  end if
1039  this%hncg = dzero
1040  this%lrch = 0
1041 
1042  ! allocate space for saving solver convergence history
1043  if (this%iprims == 2 .or. this%icsvinnerout > 0) then
1044  this%nitermax = this%linear_settings%iter1 * this%mxiter
1045  else
1046  this%nitermax = 1
1047  end if
1048 
1049  allocate (this%caccel(this%nitermax))
1050 
1051  !
1052  ! -- resize convergence report
1053  call this%cnvg_summary%reinit(this%nitermax)
1054  !
1055  ! -- check for numerical solution errors
1056  ierr = count_errors()
1057  if (ierr > 0) then
1058  call this%parser%StoreErrorUnit()
1059  end if
1060  !
1061  ! -- close ims input file
1062  call this%parser%Clear()
1063  end subroutine sln_ar
1064 
1065  !> @ brief Calculate delt
1066  !!
1067  !! Calculate time step length.
1068  !!
1069  !<
1070  subroutine sln_dt(this)
1071  ! -- modules
1072  use tdismodule, only: kstp, kper, delt
1074  use constantsmodule, only: dtwo, dthree
1075  ! -- dummy variables
1076  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1077  ! -- local variables
1078  integer(I4B) :: idir
1079  real(DP) :: delt_temp
1080  real(DP) :: fact_lower
1081  real(DP) :: fact_upper
1082  !
1083  ! -- increase or decrease delt based on kiter fraction. atsfrac should be
1084  ! a value of about 1/3. If the number of outer iterations is less than
1085  ! 1/3 of mxiter, then increase step size. If the number of outer
1086  ! iterations is greater than 2/3 of mxiter, then decrease step size.
1087  if (this%atsfrac > dzero) then
1088  delt_temp = delt
1089  fact_lower = this%mxiter * this%atsfrac
1090  fact_upper = this%mxiter - fact_lower
1091  if (this%iouttot_timestep < int(fact_lower)) then
1092  ! -- increase delt according to tsfactats
1093  idir = 1
1094  else if (this%iouttot_timestep > int(fact_upper)) then
1095  ! -- decrease delt according to tsfactats
1096  idir = -1
1097  else
1098  ! -- do not change delt
1099  idir = 0
1100  end if
1101  !
1102  ! -- submit stable dt for upcoming step
1103  call ats_submit_delt(kstp, kper, delt_temp, this%memory_path, idir=idir)
1104  end if
1105  end subroutine sln_dt
1106 
1107  !> @ brief Advance solution
1108  !!
1109  !! Advance solution.
1110  !!
1111  !<
1112  subroutine sln_ad(this)
1113  ! -- modules
1114  use tdismodule, only: kstp, kper
1115  ! -- dummy variables
1116  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1117  !
1118  ! -- write headers to CSV file
1119  if (kper == 1 .and. kstp == 1) then
1120  call this%writeCSVHeader()
1121  end if
1122 
1123  ! write PTC info on models to iout
1124  call this%writePTCInfoToFile(kper)
1125 
1126  ! reset convergence flag and inner solve counter
1127  this%icnvg = 0
1128  this%itertot_timestep = 0
1129  this%iouttot_timestep = 0
1130  end subroutine sln_ad
1131 
1132  !> @ brief Output solution
1133  !!
1134  !! Output solution data. Currently does nothing.
1135  !!
1136  !<
1137  subroutine sln_ot(this)
1138  ! -- dummy variables
1139  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1140  !
1141  ! -- Nothing to do here
1142  end subroutine sln_ot
1143 
1144  !> @ brief Finalize solution
1145  !!
1146  !! Finalize a solution.
1147  !!
1148  !<
1149  subroutine sln_fp(this)
1150  use simvariablesmodule, only: iout
1151  ! -- dummy variables
1152  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1153  !
1154  ! -- write timer output
1155  if (idevelopmode == 1) then
1156  write (iout, '(//1x,a,1x,a,1x,a)') &
1157  'Solution', trim(adjustl(this%name)), 'summary'
1158  write (iout, "(1x,70('-'))")
1159  write (iout, '(1x,a,1x,g0,1x,a)') &
1160  'Total formulate time: ', this%ttform, 'seconds'
1161  write (iout, '(1x,a,1x,g0,1x,a,/)') &
1162  'Total solution time: ', this%ttsoln, 'seconds'
1163  end if
1164  end subroutine sln_fp
1165 
1166  !> @ brief Deallocate solution
1167  !!
1168  !! Deallocate a solution.
1169  !!
1170  !<
1171  subroutine sln_da(this)
1172  ! -- modules
1174  ! -- dummy variables
1175  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1176  !
1177  ! -- IMSLinearModule
1178  if (this%linsolver == ims_solver) then
1179  call this%imslinear%imslinear_da()
1180  deallocate (this%imslinear)
1181  end if
1182  !
1183  ! -- lists
1184  call this%modellist%Clear()
1185  call this%exchangelist%Clear()
1186  deallocate (this%modellist)
1187  deallocate (this%exchangelist)
1188 
1189  call this%system_matrix%destroy()
1190  deallocate (this%system_matrix)
1191  call this%vec_x%destroy()
1192  deallocate (this%vec_x)
1193  call this%vec_rhs%destroy()
1194  deallocate (this%vec_rhs)
1195 
1196  !
1197  ! -- character arrays
1198  deallocate (this%caccel)
1199  !
1200  ! -- inner iteration table object
1201  if (associated(this%innertab)) then
1202  call this%innertab%table_da()
1203  deallocate (this%innertab)
1204  nullify (this%innertab)
1205  end if
1206  !
1207  ! -- outer iteration table object
1208  if (associated(this%outertab)) then
1209  call this%outertab%table_da()
1210  deallocate (this%outertab)
1211  nullify (this%outertab)
1212  end if
1213  !
1214  ! -- arrays
1215  call mem_deallocate(this%active)
1216  call mem_deallocate(this%xtemp)
1217  call mem_deallocate(this%dxold)
1218  call mem_deallocate(this%hncg)
1219  call mem_deallocate(this%lrch)
1220  call mem_deallocate(this%wsave)
1221  call mem_deallocate(this%hchold)
1222  call mem_deallocate(this%deold)
1223  call mem_deallocate(this%convmodstart)
1224  !
1225  ! -- convergence report
1226  call this%cnvg_summary%destroy()
1227  deallocate (this%cnvg_summary)
1228  !
1229  ! -- linear solver
1230  call this%linear_solver%destroy()
1231  deallocate (this%linear_solver)
1232  !
1233  ! -- linear solver settings
1234  call this%linear_settings%destroy()
1235  deallocate (this%linear_settings)
1236  !
1237  ! -- Scalars
1238  call mem_deallocate(this%id)
1239  call mem_deallocate(this%iu)
1240  call mem_deallocate(this%ttform)
1241  call mem_deallocate(this%ttsoln)
1242  call mem_deallocate(this%isymmetric)
1243  call mem_deallocate(this%neq)
1244  call mem_deallocate(this%matrix_offset)
1245  call mem_deallocate(this%dvclose)
1246  call mem_deallocate(this%bigchold)
1247  call mem_deallocate(this%bigch)
1248  call mem_deallocate(this%relaxold)
1249  call mem_deallocate(this%res_prev)
1250  call mem_deallocate(this%res_new)
1251  call mem_deallocate(this%icnvg)
1252  call mem_deallocate(this%itertot_timestep)
1253  call mem_deallocate(this%iouttot_timestep)
1254  call mem_deallocate(this%itertot_sim)
1255  call mem_deallocate(this%mxiter)
1256  call mem_deallocate(this%linsolver)
1257  call mem_deallocate(this%nonmeth)
1258  call mem_deallocate(this%iprims)
1259  call mem_deallocate(this%theta)
1260  call mem_deallocate(this%akappa)
1261  call mem_deallocate(this%gamma)
1262  call mem_deallocate(this%amomentum)
1263  call mem_deallocate(this%breduc)
1264  call mem_deallocate(this%btol)
1265  call mem_deallocate(this%res_lim)
1266  call mem_deallocate(this%numtrack)
1267  call mem_deallocate(this%ibflag)
1268  call mem_deallocate(this%icsvouterout)
1269  call mem_deallocate(this%icsvinnerout)
1270  call mem_deallocate(this%nitermax)
1271  call mem_deallocate(this%convnmod)
1272  call mem_deallocate(this%iallowptc)
1273  call mem_deallocate(this%iptcopt)
1274  call mem_deallocate(this%iptcout)
1275  call mem_deallocate(this%l2norm0)
1276  call mem_deallocate(this%ptcdel)
1277  call mem_deallocate(this%ptcdel0)
1278  call mem_deallocate(this%ptcexp)
1279  call mem_deallocate(this%atsfrac)
1280  end subroutine sln_da
1281 
1282  !> @ brief Solve solution
1283  !!
1284  !! Solve the models in this solution for kper and kstp.
1285  !!
1286  !<
1287  subroutine sln_ca(this, isgcnvg, isuppress_output)
1288  ! -- dummy variables
1289  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1290  integer(I4B), intent(inout) :: isgcnvg !< solution group convergence flag
1291  integer(I4B), intent(in) :: isuppress_output !< flag for suppressing output
1292  ! -- local variables
1293  class(numericalmodeltype), pointer :: mp => null()
1294  character(len=LINELENGTH) :: line
1295  character(len=LINELENGTH) :: fmt
1296  integer(I4B) :: im
1297  integer(I4B) :: kiter ! non-linear iteration counter
1298 
1299  ! advance the models, exchanges, and solution
1300  call this%prepareSolve()
1301 
1302  select case (isim_mode)
1303  case (mvalidate)
1304  line = 'mode="validation" -- Skipping matrix assembly and solution.'
1305  fmt = "(/,1x,a,/)"
1306  do im = 1, this%modellist%Count()
1307  mp => getnumericalmodelfromlist(this%modellist, im)
1308  call mp%model_message(line, fmt=fmt)
1309  end do
1310  case (mnormal)
1311  ! nonlinear iteration loop for this solution
1312  outerloop: do kiter = 1, this%mxiter
1313 
1314  ! perform a single iteration
1315  call this%solve(kiter)
1316 
1317  ! exit if converged
1318  if (this%icnvg == 1) then
1319  exit outerloop
1320  end if
1321 
1322  end do outerloop
1323 
1324  ! finish up, write convergence info, CSV file, budgets and flows, ...
1325  call this%finalizeSolve(kiter, isgcnvg, isuppress_output)
1326  end select
1327  end subroutine sln_ca
1328 
1329  !> @ brief CSV header
1330  !!
1331  !! Write header for solver output to comma-separated value files.
1332  !!
1333  !<
1334  subroutine writecsvheader(this)
1335  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1336  ! local variables
1337  integer(I4B) :: im
1338  class(numericalmodeltype), pointer :: mp => null()
1339  !
1340  ! -- outer iteration csv header
1341  if (this%icsvouterout > 0) then
1342  write (this%icsvouterout, '(*(G0,:,","))') &
1343  'total_inner_iterations', 'totim', 'kper', 'kstp', 'nouter', &
1344  'inner_iterations', 'solution_outer_dvmax', &
1345  'solution_outer_dvmax_model', 'solution_outer_dvmax_package', &
1346  'solution_outer_dvmax_node'
1347  end if
1348  !
1349  ! -- inner iteration csv header
1350  if (this%icsvinnerout > 0) then
1351  write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') &
1352  'total_inner_iterations', 'totim', 'kper', 'kstp', 'nouter', &
1353  'ninner', 'solution_inner_dvmax', 'solution_inner_dvmax_model', &
1354  'solution_inner_dvmax_node'
1355  write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') &
1356  '', 'solution_inner_rmax', 'solution_inner_rmax_model', &
1357  'solution_inner_rmax_node'
1358  ! solver items specific to ims solver
1359  if (this%linsolver == ims_solver) then
1360  write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') &
1361  '', 'solution_inner_alpha'
1362  if (this%imslinear%ilinmeth == 2) then
1363  write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') &
1364  '', 'solution_inner_omega'
1365  end if
1366  end if
1367  ! -- check for more than one model - ims only
1368  if (this%convnmod > 1 .or. simulation_mode == "PARALLEL") then
1369  do im = 1, this%modellist%Count()
1370  mp => getnumericalmodelfromlist(this%modellist, im)
1371  write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') &
1372  '', trim(adjustl(mp%name))//'_inner_dvmax', &
1373  trim(adjustl(mp%name))//'_inner_dvmax_node', &
1374  trim(adjustl(mp%name))//'_inner_rmax', &
1375  trim(adjustl(mp%name))//'_inner_rmax_node'
1376  end do
1377  end if
1378  write (this%icsvinnerout, '(a)') ''
1379  end if
1380  end subroutine writecsvheader
1381 
1382  !> @ brief PTC header
1383  !!
1384  !! Write header for pseudo-transient continuation information to a file.
1385  !!
1386  !<
1387  subroutine writeptcinfotofile(this, kper)
1388  ! -- dummy variables
1389  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1390  integer(I4B), intent(in) :: kper !< current stress period number
1391  ! -- local variable
1392  integer(I4B) :: n, im, iallowptc, iptc
1393  class(numericalmodeltype), pointer :: mp => null()
1394 
1395  ! -- determine if PTC will be used in any model
1396  n = 1
1397  do im = 1, this%modellist%Count()
1398  !
1399  ! -- set iallowptc
1400  ! -- no_ptc_option is FIRST
1401  if (this%iallowptc < 0) then
1402  if (kper > 1) then
1403  iallowptc = 1
1404  else
1405  iallowptc = 0
1406  end if
1407  ! -- no_ptc_option is ALL (0) or using PTC (1)
1408  else
1409  iallowptc = this%iallowptc
1410  end if
1411 
1412  if (iallowptc > 0) then
1413  mp => getnumericalmodelfromlist(this%modellist, im)
1414  call mp%model_ptcchk(iptc)
1415  else
1416  iptc = 0
1417  end if
1418 
1419  if (iptc /= 0) then
1420  if (n == 1) then
1421  write (iout, '(//)')
1422  n = 0
1423  end if
1424  write (iout, '(1x,a,1x,i0,1x,3a)') &
1425  'PSEUDO-TRANSIENT CONTINUATION WILL BE APPLIED TO MODEL', im, '("', &
1426  trim(adjustl(mp%name)), '") DURING THIS TIME STEP'
1427  end if
1428  end do
1429 
1430  end subroutine writeptcinfotofile
1431 
1432  !> @ brief prepare to solve
1433  !!
1434  !! Prepare for the system solve by advancing the simulation.
1435  !!
1436  !<
1437  subroutine preparesolve(this)
1438  ! -- dummy variables
1439  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1440  ! -- local variables
1441  integer(I4B) :: ic
1442  integer(I4B) :: im
1443  class(numericalexchangetype), pointer :: cp => null()
1444  class(numericalmodeltype), pointer :: mp => null()
1445 
1446  ! start timer
1447  call g_prof%start("Prepare solve"//this%id_postfix, this%tmr_prep_solve)
1448 
1449  ! synchronize for AD
1450  call this%synchronize(stg_bfr_exg_ad, this%synchronize_ctx)
1451 
1452  ! -- Exchange advance
1453  do ic = 1, this%exchangelist%Count()
1454  cp => getnumericalexchangefromlist(this%exchangelist, ic)
1455  call cp%exg_ad()
1456  end do
1457 
1458  ! -- Model advance
1459  do im = 1, this%modellist%Count()
1460  mp => getnumericalmodelfromlist(this%modellist, im)
1461  call mp%model_ad()
1462  end do
1463 
1464  ! advance solution
1465  call this%sln_ad()
1466 
1467  ! stop timer
1468  call g_prof%stop(this%tmr_prep_solve)
1469 
1470  end subroutine preparesolve
1471 
1472  !> @ brief Build and solve the simulation
1473  !!
1474  !! Builds and solve the system for this numerical solution.
1475  !! It roughly consists of the following steps
1476  !! (1) backtracking, (2) reset amat and rhs (3) calculate matrix
1477  !! terms (*_cf), (4) add coefficients to matrix (*_fc), (6) newton-raphson,
1478  !! (6) PTC, (7) linear solve, (8) convergence checks, (9) write output,
1479  !! and (10) underrelaxation
1480  !!
1481  !<
1482  subroutine solve(this, kiter)
1483  ! -- modules
1484  use tdismodule, only: kstp, kper, totim
1485  ! -- dummy variables
1486  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1487  integer(I4B), intent(in) :: kiter !< Picard iteration number
1488  ! -- local variables
1489  class(numericalmodeltype), pointer :: mp => null()
1490  class(numericalexchangetype), pointer :: cp => null()
1491  character(len=LINELENGTH) :: title
1492  character(len=LINELENGTH) :: tag
1493  character(len=LENPAKLOC) :: cmod
1494  character(len=LENPAKLOC) :: cpak
1495  character(len=LENPAKLOC) :: cpakout
1496  character(len=LENPAKLOC) :: strh
1497  character(len=25) :: cval
1498  character(len=7) :: cmsg
1499  integer(I4B) :: ic
1500  integer(I4B) :: im, m_idx, model_id
1501  integer(I4B) :: icsv0
1502  integer(I4B) :: kcsv0
1503  integer(I4B) :: ntabrows
1504  integer(I4B) :: ntabcols
1505  integer(I4B) :: i0, i1
1506  integer(I4B) :: itestmat, n
1507  integer(I4B) :: iter
1508  integer(I4B) :: inewtonur
1509  integer(I4B) :: locmax_nur
1510  integer(I4B) :: iend
1511  integer(I4B) :: icnvgmod
1512  integer(I4B) :: iptc
1513  integer(I4B) :: node_user
1514  integer(I4B) :: ipak
1515  integer(I4B) :: ipos0
1516  integer(I4B) :: ipos1
1517  real(DP) :: dxmax_nur
1518  real(DP) :: dxold_max
1519  real(DP) :: ptcf
1520  real(DP) :: ttform
1521  real(DP) :: ttsoln
1522  real(DP) :: dpak
1523  real(DP) :: outer_hncg
1524 
1525  ! start timer
1526  call g_prof%start("Solve"//this%id_postfix, this%tmr_solve)
1527 
1528  !
1529  ! -- initialize local variables
1530  icsv0 = max(1, this%itertot_sim + 1)
1531  kcsv0 = max(1, this%itertot_timestep + 1)
1532  !
1533  ! -- create header for outer iteration table
1534  if (this%iprims > 0) then
1535  if (.not. associated(this%outertab)) then
1536  !
1537  ! -- create outer iteration table
1538  ! -- table dimensions
1539  ntabrows = 1
1540  ntabcols = 6
1541  if (this%numtrack > 0) then
1542  ntabcols = ntabcols + 4
1543  end if
1544  !
1545  ! -- initialize table and define columns
1546  title = trim(this%memory_path)//' OUTER ITERATION SUMMARY'
1547  call table_cr(this%outertab, this%name, title)
1548  call this%outertab%table_df(ntabrows, ntabcols, iout, &
1549  finalize=.false.)
1550  tag = 'OUTER ITERATION STEP'
1551  call this%outertab%initialize_column(tag, 25, alignment=tableft)
1552  tag = 'OUTER ITERATION'
1553  call this%outertab%initialize_column(tag, 10, alignment=tabright)
1554  tag = 'INNER ITERATION'
1555  call this%outertab%initialize_column(tag, 10, alignment=tabright)
1556  if (this%numtrack > 0) then
1557  tag = 'BACKTRACK FLAG'
1558  call this%outertab%initialize_column(tag, 10, alignment=tabright)
1559  tag = 'BACKTRACK ITERATIONS'
1560  call this%outertab%initialize_column(tag, 10, alignment=tabright)
1561  tag = 'INCOMING RESIDUAL'
1562  call this%outertab%initialize_column(tag, 15, alignment=tabright)
1563  tag = 'OUTGOING RESIDUAL'
1564  call this%outertab%initialize_column(tag, 15, alignment=tabright)
1565  end if
1566  tag = 'MAXIMUM CHANGE'
1567  call this%outertab%initialize_column(tag, 15, alignment=tabright)
1568  tag = 'STEP SUCCESS'
1569  call this%outertab%initialize_column(tag, 7, alignment=tabright)
1570  tag = 'MAXIMUM CHANGE MODEL-(CELLID) OR MODEL-PACKAGE-(NUMBER)'
1571  call this%outertab%initialize_column(tag, 34, alignment=tabright)
1572  end if
1573  end if
1574  !
1575  ! -- backtracking
1576  if (this%numtrack > 0) then
1577  call this%sln_backtracking(mp, cp, kiter)
1578  end if
1579  !
1580  call code_timer(0, ttform, this%ttform)
1581  call g_prof%start("Formulate", this%tmr_formulate)
1582  !
1583  ! -- (re)build the solution matrix
1584  call this%sln_buildsystem(kiter, inewton=1)
1585  !
1586  ! -- Calculate pseudo-transient continuation factor for each model
1587  call this%sln_calc_ptc(iptc, ptcf)
1588  !
1589  ! -- Add model Newton-Raphson terms to solution
1590  do im = 1, this%modellist%Count()
1591  mp => getnumericalmodelfromlist(this%modellist, im)
1592  call mp%model_nr(kiter, this%system_matrix, 1)
1593  end do
1594  call code_timer(1, ttform, this%ttform)
1595  call g_prof%stop(this%tmr_formulate)
1596  !
1597  ! -- linear solve
1598  call code_timer(0, ttsoln, this%ttsoln)
1599  call g_prof%start("Linear solve", this%tmr_linsolve)
1600  call this%sln_ls(kiter, kstp, kper, iter, iptc, ptcf)
1601  call g_prof%stop(this%tmr_linsolve)
1602  call code_timer(1, ttsoln, this%ttsoln)
1603  !
1604  ! -- increment counters storing the total number of linear and
1605  ! non-linear iterations for this timestep and the total
1606  ! number of linear iterations for all timesteps
1607  this%itertot_timestep = this%itertot_timestep + iter
1608  this%iouttot_timestep = this%iouttot_timestep + 1
1609  this%itertot_sim = this%itertot_sim + iter
1610  !
1611  ! -- save matrix to a file
1612  ! to enable set itestmat to 1 and recompile
1613  !-------------------------------------------------------
1614  itestmat = 0
1615  if (itestmat /= 0) then
1616  open (99, file='sol_MF6.TXT')
1617  WRITE (99, *) 'MATRIX SOLUTION FOLLOWS'
1618  WRITE (99, '(10(I8,G15.4))') (n, this%x(n), n=1, this%NEQ)
1619  close (99)
1620  call pstop()
1621  end if
1622  !-------------------------------------------------------
1623  !
1624  ! -- check convergence of solution
1625  call this%sln_get_dxmax(this%hncg(kiter), this%lrch(1, kiter))
1626  this%icnvg = 0
1627  if (this%sln_has_converged(this%hncg(kiter))) then
1628  this%icnvg = 1
1629  end if
1630  !
1631  ! -- set failure flag
1632  if (this%icnvg == 0) then
1633  cmsg = ' '
1634  else
1635  cmsg = '*'
1636  end if
1637  !
1638  ! -- set flag if this is the last outer iteration
1639  iend = 0
1640  if (kiter == this%mxiter) then
1641  iend = 1
1642  end if
1643  !
1644  ! -- write maximum dependent-variable change from linear solver to list file
1645  if (this%iprims > 0) then
1646  cval = 'Model'
1647  call this%sln_get_loc(this%lrch(1, kiter), strh)
1648  !
1649  ! -- add data to outertab
1650  call this%outertab%add_term(cval)
1651  call this%outertab%add_term(kiter)
1652  call this%outertab%add_term(iter)
1653  if (this%numtrack > 0) then
1654  call this%outertab%add_term(' ')
1655  call this%outertab%add_term(' ')
1656  call this%outertab%add_term(' ')
1657  call this%outertab%add_term(' ')
1658  end if
1659  call this%outertab%add_term(this%hncg(kiter))
1660  call this%outertab%add_term(cmsg)
1661  call this%outertab%add_term(trim(strh))
1662  end if
1663  !
1664  ! -- Additional convergence check for exchanges
1665  do ic = 1, this%exchangelist%Count()
1666  cp => getnumericalexchangefromlist(this%exchangelist, ic)
1667  call cp%exg_cc(this%icnvg)
1668  end do
1669  !
1670  ! -- additional convergence check for model packages
1671  icnvgmod = this%icnvg
1672  cpak = ' '
1673  ipak = 0
1674  dpak = dzero
1675  do im = 1, this%modellist%Count()
1676  mp => getnumericalmodelfromlist(this%modellist, im)
1677  call mp%get_mcellid(0, cmod)
1678  call mp%model_cc(this%itertot_sim, kiter, iend, icnvgmod, &
1679  cpak, ipak, dpak)
1680  if (ipak /= 0) then
1681  ipos0 = index(cpak, '-', back=.true.)
1682  ipos1 = len_trim(cpak)
1683  write (cpakout, '(a,a,"-(",i0,")",a)') &
1684  trim(cmod), cpak(1:ipos0 - 1), ipak, cpak(ipos0:ipos1)
1685  else
1686  cpakout = ' '
1687  end if
1688  end do
1689  !
1690  ! -- evaluate package convergence - only done if convergence is achieved
1691  if (this%icnvg == 1) then
1692  this%icnvg = this%sln_package_convergence(dpak, cpakout, iend)
1693  !
1694  ! -- write maximum change in package convergence check
1695  if (this%iprims > 0) then
1696  cval = 'Package'
1697  if (this%icnvg /= 1) then
1698  cmsg = ' '
1699  else
1700  cmsg = '*'
1701  end if
1702  if (len_trim(cpakout) > 0) then
1703  !
1704  ! -- add data to outertab
1705  call this%outertab%add_term(cval)
1706  call this%outertab%add_term(kiter)
1707  call this%outertab%add_term(' ')
1708  if (this%numtrack > 0) then
1709  call this%outertab%add_term(' ')
1710  call this%outertab%add_term(' ')
1711  call this%outertab%add_term(' ')
1712  call this%outertab%add_term(' ')
1713  end if
1714  call this%outertab%add_term(dpak)
1715  call this%outertab%add_term(cmsg)
1716  call this%outertab%add_term(cpakout)
1717  end if
1718  end if
1719  end if
1720  !
1721  ! -- under-relaxation - only done if convergence not achieved
1722  if (this%icnvg /= 1) then
1723  if (this%nonmeth > 0) then
1724  call this%sln_underrelax(kiter, this%hncg(kiter), this%neq, &
1725  this%active, this%x, this%xtemp)
1726  else
1727  call this%sln_calcdx(this%neq, this%active, &
1728  this%x, this%xtemp, this%dxold)
1729  end if
1730  !
1731  ! -- adjust heads by newton under-relaxation, if necessary
1732  inewtonur = 0
1733  dxmax_nur = dzero
1734  locmax_nur = 0
1735  do im = 1, this%modellist%Count()
1736  mp => getnumericalmodelfromlist(this%modellist, im)
1737  i0 = mp%moffset + 1 - this%matrix_offset
1738  i1 = i0 + mp%neq - 1
1739  call mp%model_nur(mp%neq, this%x(i0:i1), this%xtemp(i0:i1), &
1740  this%dxold(i0:i1), inewtonur, dxmax_nur, locmax_nur)
1741  end do
1742  !
1743  ! -- synchronize Newton Under-relaxation flag
1744  inewtonur = this%sln_sync_newtonur_flag(inewtonur)
1745  !
1746  ! -- check for convergence if newton under-relaxation applied
1747  if (inewtonur /= 0) then
1748  !
1749  ! -- calculate maximum change in heads in cells that have
1750  ! not been adjusted by newton under-relxation
1751  call this%sln_maxval(this%neq, this%dxold, dxold_max)
1752  !
1753  ! -- evaluate convergence
1754  if (this%sln_nur_has_converged(dxold_max, this%hncg(kiter))) then
1755  !
1756  ! -- converged
1757  this%icnvg = 1
1758  !
1759  ! -- reset outer dependent-variable change and location for output
1760  call this%sln_get_dxmax(this%hncg(kiter), this%lrch(1, kiter))
1761  !
1762  ! -- write revised dependent-variable change data after
1763  ! newton under-relaxation
1764  if (this%iprims > 0) then
1765  cval = 'Newton under-relaxation'
1766  cmsg = '*'
1767  call this%sln_get_loc(this%lrch(1, kiter), strh)
1768  !
1769  ! -- add data to outertab
1770  call this%outertab%add_term(cval)
1771  call this%outertab%add_term(kiter)
1772  call this%outertab%add_term(iter)
1773  if (this%numtrack > 0) then
1774  call this%outertab%add_term(' ')
1775  call this%outertab%add_term(' ')
1776  call this%outertab%add_term(' ')
1777  call this%outertab%add_term(' ')
1778  end if
1779  call this%outertab%add_term(this%hncg(kiter))
1780  call this%outertab%add_term(cmsg)
1781  call this%outertab%add_term(trim(strh))
1782  end if
1783  end if
1784  end if
1785  end if
1786  !
1787  ! -- write to outer iteration csv file
1788  if (this%icsvouterout > 0) then
1789  !
1790  ! -- set outer dependent-variable change variable
1791  outer_hncg = this%hncg(kiter)
1792  !
1793  ! -- model convergence error
1794  if (abs(outer_hncg) > abs(dpak)) then
1795  !
1796  ! -- get model number and user node number
1797  call this%sln_get_nodeu(this%lrch(1, kiter), m_idx, node_user)
1798  cpakout = ''
1799  else if (outer_hncg == dzero .and. dpak == dzero) then ! zero change, location could be any
1800  m_idx = 0
1801  node_user = 0
1802  !
1803  ! -- then it's a package convergence error
1804  else
1805  !
1806  ! -- set convergence error, model number, user node number,
1807  ! and package name
1808  outer_hncg = dpak
1809  ipos0 = index(cmod, '_')
1810  read (cmod(1:ipos0 - 1), *) m_idx
1811  node_user = ipak
1812  ipos0 = index(cpak, '-', back=.true.)
1813  cpakout = cpak(1:ipos0 - 1)
1814  end if
1815  !
1816  ! -- write line to outer iteration csv file
1817  if (m_idx > 0) then
1818  mp => getnumericalmodelfromlist(this%modellist, m_idx) ! TODO_MJR: right list?
1819  model_id = mp%id
1820  else
1821  model_id = 0
1822  end if
1823  write (this%icsvouterout, '(*(G0,:,","))') &
1824  this%itertot_sim, totim, kper, kstp, kiter, iter, &
1825  outer_hncg, model_id, trim(cpakout), node_user
1826  end if
1827  !
1828  ! -- write to inner iteration csv file
1829  if (this%icsvinnerout > 0) then
1830  call this%csv_convergence_summary(this%icsvinnerout, totim, kper, kstp, &
1831  kiter, iter, icsv0, kcsv0)
1832  end if
1833 
1834  ! stop timer
1835  call g_prof%stop(this%tmr_solve)
1836 
1837  end subroutine solve
1838 
1839  !> @ brief finalize a solution
1840  !!
1841  !! Finalize the solution. Called after the outer iteration loop.
1842  !!
1843  !<
1844  subroutine finalizesolve(this, kiter, isgcnvg, isuppress_output)
1845  ! -- modules
1846  use tdismodule, only: kper, kstp
1847  ! -- dummy variables
1848  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
1849  integer(I4B), intent(in) :: kiter !< Picard iteration number after convergence or failure
1850  integer(I4B), intent(inout) :: isgcnvg !< solution group convergence flag
1851  integer(I4B), intent(in) :: isuppress_output !< flag for suppressing output
1852  ! -- local variables
1853  integer(I4B) :: ic, im
1854  class(numericalmodeltype), pointer :: mp => null()
1855  class(numericalexchangetype), pointer :: cp => null()
1856  ! -- formats for convergence info
1857  character(len=*), parameter :: fmtnocnvg = &
1858  "(1X,'Solution ', i0, ' did not converge for stress period ', i0, &
1859  &' and time step ', i0)"
1860  character(len=*), parameter :: fmtcnvg = &
1861  "(1X, I0, ' CALLS TO NUMERICAL SOLUTION ', 'IN TIME STEP ', I0, &
1862  &' STRESS PERIOD ',I0,/1X,I0,' TOTAL ITERATIONS')"
1863 
1864  ! start timer
1865  call g_prof%start("Finalize solve"//this%id_postfix, this%tmr_final_solve)
1866 
1867  !
1868  ! -- finalize the outer iteration table
1869  if (this%iprims > 0) then
1870  call this%outertab%finalize_table()
1871  end if
1872  !
1873  ! -- write convergence info
1874  !
1875  ! -- convergence was achieved
1876  if (this%icnvg /= 0) then
1877  if (this%iprims > 0) then
1878  write (iout, fmtcnvg) kiter, kstp, kper, this%itertot_timestep
1879  end if
1880  !
1881  ! -- convergence was not achieved
1882  else
1883  write (iout, fmtnocnvg) this%id, kper, kstp
1884  end if
1885  !
1886  ! -- write inner iteration convergence summary
1887  if (this%iprims == 2) then
1888  !
1889  ! -- write summary for each model
1890  do im = 1, this%modellist%Count()
1891  mp => getnumericalmodelfromlist(this%modellist, im)
1892  call this%convergence_summary(mp%iout, im, this%itertot_timestep)
1893  end do
1894  !
1895  ! -- write summary for entire solution
1896  call this%convergence_summary(iout, this%convnmod + 1, &
1897  this%itertot_timestep)
1898  end if
1899  !
1900  ! -- set solution group convergence flag
1901  if (this%icnvg == 0) isgcnvg = 0
1902 
1903  call g_prof%start("Calculate flows", this%tmr_flows)
1904 
1905  !
1906  ! -- Calculate flow for each model
1907  do im = 1, this%modellist%Count()
1908  mp => getnumericalmodelfromlist(this%modellist, im)
1909  call mp%model_cq(this%icnvg, isuppress_output)
1910  end do
1911  !
1912  ! -- Calculate flow for each exchange
1913  do ic = 1, this%exchangelist%Count()
1914  cp => getnumericalexchangefromlist(this%exchangelist, ic)
1915  call cp%exg_cq(isgcnvg, isuppress_output, this%id)
1916  end do
1917 
1918  call g_prof%stop(this%tmr_flows)
1919  call g_prof%start("Calculate budgets", this%tmr_budgets)
1920 
1921  !
1922  ! -- Budget terms for each model
1923  do im = 1, this%modellist%Count()
1924  mp => getnumericalmodelfromlist(this%modellist, im)
1925  call mp%model_bd(this%icnvg, isuppress_output)
1926  end do
1927  !
1928  ! -- Budget terms for each exchange
1929  do ic = 1, this%exchangelist%Count()
1930  cp => getnumericalexchangefromlist(this%exchangelist, ic)
1931  call cp%exg_bd(isgcnvg, isuppress_output, this%id)
1932  end do
1933 
1934  ! stop timer
1935  call g_prof%stop(this%tmr_budgets)
1936  call g_prof%stop(this%tmr_final_solve)
1937 
1938  end subroutine finalizesolve
1939 
1940  ! helper routine to calculate coefficients and setup the solution matrix
1941  subroutine sln_buildsystem(this, kiter, inewton)
1942  class(numericalsolutiontype) :: this
1943  integer(I4B), intent(in) :: kiter
1944  integer(I4B), intent(in) :: inewton
1945  ! local
1946  integer(I4B) :: im, ic
1947  class(numericalmodeltype), pointer :: mp
1948  class(numericalexchangetype), pointer :: cp
1949  !
1950  ! -- Set amat and rhs to zero
1951  call this%sln_reset()
1952 
1953  ! reset models
1954  do im = 1, this%modellist%Count()
1955  mp => getnumericalmodelfromlist(this%modellist, im)
1956  call mp%model_reset()
1957  end do
1958 
1959  ! synchronize for CF
1960  call this%synchronize(stg_bfr_exg_cf, this%synchronize_ctx)
1961 
1962  !
1963  ! -- Calculate the matrix terms for each exchange
1964  do ic = 1, this%exchangelist%Count()
1965  cp => getnumericalexchangefromlist(this%exchangelist, ic)
1966  call cp%exg_cf(kiter)
1967  end do
1968  !
1969  ! -- Calculate the matrix terms for each model
1970  do im = 1, this%modellist%Count()
1971  mp => getnumericalmodelfromlist(this%modellist, im)
1972  call mp%model_cf(kiter)
1973  end do
1974 
1975  ! synchronize for FC
1976  call this%synchronize(stg_bfr_exg_fc, this%synchronize_ctx)
1977 
1978  !
1979  ! -- Add exchange coefficients to the solution
1980  do ic = 1, this%exchangelist%Count()
1981  cp => getnumericalexchangefromlist(this%exchangelist, ic)
1982  call cp%exg_fc(kiter, this%system_matrix, this%rhs, inewton)
1983  end do
1984  !
1985  ! -- Add model coefficients to the solution
1986  do im = 1, this%modellist%Count()
1987  mp => getnumericalmodelfromlist(this%modellist, im)
1988  call mp%model_fc(kiter, this%system_matrix, inewton)
1989  end do
1990 
1991  end subroutine sln_buildsystem
1992 
1993  !> @ brief Solution convergence summary
1994  !!
1995  !! Save convergence summary to a File.
1996  !!
1997  !<
1998  subroutine convergence_summary(this, iu, im, itertot_timestep)
1999  ! -- modules
2000  use inputoutputmodule, only: getunit
2001  ! -- dummy variables
2002  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2003  integer(I4B), intent(in) :: iu !< file unit number for summary file
2004  integer(I4B), intent(in) :: im !< model number
2005  integer(I4B), intent(in) :: itertot_timestep !< total iteration for the time step
2006  ! -- local variables
2007  character(len=LINELENGTH) :: title
2008  character(len=LINELENGTH) :: tag
2009  character(len=LENPAKLOC) :: loc_dvmax_str
2010  character(len=LENPAKLOC) :: loc_rmax_str
2011  integer(I4B) :: ntabrows
2012  integer(I4B) :: ntabcols
2013  integer(I4B) :: iinner
2014  integer(I4B) :: i0
2015  integer(I4B) :: iouter
2016  integer(I4B) :: j
2017  integer(I4B) :: k
2018  integer(I4B) :: locdv
2019  integer(I4B) :: locdr
2020  real(DP) :: dv !< the maximum change in the dependent variable
2021  real(DP) :: res !< the maximum value of the residual vector
2022  !
2023  ! -- initialize local variables
2024  loc_dvmax_str = ''
2025  loc_rmax_str = ''
2026  iouter = 1
2027  locdv = 0
2028  locdr = 0
2029  !
2030  ! -- initialize inner iteration summary table
2031  if (.not. associated(this%innertab)) then
2032  !
2033  ! -- create outer iteration table
2034  ! -- table dimensions
2035  ntabrows = itertot_timestep
2036  ntabcols = 7
2037  !
2038  ! -- initialize table and define columns
2039  title = trim(this%memory_path)//' INNER ITERATION SUMMARY'
2040  call table_cr(this%innertab, this%name, title)
2041  call this%innertab%table_df(ntabrows, ntabcols, iu)
2042  tag = 'TOTAL ITERATION'
2043  call this%innertab%initialize_column(tag, 10, alignment=tabright)
2044  tag = 'OUTER ITERATION'
2045  call this%innertab%initialize_column(tag, 10, alignment=tabright)
2046  tag = 'INNER ITERATION'
2047  call this%innertab%initialize_column(tag, 10, alignment=tabright)
2048  tag = 'MAXIMUM CHANGE'
2049  call this%innertab%initialize_column(tag, 15, alignment=tabright)
2050  tag = 'MAXIMUM CHANGE MODEL-(CELLID)'
2051  call this%innertab%initialize_column(tag, lenpakloc, alignment=tabright)
2052  tag = 'MAXIMUM RESIDUAL'
2053  call this%innertab%initialize_column(tag, 15, alignment=tabright)
2054  tag = 'MAXIMUM RESIDUAL MODEL-(CELLID)'
2055  call this%innertab%initialize_column(tag, lenpakloc, alignment=tabright)
2056  !
2057  ! -- reset the output unit and the number of rows (maxbound)
2058  else
2059  call this%innertab%set_maxbound(itertot_timestep)
2060  call this%innertab%set_iout(iu)
2061  end if
2062  !
2063  ! -- write the inner iteration summary to unit iu
2064  i0 = 0
2065  do k = 1, itertot_timestep
2066  iinner = this%cnvg_summary%itinner(k)
2067  if (iinner <= i0) then
2068  iouter = iouter + 1
2069  end if
2070  if (im > this%convnmod) then
2071  dv = dzero
2072  res = dzero
2073  do j = 1, this%convnmod
2074  if (abs(this%cnvg_summary%convdvmax(j, k)) > abs(dv)) then
2075  locdv = this%cnvg_summary%convlocdv(j, k)
2076  dv = this%cnvg_summary%convdvmax(j, k)
2077  end if
2078  if (abs(this%cnvg_summary%convrmax(j, k)) > abs(res)) then
2079  locdr = this%cnvg_summary%convlocr(j, k)
2080  res = this%cnvg_summary%convrmax(j, k)
2081  end if
2082  end do
2083  else
2084  locdv = this%cnvg_summary%convlocdv(im, k)
2085  locdr = this%cnvg_summary%convlocr(im, k)
2086  dv = this%cnvg_summary%convdvmax(im, k)
2087  res = this%cnvg_summary%convrmax(im, k)
2088  end if
2089  call this%sln_get_loc(locdv, loc_dvmax_str)
2090  call this%sln_get_loc(locdr, loc_rmax_str)
2091  !
2092  ! -- add data to innertab
2093  call this%innertab%add_term(k)
2094  call this%innertab%add_term(iouter)
2095  call this%innertab%add_term(iinner)
2096  call this%innertab%add_term(dv)
2097  call this%innertab%add_term(adjustr(trim(loc_dvmax_str)))
2098  call this%innertab%add_term(res)
2099  call this%innertab%add_term(adjustr(trim(loc_rmax_str)))
2100  !
2101  ! -- update i0
2102  i0 = iinner
2103  end do
2104  end subroutine convergence_summary
2105 
2106  !> @ brief Solution convergence CSV summary
2107  !!
2108  !! Save convergence summary to a comma-separated value file.
2109  !!
2110  !<
2111  subroutine csv_convergence_summary(this, iu, totim, kper, kstp, kouter, &
2112  niter, istart, kstart)
2113  ! -- modules
2114  use inputoutputmodule, only: getunit
2115  ! -- dummy variables
2116  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2117  integer(I4B), intent(in) :: iu !< file unit number
2118  real(DP), intent(in) :: totim !< total simulation time
2119  integer(I4B), intent(in) :: kper !< stress period number
2120  integer(I4B), intent(in) :: kstp !< time step number
2121  integer(I4B), intent(in) :: kouter !< number of outer (Picard) iterations
2122  integer(I4B), intent(in) :: niter !< number of inner iteration in this time step
2123  integer(I4B), intent(in) :: istart !< starting iteration number for this time step
2124  integer(I4B), intent(in) :: kstart !< starting position in the conv* arrays
2125  ! -- local
2126  integer(I4B) :: itot
2127  integer(I4B) :: m_idx, j, k
2128  integer(I4B) :: kpos
2129  integer(I4B) :: loc_dvmax !< solution node number (row) of max. dep. var. change
2130  integer(I4B) :: loc_rmax !< solution node number (row) of max. residual
2131  integer(I4B) :: model_id, node_user
2132  real(DP) :: dvmax !< maximum dependent variable change
2133  real(DP) :: rmax !< maximum residual
2134  class(numericalmodeltype), pointer :: num_mod => null()
2135  !
2136  ! -- initialize local variables
2137  itot = istart
2138  !
2139  ! -- write inner iteration results to the inner csv output file
2140  do k = 1, niter
2141  kpos = kstart + k - 1
2142  write (iu, '(*(G0,:,","))', advance='NO') &
2143  itot, totim, kper, kstp, kouter, k
2144  !
2145  ! -- solution summary
2146  dvmax = dzero
2147  rmax = dzero
2148  do j = 1, this%convnmod
2149  if (abs(this%cnvg_summary%convdvmax(j, kpos)) > abs(dvmax)) then
2150  loc_dvmax = this%cnvg_summary%convlocdv(j, kpos)
2151  dvmax = this%cnvg_summary%convdvmax(j, kpos)
2152  end if
2153  if (abs(this%cnvg_summary%convrmax(j, kpos)) > abs(rmax)) then
2154  loc_rmax = this%cnvg_summary%convlocr(j, kpos)
2155  rmax = this%cnvg_summary%convrmax(j, kpos)
2156  end if
2157  end do
2158  !
2159  ! -- no change, could be anywhere
2160  if (dvmax == dzero) loc_dvmax = 0
2161  if (rmax == dzero) loc_rmax = 0
2162  !
2163  ! -- get model number and user node number for max. dep. var. change
2164  if (loc_dvmax > 0) then
2165  call this%sln_get_nodeu(loc_dvmax, m_idx, node_user)
2166  num_mod => getnumericalmodelfromlist(this%modellist, m_idx)
2167  model_id = num_mod%id
2168  else
2169  model_id = 0
2170  node_user = 0
2171  end if
2172  write (iu, '(*(G0,:,","))', advance='NO') '', dvmax, model_id, node_user
2173  !
2174  ! -- get model number and user node number for max. residual
2175  if (loc_rmax > 0) then
2176  call this%sln_get_nodeu(loc_rmax, m_idx, node_user)
2177  num_mod => getnumericalmodelfromlist(this%modellist, m_idx)
2178  model_id = num_mod%id
2179  else
2180  model_id = 0
2181  node_user = 0
2182  end if
2183  write (iu, '(*(G0,:,","))', advance='NO') '', rmax, model_id, node_user
2184  !
2185  ! -- write ims acceleration parameters
2186  if (this%linsolver == ims_solver) then
2187  write (iu, '(*(G0,:,","))', advance='NO') &
2188  '', trim(adjustl(this%caccel(kpos)))
2189  end if
2190  !
2191  ! -- write information for each model
2192  if (this%convnmod > 1 .or. simulation_mode == "PARALLEL") then
2193  do j = 1, this%cnvg_summary%convnmod
2194  loc_dvmax = this%cnvg_summary%convlocdv(j, kpos)
2195  dvmax = this%cnvg_summary%convdvmax(j, kpos)
2196  loc_rmax = this%cnvg_summary%convlocr(j, kpos)
2197  rmax = this%cnvg_summary%convrmax(j, kpos)
2198  !
2199  ! -- get model number and user node number for max. dep. var. change
2200  if (loc_dvmax > 0) then
2201  call this%sln_get_nodeu(loc_dvmax, m_idx, node_user)
2202  else
2203  node_user = 0
2204  end if
2205  write (iu, '(*(G0,:,","))', advance='NO') '', dvmax, node_user
2206  !
2207  ! -- get model number and user node number for max. residual
2208  if (loc_rmax > 0) then
2209  call this%sln_get_nodeu(loc_rmax, m_idx, node_user)
2210  else
2211  node_user = 0
2212  end if
2213  write (iu, '(*(G0,:,","))', advance='NO') '', rmax, node_user
2214  end do
2215  end if
2216  !
2217  ! -- write line
2218  write (iu, '(a)') ''
2219  !
2220  ! -- update itot
2221  itot = itot + 1
2222  end do
2223  !
2224  ! -- flush file
2225  flush (iu)
2226  end subroutine csv_convergence_summary
2227 
2228  !> @ brief Save solution data to a file
2229  !!
2230  !! Save solution ia vector, ja vector , coefficient matrix, right-hand side
2231  !! vector, and the dependent-variable vector to a file.
2232  !!
2233  !<
2234  subroutine save(this, filename)
2235  use sparsematrixmodule
2236  ! -- modules
2237  use inputoutputmodule, only: getunit
2238  ! -- dummy variables
2239  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2240  character(len=*), intent(in) :: filename !< filename to save solution data
2241  ! -- local variables
2242  integer(I4B) :: inunit
2243  !
2244  select type (spm => this%system_matrix)
2245  class is (sparsematrixtype)
2246  inunit = getunit()
2247  open (unit=inunit, file=filename, status='unknown')
2248  write (inunit, *) 'ia'
2249  write (inunit, *) spm%ia
2250  write (inunit, *) 'ja'
2251  write (inunit, *) spm%ja
2252  write (inunit, *) 'amat'
2253  write (inunit, *) spm%amat
2254  write (inunit, *) 'rhs'
2255  write (inunit, *) this%rhs
2256  write (inunit, *) 'x'
2257  write (inunit, *) this%x
2258  close (inunit)
2259  end select
2260  end subroutine save
2261 
2262  !> @ brief Add a model
2263  !!
2264  !! Add a model to solution.
2265  !!
2266  !<
2267  subroutine add_model(this, mp)
2268  ! -- dummy variables
2269  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2270  class(basemodeltype), pointer, intent(in) :: mp !< model instance
2271  ! -- local variables
2272  class(numericalmodeltype), pointer :: m => null()
2273  !
2274  ! -- add a model
2275  select type (mp)
2276  class is (numericalmodeltype)
2277  m => mp
2278  call addnumericalmodeltolist(this%modellist, m)
2279  end select
2280  end subroutine add_model
2281 
2282  !> @brief Get a list of models
2283  !!
2284  !! Returns a pointer to the list of models in this solution.
2285  !!
2286  !<
2287  function get_models(this) result(models)
2288  ! -- return variable
2289  type(listtype), pointer :: models !< pointer to the model list
2290  ! -- dummy variables
2291  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2292 
2293  models => this%modellist
2294 
2295  end function get_models
2296 
2297  !> @brief Add exchange
2298  !!
2299  !! Add and exchange to this%exchangelist.
2300  !!
2301  !<
2302  subroutine add_exchange(this, exchange)
2303  ! -- dummy variables
2304  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2305  class(baseexchangetype), pointer, intent(in) :: exchange !< model exchange instance
2306  ! -- local variables
2307  class(numericalexchangetype), pointer :: num_ex => null()
2308  !
2309  ! -- add exchange
2310  select type (exchange)
2311  class is (numericalexchangetype)
2312  num_ex => exchange
2313  call addnumericalexchangetolist(this%exchangelist, num_ex)
2314  end select
2315  end subroutine add_exchange
2316 
2317  !> @brief Returns a pointer to the list of exchanges in this solution
2318  !<
2319  function get_exchanges(this) result(exchanges)
2320  class(numericalsolutiontype) :: this !< instance of the numerical solution
2321  type(listtype), pointer :: exchanges !< pointer to the exchange list
2322 
2323  exchanges => this%exchangelist
2324 
2325  end function get_exchanges
2326 
2327  !> @ brief Assign solution connections
2328  !!
2329  !! Assign solution connections. This is the main workhorse method for a
2330  !! solution. The method goes through all the models and all the connections
2331  !! and builds up the sparse matrix. Steps are (1) add internal model
2332  !! connections, (2) add cross terms, (3) allocate solution arrays, (4) create
2333  !! mapping arrays, and (5) fill cross term values if necessary.
2334  !!
2335  !<
2336  subroutine sln_connect(this)
2337  ! -- modules
2339  ! -- dummy variables
2340  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2341  ! -- local variables
2342  class(numericalmodeltype), pointer :: mp => null()
2343  class(numericalexchangetype), pointer :: cp => null()
2344  integer(I4B) :: im
2345  integer(I4B) :: ic
2346  !
2347  ! -- Add internal model connections to sparse
2348  do im = 1, this%modellist%Count()
2349  mp => getnumericalmodelfromlist(this%modellist, im)
2350  call mp%model_ac(this%sparse)
2351  end do
2352  !
2353  ! -- synchronize before AC
2354  call this%synchronize(stg_bfr_exg_ac, this%synchronize_ctx)
2355  !
2356  ! -- Add the cross terms to sparse
2357  do ic = 1, this%exchangelist%Count()
2358  cp => getnumericalexchangefromlist(this%exchangelist, ic)
2359  call cp%exg_ac(this%sparse)
2360  end do
2361  !
2362  ! -- The number of non-zero array values are now known so
2363  ! -- ia and ja can be created from sparse. then destroy sparse
2364  call this%sparse%sort()
2365  call this%system_matrix%init(this%sparse, this%name)
2366  call this%sparse%destroy()
2367  !
2368  ! -- Create mapping arrays for each model. Mapping assumes
2369  ! -- that each row has the diagonal in the first position,
2370  ! -- however, rows do not need to be sorted.
2371  do im = 1, this%modellist%Count()
2372  mp => getnumericalmodelfromlist(this%modellist, im)
2373  call mp%model_mc(this%system_matrix)
2374  end do
2375  !
2376  ! -- Create arrays for mapping exchange connections to global solution
2377  do ic = 1, this%exchangelist%Count()
2378  cp => getnumericalexchangefromlist(this%exchangelist, ic)
2379  call cp%exg_mc(this%system_matrix)
2380  end do
2381  end subroutine sln_connect
2382 
2383  !> @ brief Reset the solution
2384  !!
2385  !! Reset the solution by setting the coefficient matrix and right-hand side
2386  !! vectors to zero.
2387  !!
2388  !<
2389  subroutine sln_reset(this)
2390  ! -- dummy variables
2391  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2392  !
2393  ! -- reset the solution
2394  call this%system_matrix%zero_entries()
2395  call this%vec_rhs%zero_entries()
2396  end subroutine sln_reset
2397 
2398  !> @ brief Solve the linear system of equations
2399  !!
2400  !! Solve the linear system of equations. Steps include (1) matrix cleanup,
2401  !! (2) add pseudo-transient continuation terms, and (3) residual reduction.
2402  !!
2403  !<
2404  subroutine sln_ls(this, kiter, kstp, kper, in_iter, iptc, ptcf)
2405  ! -- dummy variables
2406  class(numericalsolutiontype), intent(inout) :: this !< NumericalSolutionType instance
2407  integer(I4B), intent(in) :: kiter
2408  integer(I4B), intent(in) :: kstp
2409  integer(I4B), intent(in) :: kper
2410  integer(I4B), intent(inout) :: in_iter
2411  integer(I4B), intent(inout) :: iptc
2412  real(DP), intent(in) :: ptcf
2413  ! -- local variables
2414  logical(LGP) :: lsame
2415  integer(I4B) :: ieq
2416  integer(I4B) :: irow_glo
2417  integer(I4B) :: itestmat
2418  integer(I4B) :: ipos
2419  integer(I4B) :: icol_s
2420  integer(I4B) :: icol_e
2421  integer(I4B) :: jcol
2422  integer(I4B) :: iptct
2423  integer(I4B) :: iallowptc
2424  real(DP) :: adiag
2425  real(DP) :: diagval
2426  real(DP) :: l2norm
2427  real(DP) :: ptcval
2428  real(DP) :: bnorm
2429  character(len=50) :: fname
2430  character(len=*), parameter :: fmtfname = "('mf6mat_', i0, '_', i0, &
2431  &'_', i0, '_', i0, '.txt')"
2432  !
2433  ! -- take care of loose ends for all nodes before call to solver
2434  do ieq = 1, this%neq
2435  !
2436  ! -- get (global) cell id
2437  irow_glo = ieq + this%matrix_offset
2438  !
2439  ! -- store x in temporary location
2440  this%xtemp(ieq) = this%x(ieq)
2441  !
2442  ! -- make adjustments to the continuity equation for the node
2443  ! -- adjust small diagonal coefficient in an active cell
2444  if (this%active(ieq) > 0) then
2445  diagval = -done
2446  adiag = abs(this%system_matrix%get_diag_value(irow_glo))
2447  if (adiag < dem15) then
2448  call this%system_matrix%set_diag_value(irow_glo, diagval)
2449  this%rhs(ieq) = this%rhs(ieq) + diagval * this%x(ieq)
2450  end if
2451  ! -- Dirichlet boundary or no-flow cell
2452  else
2453  call this%system_matrix%set_diag_value(irow_glo, done)
2454  call this%system_matrix%zero_row_offdiag(irow_glo)
2455  this%rhs(ieq) = this%x(ieq)
2456  end if
2457  end do
2458  !
2459  ! -- complete adjustments for Dirichlet boundaries for a symmetric matrix
2460  if (this%isymmetric == 1 .and. simulation_mode == "SEQUENTIAL") then
2461  do ieq = 1, this%neq
2462  if (this%active(ieq) > 0) then
2463  icol_s = this%system_matrix%get_first_col_pos(ieq)
2464  icol_e = this%system_matrix%get_last_col_pos(ieq)
2465  do ipos = icol_s, icol_e
2466  jcol = this%system_matrix%get_column(ipos)
2467  if (jcol == ieq) cycle
2468  if (this%active(jcol) < 0) then
2469  this%rhs(ieq) = this%rhs(ieq) - &
2470  (this%system_matrix%get_value_pos(ipos) * &
2471  this%x(jcol))
2472  call this%system_matrix%set_value_pos(ipos, dzero)
2473  end if
2474 
2475  end do
2476  end if
2477  end do
2478  end if
2479  !
2480  ! -- pseudo transient continuation
2481  !
2482  ! -- set iallowptc
2483  ! -- no_ptc_option is FIRST
2484  if (this%iallowptc < 0) then
2485  if (kper > 1) then
2486  iallowptc = 1
2487  else
2488  iallowptc = 0
2489  end if
2490  !
2491  ! -- no_ptc_option is ALL (0) or using PTC (1)
2492  else
2493  iallowptc = this%iallowptc
2494  end if
2495  !
2496  ! -- set iptct
2497  iptct = iptc * iallowptc
2498  !
2499  ! -- calculate or modify pseudo transient continuation terms and add
2500  ! to amat diagonals
2501  if (iptct /= 0) then
2502  call this%sln_l2norm(l2norm)
2503  ! -- confirm that the l2norm exceeds previous l2norm
2504  ! if not, there is no need to add ptc terms
2505  if (kiter == 1) then
2506  if (kper > 1 .or. kstp > 1) then
2507  if (l2norm <= this%l2norm0) then
2508  iptc = 0
2509  end if
2510  end if
2511  else
2512  lsame = is_close(l2norm, this%l2norm0)
2513  if (lsame) then
2514  iptc = 0
2515  end if
2516  end if
2517  end if
2518  iptct = iptc * iallowptc
2519  if (iptct /= 0) then
2520  if (kiter == 1) then
2521  if (this%iptcout > 0) then
2522  write (this%iptcout, '(A10,6(1x,A15))') 'OUTER ITER', &
2523  ' PTCDEL', ' L2NORM0', ' L2NORM', &
2524  ' RHSNORM', ' 1/PTCDEL', ' RHSNORM/L2NORM'
2525  end if
2526  if (this%ptcdel0 > dzero) then
2527  this%ptcdel = this%ptcdel0
2528  else
2529  if (this%iptcopt == 0) then
2530  !
2531  ! -- ptcf is the reciprocal of the pseudo-time step
2532  this%ptcdel = done / ptcf
2533  else
2534  bnorm = dzero
2535  do ieq = 1, this%neq
2536  if (this%active(ieq) .gt. 0) then
2537  bnorm = bnorm + this%rhs(ieq) * this%rhs(ieq)
2538  end if
2539  end do
2540  bnorm = sqrt(bnorm)
2541  this%ptcdel = bnorm / l2norm
2542  end if
2543  end if
2544  else
2545  if (l2norm > dzero) then
2546  this%ptcdel = this%ptcdel * (this%l2norm0 / l2norm)**this%ptcexp
2547  else
2548  this%ptcdel = dzero
2549  end if
2550  end if
2551  if (this%ptcdel > dzero) then
2552  ptcval = done / this%ptcdel
2553  else
2554  ptcval = done
2555  end if
2556  bnorm = dzero
2557  do ieq = 1, this%neq
2558  irow_glo = ieq + this%matrix_offset
2559  if (this%active(ieq) > 0) then
2560  diagval = abs(this%system_matrix%get_diag_value(irow_glo))
2561  bnorm = bnorm + this%rhs(ieq) * this%rhs(ieq)
2562  call this%system_matrix%add_diag_value(irow_glo, -ptcval)
2563  this%rhs(ieq) = this%rhs(ieq) - ptcval * this%x(ieq)
2564  end if
2565  end do
2566  bnorm = sqrt(bnorm)
2567  if (this%iptcout > 0) then
2568  write (this%iptcout, '(i10,5(1x,e15.7),1(1x,f15.6))') &
2569  kiter, this%ptcdel, this%l2norm0, l2norm, bnorm, &
2570  ptcval, bnorm / l2norm
2571  end if
2572  this%l2norm0 = l2norm
2573  end if
2574  !
2575  ! -- save rhs, amat to a file
2576  ! to enable set itestmat to 1 and recompile
2577  !-------------------------------------------------------
2578  itestmat = 0
2579  if (itestmat == 1) then
2580  write (fname, fmtfname) this%id, kper, kstp, kiter
2581  print *, 'Saving amat to: ', trim(adjustl(fname))
2582 
2583  itestmat = getunit()
2584  open (itestmat, file=trim(adjustl(fname)))
2585  write (itestmat, *) 'NODE, RHS, AMAT FOLLOW'
2586  do ieq = 1, this%neq
2587  irow_glo = ieq + this%matrix_offset
2588  icol_s = this%system_matrix%get_first_col_pos(irow_glo)
2589  icol_e = this%system_matrix%get_last_col_pos(irow_glo)
2590  write (itestmat, '(*(G0,:,","))') &
2591  irow_glo, &
2592  this%rhs(ieq), &
2593  (this%system_matrix%get_column(ipos), ipos=icol_s, icol_e), &
2594  (this%system_matrix%get_value_pos(ipos), ipos=icol_s, icol_e)
2595  end do
2596  close (itestmat)
2597  !stop
2598  end if
2599  !-------------------------------------------------------
2600  !
2601  ! -- call appropriate linear solver
2602  !
2603  ! -- ims linear solver - linmeth option 1
2604  if (this%linsolver == ims_solver) then
2605  call this%imslinear%imslinear_apply(this%icnvg, kstp, kiter, in_iter, &
2606  this%nitermax, this%convnmod, &
2607  this%convmodstart, this%caccel, &
2608  this%cnvg_summary)
2609  else if (this%linsolver == petsc_solver) then
2610  call this%linear_solver%solve(kiter, this%vec_rhs, &
2611  this%vec_x, this%cnvg_summary)
2612  in_iter = this%linear_solver%iteration_number
2613  this%icnvg = this%linear_solver%is_converged
2614  end if
2615  end subroutine sln_ls
2616 
2617  !
2618  !> @ brief Set default Picard iteration variables
2619  !!
2620  !! Set default Picard iteration variables based on passed complexity option.
2621  !!
2622  !<
2623  subroutine sln_setouter(this, ifdparam)
2624  ! -- dummy variables
2625  class(numericalsolutiontype), intent(inout) :: this !< NumericalSolutionType instance
2626  integer(I4B), intent(in) :: ifdparam !< complexity option (1) simple (2) moderate (3) complex
2627  !
2628  ! -- simple option
2629  select case (ifdparam)
2630  case (1)
2631  this%dvclose = dem3
2632  this%mxiter = 25
2633  this%nonmeth = 0
2634  this%theta = done
2635  this%akappa = dzero
2636  this%gamma = done
2637  this%amomentum = dzero
2638  this%numtrack = 0
2639  this%btol = dzero
2640  this%breduc = dzero
2641  this%res_lim = dzero
2642  !
2643  ! -- moderate
2644  case (2)
2645  this%dvclose = dem2
2646  this%mxiter = 50
2647  this%nonmeth = 3
2648  this%theta = 0.9d0
2649  this%akappa = 0.0001d0
2650  this%gamma = dzero
2651  this%amomentum = dzero
2652  this%numtrack = 0
2653  this%btol = dzero
2654  this%breduc = dzero
2655  this%res_lim = dzero
2656  !
2657  ! -- complex
2658  case (3)
2659  this%dvclose = dem1
2660  this%mxiter = 100
2661  this%nonmeth = 3
2662  this%theta = 0.8d0
2663  this%akappa = 0.0001d0
2664  this%gamma = dzero
2665  this%amomentum = dzero
2666  this%numtrack = 20
2667  this%btol = 1.05d0
2668  this%breduc = 0.1d0
2669  this%res_lim = 0.002d0
2670  end select
2671  end subroutine sln_setouter
2672 
2673  !> @ brief Perform backtracking
2674  !!
2675  !! Perform backtracking on the solution in the maximum number of backtrack
2676  !! iterations (nbtrack) is greater than 0 and the backtracking criteria
2677  !! are exceeded.
2678  !!
2679  !<
2680  subroutine sln_backtracking(this, mp, cp, kiter)
2681  ! -- dummy variables
2682  class(numericalsolutiontype), intent(inout) :: this !< NumericalSolutionType instance
2683  class(numericalmodeltype), pointer :: mp !< model pointer (currently null())
2684  class(numericalexchangetype), pointer :: cp !< exchange pointer (currently null())
2685  integer(I4B), intent(in) :: kiter !< Picard iteration number
2686  ! -- local variables
2687  character(len=7) :: cmsg
2688  integer(I4B) :: nb
2689  integer(I4B) :: btflag
2690  integer(I4B) :: ibflag
2691  integer(I4B) :: ibtcnt
2692  real(DP) :: resin
2693  !
2694  ! -- initialize local variables
2695  ibflag = 0
2696 
2697  !
2698  ! -- refill amat and rhs with standard conductance
2699  call this%sln_buildsystem(kiter, inewton=0)
2700 
2701  !
2702  ! -- calculate initial l2 norm
2703  if (kiter == 1) then
2704  call this%sln_l2norm(this%res_prev)
2705  resin = this%res_prev
2706  ibflag = 0
2707  else
2708  call this%sln_l2norm(this%res_new)
2709  resin = this%res_new
2710  end if
2711  ibtcnt = 0
2712  if (kiter > 1) then
2713  if (this%res_new > this%res_prev * this%btol) then
2714  !
2715  ! -- iterate until backtracking complete
2716  btloop: do nb = 1, this%numtrack
2717  !
2718  ! -- backtrack the dependent variable
2719  call this%sln_backtracking_xupdate(btflag)
2720  !
2721  ! -- dependent-variable change less than dvclose
2722  if (btflag == 0) then
2723  ibflag = 4
2724  exit btloop
2725  end if
2726  !
2727  ibtcnt = nb
2728 
2729  ! recalculate linear system (amat and rhs)
2730  call this%sln_buildsystem(kiter, inewton=0)
2731 
2732  !
2733  ! -- calculate updated l2norm
2734  call this%sln_l2norm(this%res_new)
2735  !
2736  ! -- evaluate if back tracking can be terminated
2737  if (nb == this%numtrack) then
2738  ibflag = 2
2739  exit btloop
2740  end if
2741  if (this%res_new < this%res_prev * this%btol) then
2742  ibflag = 1
2743  exit btloop
2744  end if
2745  if (this%res_new < this%res_lim) then
2746  exit btloop
2747  end if
2748  end do btloop
2749  end if
2750  ! -- save new residual
2751  this%res_prev = this%res_new
2752  end if
2753  !
2754  ! -- write back backtracking results
2755  if (this%iprims > 0) then
2756  if (ibtcnt > 0) then
2757  cmsg = ' '
2758  else
2759  cmsg = '*'
2760  end if
2761  !
2762  ! -- add data to outertab
2763  call this%outertab%add_term('Backtracking')
2764  call this%outertab%add_term(kiter)
2765  call this%outertab%add_term(' ')
2766  if (this%numtrack > 0) then
2767  call this%outertab%add_term(ibflag)
2768  call this%outertab%add_term(ibtcnt)
2769  call this%outertab%add_term(resin)
2770  call this%outertab%add_term(this%res_prev)
2771  end if
2772  call this%outertab%add_term(' ')
2773  call this%outertab%add_term(cmsg)
2774  call this%outertab%add_term(' ')
2775  end if
2776  end subroutine sln_backtracking
2777 
2778  !> @ brief Backtracking update of the dependent variable
2779  !!
2780  !! Backtracking update of the dependent variable if the calculated backtracking
2781  !! update exceeds the dependent variable closure criteria.
2782  !!
2783  !<
2784  subroutine sln_backtracking_xupdate(this, bt_flag)
2785  ! -- dummy variables
2786  class(numericalsolutiontype), intent(inout) :: this !< NumericalSolutionType instance
2787  integer(I4B), intent(inout) :: bt_flag !< backtracking flag (1) backtracking performed (0) backtracking not performed
2788 
2789  bt_flag = this%get_backtracking_flag()
2790 
2791  ! perform backtracking if ...
2792  if (bt_flag > 0) then
2793  call this%apply_backtracking()
2794  end if
2795 
2796  end subroutine sln_backtracking_xupdate
2797 
2798  !> @brief Check if backtracking should be applied for this solution,
2799  !< returns 1: yes, 0: no
2800  function get_backtracking_flag(this) result(bt_flag)
2801  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2802  integer(I4B) :: bt_flag !< backtracking flag (1) backtracking performed (0) backtracking not performed
2803  ! local
2804  integer(I4B) :: n
2805  real(dp) :: dx
2806  real(dp) :: dx_abs
2807  real(dp) :: dx_abs_max
2808 
2809  ! default is off
2810  bt_flag = 0
2811 
2812  ! find max. change
2813  dx_abs_max = 0.0
2814  do n = 1, this%neq
2815  if (this%active(n) < 1) cycle
2816  dx = this%x(n) - this%xtemp(n)
2817  dx_abs = abs(dx)
2818  if (dx_abs > dx_abs_max) dx_abs_max = dx_abs
2819  end do
2820 
2821  ! if backtracking, set flag
2822  if (this%breduc * dx_abs_max >= this%dvclose) then
2823  bt_flag = 1
2824  end if
2825 
2826  end function get_backtracking_flag
2827 
2828  !> @brief Update x with backtracking
2829  !<
2830  subroutine apply_backtracking(this)
2831  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2832  ! local
2833  integer(I4B) :: n
2834  real(DP) :: delx
2835 
2836  do n = 1, this%neq
2837  if (this%active(n) < 1) cycle
2838  delx = this%breduc * (this%x(n) - this%xtemp(n))
2839  this%x(n) = this%xtemp(n) + delx
2840  end do
2841 
2842  end subroutine
2843 
2844  !> @ brief Calculate the solution L-2 norm for all
2845  !! active cells using
2846  !!
2847  !! A = the linear system matrix
2848  !! x = the dependent variable vector
2849  !! b = the right-hand side vector
2850  !!
2851  !! r = A * x - b
2852  !! r_i = 0 if cell i is inactive
2853  !! L2norm = || r ||_2
2854  !<
2855  subroutine sln_l2norm(this, l2norm)
2856  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2857  real(DP) :: l2norm !< calculated L-2 norm
2858  ! local
2859  class(vectorbasetype), pointer :: vec_resid
2860 
2861  ! calc. residual vector
2862  vec_resid => this%system_matrix%create_vec(this%neq)
2863  call this%sln_calc_residual(vec_resid)
2864 
2865  ! 2-norm
2866  l2norm = vec_resid%norm2()
2867 
2868  ! clean up temp. vector
2869  call vec_resid%destroy()
2870  deallocate (vec_resid)
2871  end subroutine sln_l2norm
2872 
2873  !> @ brief Get the maximum value from a vector
2874  !!
2875  !! Return the maximum value in a vector using a normalized form.
2876  !!
2877  !<
2878  subroutine sln_maxval(this, nsize, v, vmax)
2879  ! -- dummy variables
2880  class(numericalsolutiontype), intent(inout) :: this !< NumericalSolutionType instance
2881  integer(I4B), intent(in) :: nsize !< length of vector
2882  real(DP), dimension(nsize), intent(in) :: v !< input vector
2883  real(DP), intent(inout) :: vmax !< maximum value
2884  ! -- local variables
2885  integer(I4B) :: n
2886  real(DP) :: d
2887  real(DP) :: denom
2888  real(DP) :: dnorm
2889  !
2890  ! -- determine maximum value
2891  vmax = v(1)
2892  do n = 2, nsize
2893  d = v(n)
2894  denom = abs(vmax)
2895  if (denom == dzero) then
2896  denom = dprec
2897  end if
2898  !
2899  ! -- calculate normalized value
2900  dnorm = abs(d) / denom
2901  if (dnorm > done) then
2902  vmax = d
2903  end if
2904  end do
2905  end subroutine sln_maxval
2906 
2907  !> @ brief Calculate dependent-variable change
2908  !!
2909  !! Calculate the dependent-variable change for every cell.
2910  !!
2911  !<
2912  subroutine sln_calcdx(this, neq, active, x, xtemp, dx)
2913  ! -- dummy variables
2914  class(numericalsolutiontype), intent(inout) :: this !< NumericalSolutionType instance
2915  integer(I4B), intent(in) :: neq !< number of equations
2916  integer(I4B), dimension(neq), intent(in) :: active !< active cell flag (1)
2917  real(DP), dimension(neq), intent(in) :: x !< current dependent-variable
2918  real(DP), dimension(neq), intent(in) :: xtemp !< previous dependent-variable
2919  real(DP), dimension(neq), intent(inout) :: dx !< dependent-variable change
2920  ! -- local
2921  integer(I4B) :: n
2922  !
2923  ! -- calculate dependent-variable change
2924  do n = 1, neq
2925  ! -- skip inactive nodes
2926  if (active(n) < 1) then
2927  dx(n) = dzero
2928  else
2929  dx(n) = x(n) - xtemp(n)
2930  end if
2931  end do
2932  end subroutine sln_calcdx
2933 
2934  !> @brief Calculate pseudo-transient continuation factor
2935  !< from the models in the solution
2936  subroutine sln_calc_ptc(this, iptc, ptcf)
2937  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2938  integer(I4B) :: iptc !< PTC (1) or not (0)
2939  real(DP) :: ptcf !< the PTC factor calculated
2940  ! local
2941  integer(I4B) :: im
2942  class(numericalmodeltype), pointer :: mp
2943  class(vectorbasetype), pointer :: vec_resid
2944 
2945  iptc = 0
2946  ptcf = dzero
2947 
2948  ! calc. residual vector
2949  vec_resid => this%system_matrix%create_vec(this%neq)
2950  call this%sln_calc_residual(vec_resid)
2951 
2952  ! determine ptc
2953  do im = 1, this%modellist%Count()
2954  mp => getnumericalmodelfromlist(this%modellist, im)
2955  call mp%model_ptc(vec_resid, iptc, ptcf)
2956  end do
2957 
2958  ! clean up temp. vector
2959  call vec_resid%destroy()
2960  deallocate (vec_resid)
2961 
2962  end subroutine sln_calc_ptc
2963 
2964  !> @brief Calculate the current residual vector r = A*x - b,
2965  !< zeroes out for inactive cells
2966  subroutine sln_calc_residual(this, vec_resid)
2967  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2968  class(vectorbasetype), pointer :: vec_resid !< the residual vector
2969  ! local
2970  integer(I4B) :: n
2971 
2972  call this%system_matrix%multiply(this%vec_x, vec_resid) ! r = A*x
2973 
2974  call vec_resid%axpy(-1.0_dp, this%vec_rhs) ! r = r - b
2975 
2976  do n = 1, this%neq
2977  if (this%active(n) < 1) then
2978  call vec_resid%set_value_local(n, 0.0_dp) ! r_i = 0 if inactive
2979  end if
2980  end do
2981 
2982  end subroutine sln_calc_residual
2983 
2984  !> @ brief Under-relaxation
2985  !!
2986  !! Under relax using the simple, cooley, or delta-bar-delta methods.
2987  !!
2988  !<
2989  subroutine sln_underrelax(this, kiter, bigch, neq, active, x, xtemp)
2990  ! -- dummy variables
2991  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
2992  integer(I4B), intent(in) :: kiter !< Picard iteration number
2993  real(DP), intent(in) :: bigch !< maximum dependent-variable change
2994  integer(I4B), intent(in) :: neq !< number of equations
2995  integer(I4B), dimension(neq), intent(in) :: active !< active cell flag (1)
2996  real(DP), dimension(neq), intent(inout) :: x !< current dependent-variable
2997  real(DP), dimension(neq), intent(in) :: xtemp !< previous dependent-variable
2998  ! -- local variables
2999  integer(I4B) :: n
3000  real(DP) :: ww
3001  real(DP) :: delx
3002  real(DP) :: relax
3003  real(DP) :: es
3004  real(DP) :: aes
3005  real(DP) :: amom
3006  !
3007  ! -- option for using simple dampening (as done by MODFLOW-2005 PCG)
3008  if (this%nonmeth == 1) then
3009  do n = 1, neq
3010  !
3011  ! -- skip inactive nodes
3012  if (active(n) < 1) cycle
3013  !
3014  ! -- compute step-size (delta x)
3015  delx = x(n) - xtemp(n)
3016  this%dxold(n) = delx
3017 
3018  ! -- dampen dependent variable solution
3019  x(n) = xtemp(n) + this%gamma * delx
3020  end do
3021  !
3022  ! -- option for using cooley underrelaxation
3023  else if (this%nonmeth == 2) then
3024  !
3025  ! -- set bigch
3026  this%bigch = bigch
3027  !
3028  ! -- initialize values for first iteration
3029  if (kiter == 1) then
3030  relax = done
3031  this%relaxold = done
3032  this%bigchold = bigch
3033  else
3034  !
3035  ! -- compute relaxation factor
3036  es = this%bigch / (this%bigchold * this%relaxold)
3037  aes = abs(es)
3038  if (es < -done) then
3039  relax = dhalf / aes
3040  else
3041  relax = (dthree + es) / (dthree + aes)
3042  end if
3043  end if
3044  this%relaxold = relax
3045  !
3046  ! -- modify cooley to use weighted average of past changes
3047  this%bigchold = (done - this%gamma) * this%bigch + this%gamma * &
3048  this%bigchold
3049  !
3050  ! -- compute new dependent variable after under-relaxation
3051  if (relax < done) then
3052  do n = 1, neq
3053  !
3054  ! -- skip inactive nodes
3055  if (active(n) < 1) cycle
3056  !
3057  ! -- update dependent variable
3058  delx = x(n) - xtemp(n)
3059  this%dxold(n) = delx
3060  x(n) = xtemp(n) + relax * delx
3061  end do
3062  end if
3063  !
3064  ! -- option for using delta-bar-delta scheme to under-relax for all equations
3065  else if (this%nonmeth == 3) then
3066  do n = 1, neq
3067  !
3068  ! -- skip inactive nodes
3069  if (active(n) < 1) cycle
3070  !
3071  ! -- compute step-size (delta x) and initialize d-b-d parameters
3072  delx = x(n) - xtemp(n)
3073  !
3074  ! -- initialize values for first iteration
3075  if (kiter == 1) then
3076  this%wsave(n) = done
3077  this%hchold(n) = dem20
3078  this%deold(n) = dzero
3079  end if
3080  !
3081  ! -- compute new relaxation term as per delta-bar-delta
3082  ww = this%wsave(n)
3083  !
3084  ! for flip-flop condition, decrease factor
3085  if (this%deold(n) * delx < dzero) then
3086  ww = this%theta * this%wsave(n)
3087  ! -- when change is of same sign, increase factor
3088  else
3089  ww = this%wsave(n) + this%akappa
3090  end if
3091  if (ww > done) ww = done
3092  this%wsave(n) = ww
3093  !
3094  ! -- compute weighted average of past changes in hchold
3095  if (kiter == 1) then
3096  this%hchold(n) = delx
3097  else
3098  this%hchold(n) = (done - this%gamma) * delx + &
3099  this%gamma * this%hchold(n)
3100  end if
3101  !
3102  ! -- store slope (change) term for next iteration
3103  this%deold(n) = delx
3104  this%dxold(n) = delx
3105  !
3106  ! -- compute accepted step-size and new dependent variable
3107  amom = dzero
3108  if (kiter > 4) amom = this%amomentum
3109  delx = delx * ww + amom * this%hchold(n)
3110  x(n) = xtemp(n) + delx
3111  end do
3112  !
3113  end if
3114  end subroutine sln_underrelax
3115 
3116  !> @ brief Determine maximum dependent-variable change
3117  !!
3118  !! Determine the maximum dependent-variable change at the end of a
3119  !! Picard iteration.
3120  !!
3121  !<
3122  subroutine sln_get_dxmax(this, hncg, lrch)
3123  ! -- dummy variables
3124  class(numericalsolutiontype), intent(inout) :: this !< NumericalSolutionType instance
3125  real(DP), intent(inout) :: hncg !< maximum dependent-variable change
3126  integer(I4B), intent(inout) :: lrch !< location of the maximum dependent-variable change
3127  ! -- local variables
3128  integer(I4B) :: nb
3129  real(DP) :: bigch
3130  real(DP) :: abigch
3131  integer(I4B) :: n
3132  real(DP) :: hdif
3133  real(DP) :: ahdif
3134  !
3135  ! -- determine the maximum change
3136  nb = 0
3137  bigch = dzero
3138  abigch = dzero
3139  do n = 1, this%neq
3140  if (this%active(n) < 1) cycle
3141  hdif = this%x(n) - this%xtemp(n)
3142  ahdif = abs(hdif)
3143  if (ahdif > abigch) then
3144  bigch = hdif
3145  abigch = ahdif
3146  nb = n
3147  end if
3148  end do
3149  !
3150  !-----store maximum change value and location
3151  hncg = bigch
3152  lrch = nb
3153  end subroutine sln_get_dxmax
3154 
3155  function sln_has_converged(this, max_dvc) result(has_converged)
3156  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
3157  real(dp) :: max_dvc !< the maximum dependent variable change
3158  logical(LGP) :: has_converged !< True, when converged
3159 
3160  has_converged = .false.
3161  if (abs(max_dvc) <= this%dvclose) then
3162  has_converged = .true.
3163  end if
3164 
3165  end function sln_has_converged
3166 
3167  !> @brief Check package convergence
3168  !<
3169  function sln_package_convergence(this, dpak, cpakout, iend) result(ivalue)
3170  ! dummy
3171  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
3172  real(dp), intent(in) :: dpak !< Newton Under-relaxation flag
3173  character(len=LENPAKLOC), intent(in) :: cpakout !< string with package that caused failure
3174  integer(I4B), intent(in) :: iend !< flag indicating if last inner iteration (iend=1)
3175  ! local
3176  integer(I4B) :: ivalue
3177  ivalue = 1
3178  if (abs(dpak) > this%dvclose) then
3179  ivalue = 0
3180  ! -- write message to stdout
3181  if (iend /= 0) then
3182  write (errmsg, '(3a)') &
3183  'PACKAGE (', trim(cpakout), ') CAUSED CONVERGENCE FAILURE'
3184  call write_message(errmsg)
3185  end if
3186  end if
3187 
3188  end function sln_package_convergence
3189 
3190  !> @brief Synchronize Newton Under-relaxation flag
3191  !<
3192  function sln_sync_newtonur_flag(this, inewtonur) result(ivalue)
3193  ! dummy
3194  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
3195  integer(I4B), intent(in) :: inewtonur !< Newton Under-relaxation flag
3196  ! local
3197  integer(I4B) :: ivalue !< Default is set to current value (1 = under-relaxation applied)
3198 
3199  ivalue = inewtonur
3200 
3201  end function sln_sync_newtonur_flag
3202 
3203  !> @brief Custom convergence check for when Newton UR has been applied
3204  !<
3205  function sln_nur_has_converged(this, dxold_max, hncg) &
3206  result(has_converged)
3207  class(numericalsolutiontype) :: this !< NumericalSolutionType instance
3208  real(dp), intent(in) :: dxold_max !< the maximum dependent variable change for unrelaxed cells
3209  real(dp), intent(in) :: hncg !< largest dep. var. change at end of Picard iteration
3210  logical(LGP) :: has_converged !< True, when converged
3211 
3212  has_converged = .false.
3213  if (abs(dxold_max) <= this%dvclose .and. &
3214  abs(hncg) <= this%dvclose) then
3215  has_converged = .true.
3216  end if
3217 
3218  end function sln_nur_has_converged
3219 
3220  !> @ brief Get cell location string
3221  !!
3222  !! Get the cell location string for the provided solution node number.
3223  !< Returns empty string when node not found.
3224  subroutine sln_get_loc(this, nodesln, str)
3225  ! -- dummy variables
3226  class(numericalsolutiontype), intent(inout) :: this !< NumericalSolutionType instance
3227  integer(I4B), intent(in) :: nodesln !< solution node number
3228  character(len=*), intent(inout) :: str !< string with user node number
3229  ! -- local variables
3230  class(numericalmodeltype), pointer :: mp => null()
3231  integer(I4B) :: i
3232  integer(I4B) :: istart
3233  integer(I4B) :: iend
3234  integer(I4B) :: noder
3235  integer(I4B) :: nglo
3236  !
3237  ! -- initialize dummy variables
3238  str = ''
3239  !
3240  ! -- initialize local variables
3241  noder = 0
3242  !
3243  ! -- when parallel, account for offset
3244  nglo = nodesln + this%matrix_offset
3245  !
3246  ! -- calculate and set offsets
3247  do i = 1, this%modellist%Count()
3248  mp => getnumericalmodelfromlist(this%modellist, i)
3249  istart = 0
3250  iend = 0
3251  call mp%get_mrange(istart, iend)
3252  if (nglo >= istart .and. nglo <= iend) then
3253  noder = nglo - istart + 1
3254  call mp%get_mcellid(noder, str)
3255  exit
3256  end if
3257  end do
3258  end subroutine sln_get_loc
3259 
3260  !> @ brief Get user node number
3261  !!
3262  !! Get the user node number from a model for the provided solution node number.
3263  !!
3264  !<
3265  subroutine sln_get_nodeu(this, nodesln, im, nodeu)
3266  ! -- dummy variables
3267  class(numericalsolutiontype), intent(inout) :: this !< NumericalSolutionType instance
3268  integer(I4B), intent(in) :: nodesln !< solution node number
3269  integer(I4B), intent(inout) :: im !< solution model index (index in model list for this solution)
3270  integer(I4B), intent(inout) :: nodeu !< user node number
3271  ! -- local variables
3272  class(numericalmodeltype), pointer :: mp => null()
3273  integer(I4B) :: i
3274  integer(I4B) :: istart
3275  integer(I4B) :: iend
3276  integer(I4B) :: noder, nglo
3277  !
3278  ! -- initialize local variables
3279  noder = 0
3280  !
3281  ! -- when parallel, account for offset
3282  nglo = nodesln + this%matrix_offset
3283  !
3284  ! -- calculate and set offsets
3285  do i = 1, this%modellist%Count()
3286  mp => getnumericalmodelfromlist(this%modellist, i)
3287  istart = 0
3288  iend = 0
3289  call mp%get_mrange(istart, iend)
3290  if (nglo >= istart .and. nglo <= iend) then
3291  noder = nglo - istart + 1
3292  call mp%get_mnodeu(noder, nodeu)
3293  im = i
3294  exit
3295  end if
3296  end do
3297  end subroutine sln_get_nodeu
3298 
3299  !> @ brief Cast a object as a Numerical Solution
3300  !!
3301  !! Get a numerical solution from a list.
3302  !!
3303  !<
3304  function castasnumericalsolutionclass(obj) result(res)
3305  ! -- dummy variables
3306  class(*), pointer, intent(inout) :: obj !< generic object
3307  ! -- return variable
3308  class(numericalsolutiontype), pointer :: res !< output NumericalSolutionType
3309  !
3310  ! -- initialize return variable
3311  res => null()
3312  !
3313  ! -- determine if obj is associated
3314  if (.not. associated(obj)) return
3315  !
3316  ! -- set res
3317  select type (obj)
3318  class is (numericalsolutiontype)
3319  res => obj
3320  end select
3321  end function castasnumericalsolutionclass
3322 
3323  !> @ brief Get a numerical solution
3324  !!
3325  !! Get a numerical solution from a list.
3326  !!
3327  !<
3328  function getnumericalsolutionfromlist(list, idx) result(res)
3329  ! -- dummy variables
3330  type(listtype), intent(inout) :: list !< list of numerical solutions
3331  integer(I4B), intent(in) :: idx !< value to retrieve from the list
3332  ! -- return variables
3333  class(numericalsolutiontype), pointer :: res !< numerical solution
3334  ! -- local variables
3335  class(*), pointer :: obj
3336  !
3337  obj => list%GetItem(idx)
3338  res => castasnumericalsolutionclass(obj)
3339  end function getnumericalsolutionfromlist
3340 end module numericalsolutionmodule
subroutine, public ats_submit_delt(kstp, kper, dt, sloc, idir)
@ brief Allow and external caller to submit preferred time step
Definition: ats.f90:493
subroutine, public addbasesolutiontolist(list, solution)
This module contains block parser methods.
Definition: BlockParser.f90:7
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
@ tabright
right justified table column
Definition: Constants.f90:173
@ tableft
left justified table column
Definition: Constants.f90:171
@ mvalidate
validation mode - do not run time steps
Definition: Constants.f90:205
@ mnormal
normal output mode
Definition: Constants.f90:206
real(dp), parameter dem20
real constant 1e-20
Definition: Constants.f90:117
real(dp), parameter dep3
real constant 1000
Definition: Constants.f90:88
integer(i4b), parameter lensolutionname
maximum length of the solution name
Definition: Constants.f90:21
real(dp), parameter dep6
real constant 1000000
Definition: Constants.f90:89
real(dp), parameter donethird
real constant 1/3
Definition: Constants.f90:67
real(dp), parameter dnodata
real no data constant
Definition: Constants.f90:95
integer(i4b), parameter lenpakloc
maximum length of a package location
Definition: Constants.f90:50
real(dp), parameter dep20
real constant 1e20
Definition: Constants.f90:91
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
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 dprec
real constant machine precision
Definition: Constants.f90:120
real(dp), parameter dem15
real constant 1e-15
Definition: Constants.f90:116
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 dthree
real constant 3
Definition: Constants.f90:80
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
Definition: ErrorUtil.f90:24
subroutine allocate_scalars(this)
@ brief Allocate and initialize scalars
Definition: ImsLinear.f90:456
integer(i4b), parameter, public cg_method
integer(i4b) function, public getunit()
Get a free unit number.
subroutine, public append_processor_id(name, proc_id)
Append processor id to a string.
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
This module defines variable data types.
Definition: kind.f90:8
class(linearsolverbasetype) function, pointer, public create_linear_solver(solver_mode, sln_name)
Factory method to create the linear solver object.
type(listtype), public basesolutionlist
Definition: mf6lists.f90:19
pure logical function, public is_close(a, b, rtol, atol, symmetric)
Check if a real value is approximately equal to another.
Definition: MathUtil.f90:46
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
Store and issue logging messages to output units.
Definition: Message.f90:2
subroutine, public write_message(text, iunit, fmt, skipbefore, skipafter, advance)
Write a message to an output unit.
Definition: Message.f90:210
class(numericalexchangetype) function, pointer, public getnumericalexchangefromlist(list, idx)
Retrieve a specific numerical exchange from a list.
subroutine, public addnumericalexchangetolist(list, exchange)
Add numerical exchange to a list.
subroutine, public addnumericalmodeltolist(list, model)
class(numericalmodeltype) function, pointer, public getnumericalmodelfromlist(list, idx)
subroutine convergence_summary(this, iu, im, itertot_timestep)
@ brief Solution convergence summary
subroutine sln_get_nodeu(this, nodesln, im, nodeu)
@ brief Get user node number
integer(i4b) function sln_sync_newtonur_flag(this, inewtonur)
Synchronize Newton Under-relaxation flag.
subroutine save(this, filename)
@ brief Save solution data to a file
subroutine sln_backtracking_xupdate(this, bt_flag)
@ brief Backtracking update of the dependent variable
logical(lgp) function sln_nur_has_converged(this, dxold_max, hncg)
Custom convergence check for when Newton UR has been applied.
logical(lgp) function sln_has_converged(this, max_dvc)
integer(i4b), parameter petsc_solver
integer(i4b) function sln_package_convergence(this, dpak, cpakout, iend)
Check package convergence.
type(listtype) function, pointer get_exchanges(this)
Returns a pointer to the list of exchanges in this solution.
subroutine sln_l2norm(this, l2norm)
@ brief Calculate the solution L-2 norm for all active cells using
subroutine sln_connect(this)
@ brief Assign solution connections
subroutine sln_get_loc(this, nodesln, str)
@ brief Get cell location string
subroutine apply_backtracking(this)
Update x with backtracking.
subroutine writecsvheader(this)
@ brief CSV header
subroutine sln_buildsystem(this, kiter, inewton)
subroutine sln_maxval(this, nsize, v, vmax)
@ brief Get the maximum value from a vector
subroutine sln_calc_residual(this, vec_resid)
Calculate the current residual vector r = A*x - b,.
subroutine sln_backtracking(this, mp, cp, kiter)
@ brief Perform backtracking
subroutine sln_calc_ptc(this, iptc, ptcf)
Calculate pseudo-transient continuation factor.
subroutine sln_underrelax(this, kiter, bigch, neq, active, x, xtemp)
@ brief Under-relaxation
class(numericalsolutiontype) function, pointer, public castasnumericalsolutionclass(obj)
@ brief Cast a object as a Numerical Solution
type(listtype) function, pointer get_models(this)
Get a list of models.
subroutine finalizesolve(this, kiter, isgcnvg, isuppress_output)
@ brief finalize a solution
subroutine sln_setouter(this, ifdparam)
@ brief Set default Picard iteration variables
subroutine sln_ls(this, kiter, kstp, kper, in_iter, iptc, ptcf)
@ brief Solve the linear system of equations
subroutine sln_calcdx(this, neq, active, x, xtemp, dx)
@ brief Calculate dependent-variable change
integer(i4b), parameter ims_solver
subroutine csv_convergence_summary(this, iu, totim, kper, kstp, kouter, niter, istart, kstart)
@ brief Solution convergence CSV summary
subroutine sln_reset(this)
@ brief Reset the solution
class(numericalsolutiontype) function, pointer, public getnumericalsolutionfromlist(list, idx)
@ brief Get a numerical solution
subroutine allocate_arrays(this)
@ brief Allocate arrays
integer(i4b) function get_backtracking_flag(this)
Check if backtracking should be applied for this solution,.
subroutine preparesolve(this)
@ brief prepare to solve
subroutine writeptcinfotofile(this, kper)
@ brief PTC header
subroutine add_exchange(this, exchange)
Add exchange.
subroutine add_model(this, mp)
@ brief Add a model
subroutine, public create_numerical_solution(num_sol, filename, id)
@ brief Create a new solution
subroutine sln_get_dxmax(this, hncg, lrch)
@ brief Determine maximum dependent-variable change
type(profilertype), public g_prof
the global timer object (to reduce trivial lines of code)
Definition: Profiler.f90:65
subroutine print(this, output_unit)
Definition: Profiler.f90:191
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
integer(i4b), parameter, public stg_bfr_exg_ad
before exchange advance (per solution)
Definition: SimStages.f90:21
integer(i4b), parameter, public stg_bfr_exg_cf
before exchange calculate (per solution)
Definition: SimStages.f90:22
integer(i4b), parameter, public stg_bfr_exg_fc
before exchange formulate (per solution)
Definition: SimStages.f90:23
integer(i4b), parameter, public stg_bfr_exg_ac
before exchange add connections (per solution)
Definition: SimStages.f90:16
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
character(len=linelength) simulation_mode
integer(i4b) nr_procs
integer(i4b) iout
file unit number for simulation output
integer(i4b) proc_id
integer(i4b) isim_mode
simulation mode
subroutine, public table_cr(this, name, title)
Definition: Table.f90:87
real(dp), pointer, public totim
time relative to start of simulation
Definition: tdis.f90:32
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:29
subroutine, public code_timer(it, t1, ts)
Get end time and calculate elapsed time.
Definition: Timer.f90:159
This module contains version information.
Definition: version.f90:7
integer(i4b), parameter idevelopmode
Definition: version.f90:19
Highest level model type. All models extend this parent type.
Definition: BaseModel.f90:13
This structure stores the generic convergence info for a solution.
A generic heterogeneous doubly-linked list.
Definition: List.f90:14