57 character(len=LENMEMPATH) :: memory_path
58 character(len=LINELENGTH) :: fname
59 character(len=16) :: solver_mode
62 integer(I4B),
pointer :: id
63 integer(I4B),
pointer :: iu
64 real(dp),
pointer :: ttform
65 real(dp),
pointer :: ttsoln
66 integer(I4B),
pointer :: isymmetric => null()
67 integer(I4B),
pointer :: neq => null()
68 integer(I4B),
pointer :: matrix_offset => null()
73 real(dp),
dimension(:),
pointer,
contiguous :: rhs => null()
74 real(dp),
dimension(:),
pointer,
contiguous :: x => null()
75 integer(I4B),
dimension(:),
pointer,
contiguous :: active => null()
76 real(dp),
dimension(:),
pointer,
contiguous :: xtemp => null()
80 real(dp),
pointer :: theta => null()
81 real(dp),
pointer :: akappa => null()
82 real(dp),
pointer :: gamma => null()
83 real(dp),
pointer :: amomentum => null()
84 real(dp),
pointer :: breduc => null()
85 real(dp),
pointer :: btol => null()
86 real(dp),
pointer :: res_lim => null()
87 real(dp),
pointer :: dvclose => null()
88 real(dp),
pointer :: bigchold => null()
89 real(dp),
pointer :: bigch => null()
90 real(dp),
pointer :: relaxold => null()
91 real(dp),
pointer :: res_prev => null()
92 real(dp),
pointer :: res_new => null()
93 integer(I4B),
pointer :: icnvg => null()
94 integer(I4B),
pointer :: itertot_timestep => null()
95 integer(I4B),
pointer :: iouttot_timestep => null()
96 integer(I4B),
pointer :: itertot_sim => null()
97 integer(I4B),
pointer :: mxiter => null()
98 integer(I4B),
pointer :: linsolver => null()
99 integer(I4B),
pointer :: nonmeth => null()
100 integer(I4B),
pointer :: numtrack => null()
101 integer(I4B),
pointer :: iprims => null()
102 integer(I4B),
pointer :: ibflag => null()
103 integer(I4B),
dimension(:, :),
pointer,
contiguous :: lrch => null()
104 real(dp),
dimension(:),
pointer,
contiguous :: hncg => null()
105 real(dp),
dimension(:),
pointer,
contiguous :: dxold => null()
106 real(dp),
dimension(:),
pointer,
contiguous :: deold => null()
107 real(dp),
dimension(:),
pointer,
contiguous :: wsave => null()
108 real(dp),
dimension(:),
pointer,
contiguous :: hchold => null()
111 character(len=31),
dimension(:),
pointer,
contiguous :: caccel => null()
112 integer(I4B),
pointer :: icsvouterout => null()
113 integer(I4B),
pointer :: icsvinnerout => null()
114 integer(I4B),
pointer :: nitermax => null()
115 integer(I4B),
pointer :: convnmod => null()
116 integer(I4B),
dimension(:),
pointer,
contiguous :: convmodstart => null()
123 integer(I4B),
pointer :: iallowptc => null()
124 integer(I4B),
pointer :: iptcopt => null()
125 integer(I4B),
pointer :: iptcout => null()
126 real(dp),
pointer :: l2norm0 => null()
127 real(dp),
pointer :: ptcdel => null()
128 real(dp),
pointer :: ptcdel0 => null()
129 real(dp),
pointer :: ptcexp => null()
132 real(dp),
pointer :: atsfrac => null()
145 class(*),
pointer :: synchronize_ctx => null()
207 integer(I4B) :: stage
208 class(*),
pointer :: ctx
227 character(len=*),
intent(in) :: filename
228 integer(I4B),
intent(in) :: id
230 integer(I4B) :: inunit
232 character(len=LENSOLUTIONNAME) :: solutionname
236 write (solutionname,
'(a, i0)')
'SLN_', id
238 num_sol%name = solutionname
240 allocate (num_sol%modellist)
241 allocate (num_sol%exchangelist)
243 call num_sol%allocate_scalars()
252 inquire (file=filename, number=inunit)
254 if (inunit < 0) inunit =
getunit()
256 write (
iout,
'(/a,a)')
' Creating solution: ', num_sol%name
260 call num_sol%parser%Initialize(num_sol%iu,
iout)
280 call mem_allocate(this%ttform,
'TTFORM', this%memory_path)
281 call mem_allocate(this%ttsoln,
'TTSOLN', this%memory_path)
282 call mem_allocate(this%isymmetric,
'ISYMMETRIC', this%memory_path)
284 call mem_allocate(this%matrix_offset,
'MATRIX_OFFSET', this%memory_path)
285 call mem_allocate(this%dvclose,
'DVCLOSE', this%memory_path)
286 call mem_allocate(this%bigchold,
'BIGCHOLD', this%memory_path)
287 call mem_allocate(this%bigch,
'BIGCH', this%memory_path)
288 call mem_allocate(this%relaxold,
'RELAXOLD', this%memory_path)
289 call mem_allocate(this%res_prev,
'RES_PREV', this%memory_path)
290 call mem_allocate(this%res_new,
'RES_NEW', this%memory_path)
291 call mem_allocate(this%icnvg,
'ICNVG', this%memory_path)
292 call mem_allocate(this%itertot_timestep,
'ITERTOT_TIMESTEP', this%memory_path)
293 call mem_allocate(this%iouttot_timestep,
'IOUTTOT_TIMESTEP', this%memory_path)
294 call mem_allocate(this%itertot_sim,
'INNERTOT_SIM', this%memory_path)
295 call mem_allocate(this%mxiter,
'MXITER', this%memory_path)
296 call mem_allocate(this%linsolver,
'LINSOLVER', this%memory_path)
297 call mem_allocate(this%nonmeth,
'NONMETH', this%memory_path)
298 call mem_allocate(this%iprims,
'IPRIMS', this%memory_path)
299 call mem_allocate(this%theta,
'THETA', this%memory_path)
300 call mem_allocate(this%akappa,
'AKAPPA', this%memory_path)
301 call mem_allocate(this%gamma,
'GAMMA', this%memory_path)
302 call mem_allocate(this%amomentum,
'AMOMENTUM', this%memory_path)
303 call mem_allocate(this%breduc,
'BREDUC', this%memory_path)
305 call mem_allocate(this%res_lim,
'RES_LIM', this%memory_path)
306 call mem_allocate(this%numtrack,
'NUMTRACK', this%memory_path)
307 call mem_allocate(this%ibflag,
'IBFLAG', this%memory_path)
308 call mem_allocate(this%icsvouterout,
'ICSVOUTEROUT', this%memory_path)
309 call mem_allocate(this%icsvinnerout,
'ICSVINNEROUT', this%memory_path)
310 call mem_allocate(this%nitermax,
'NITERMAX', this%memory_path)
311 call mem_allocate(this%convnmod,
'CONVNMOD', this%memory_path)
312 call mem_allocate(this%iallowptc,
'IALLOWPTC', this%memory_path)
313 call mem_allocate(this%iptcopt,
'IPTCOPT', this%memory_path)
314 call mem_allocate(this%iptcout,
'IPTCOUT', this%memory_path)
315 call mem_allocate(this%l2norm0,
'L2NORM0', this%memory_path)
316 call mem_allocate(this%ptcdel,
'PTCDEL', this%memory_path)
317 call mem_allocate(this%ptcdel0,
'PTCDEL0', this%memory_path)
318 call mem_allocate(this%ptcexp,
'PTCEXP', this%memory_path)
319 call mem_allocate(this%atsfrac,
'ATSFRAC', this%memory_path)
329 this%bigchold =
dzero
331 this%relaxold =
dzero
332 this%res_prev =
dzero
334 this%itertot_timestep = 0
335 this%iouttot_timestep = 0
344 this%amomentum =
dzero
350 this%icsvouterout = 0
351 this%icsvinnerout = 0
383 this%convnmod = this%modellist%Count()
386 call mem_allocate(this%active, this%neq,
'IACTIVE', this%memory_path)
387 call mem_allocate(this%xtemp, this%neq,
'XTEMP', this%memory_path)
388 call mem_allocate(this%dxold, this%neq,
'DXOLD', this%memory_path)
389 call mem_allocate(this%hncg, 0,
'HNCG', this%memory_path)
390 call mem_allocate(this%lrch, 3, 0,
'LRCH', this%memory_path)
391 call mem_allocate(this%wsave, 0,
'WSAVE', this%memory_path)
392 call mem_allocate(this%hchold, 0,
'HCHOLD', this%memory_path)
393 call mem_allocate(this%deold, 0,
'DEOLD', this%memory_path)
394 call mem_allocate(this%convmodstart, this%convnmod + 1,
'CONVMODSTART', &
399 this%xtemp(i) =
dzero
400 this%dxold(i) =
dzero
406 this%convmodstart(1) = ieq
407 do i = 1, this%modellist%Count()
410 this%convmodstart(i + 1) = ieq
435 integer(I4B),
allocatable,
dimension(:) :: rowmaxnnz
436 integer(I4B) :: ncol, irow_start, irow_end
437 integer(I4B) :: mod_offset
440 do i = 1, this%modellist%Count()
442 call mp%set_idsoln(this%id)
443 this%neq = this%neq + mp%neq
448 this%solver_mode =
'PETSC'
450 this%solver_mode =
'IMS'
454 allocate (this%linear_settings)
458 this%system_matrix => this%linear_solver%create_matrix()
459 this%vec_x => this%system_matrix%create_vec_mm(this%neq,
'X', &
461 this%x => this%vec_x%get_array()
462 this%vec_rhs => this%system_matrix%create_vec_mm(this%neq,
'RHS', &
464 this%rhs => this%vec_rhs%get_array()
466 call this%vec_rhs%get_ownership_range(irow_start, irow_end)
467 ncol = this%vec_rhs%get_size()
471 this%matrix_offset = irow_start - 1
472 do i = 1, this%modellist%Count()
479 call this%allocate_arrays()
482 allocate (this%cnvg_summary)
483 call this%cnvg_summary%init(this%modellist%Count(), this%convmodstart, &
487 do i = 1, this%modellist%Count()
489 call mp%set_xptr(this%x, this%matrix_offset,
'X', this%name)
490 call mp%set_rhsptr(this%rhs, this%matrix_offset,
'RHS', this%name)
491 call mp%set_iboundptr(this%active, this%matrix_offset,
'IBOUND', this%name)
495 allocate (rowmaxnnz(this%neq))
499 call this%sparse%init(this%neq, ncol, rowmaxnnz)
500 this%sparse%offset = this%matrix_offset
501 deallocate (rowmaxnnz)
504 call this%sln_connect()
526 character(len=linelength) :: warnmsg
527 character(len=linelength) :: keyword
528 character(len=linelength) :: fname
529 character(len=linelength) :: msg
531 integer(I4B) :: ifdparam, mxvl, npp
533 logical(LGP) :: isfound, endOfBlock
536 character(len=*),
parameter :: fmtcsvout = &
537 "(4x, 'CSV OUTPUT WILL BE SAVED TO FILE: ', a, &
538 &/4x, 'OPENED ON UNIT: ', I7)"
539 character(len=*),
parameter :: fmtptcout = &
540 "(4x, 'PTC OUTPUT WILL BE SAVED TO FILE: ', a, &
541 &/4x, 'OPENED ON UNIT: ', I7)"
542 character(len=*),
parameter :: fmterrasym = &
543 "(a,' **',a,'** PRODUCES AN ASYMMETRIC COEFFICIENT MATRIX, BUT THE &
544 &CONJUGATE GRADIENT METHOD WAS SELECTED. USE BICGSTAB INSTEAD. ')"
547 WRITE (
iout, 1) this%iu
548 00001
FORMAT(1x, /1x,
'IMS -- ITERATIVE MODEL SOLUTION PACKAGE, VERSION 6', &
549 ', 4/28/2017', /, 9x,
'INPUT READ FROM UNIT', i5)
558 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, &
559 supportopenclose=.true., blockrequired=.false.)
563 write (
iout,
'(/1x,a)')
'PROCESSING IMS OPTIONS'
565 call this%parser%GetNextLine(endofblock)
567 call this%parser%GetStringCaps(keyword)
568 select case (keyword)
569 case (
'PRINT_OPTION')
570 call this%parser%GetStringCaps(keyword)
571 if (keyword .eq.
'NONE')
then
573 else if (keyword .eq.
'SUMMARY')
then
575 else if (keyword .eq.
'ALL')
then
578 write (errmsg,
'(3a)') &
579 'Unknown IMS print option (', trim(keyword),
').'
583 call this%parser%GetStringCaps(keyword)
584 if (keyword .eq.
'SIMPLE')
then
587 else if (keyword .eq.
'MODERATE')
then
590 else if (keyword .eq.
'COMPLEX')
then
594 write (errmsg,
'(3a)') &
595 'Unknown IMS COMPLEXITY option (', trim(keyword),
').'
598 case (
'CSV_OUTER_OUTPUT')
599 call this%parser%GetStringCaps(keyword)
600 if (keyword ==
'FILEOUT')
then
601 call this%parser%GetString(fname)
602 if (nr_procs > 1)
then
603 call append_processor_id(fname, proc_id)
606 call openfile(this%icsvouterout,
iout, fname,
'CSV_OUTER_OUTPUT', &
607 filstat_opt=
'REPLACE')
608 write (
iout, fmtcsvout) trim(fname), this%icsvouterout
610 write (errmsg,
'(a)')
'Optional CSV_OUTER_OUTPUT '// &
611 'keyword must be followed by FILEOUT'
614 case (
'CSV_INNER_OUTPUT')
615 call this%parser%GetStringCaps(keyword)
616 if (keyword ==
'FILEOUT')
then
617 call this%parser%GetString(fname)
618 if (nr_procs > 1)
then
619 call append_processor_id(fname, proc_id)
622 call openfile(this%icsvinnerout,
iout, fname,
'CSV_INNER_OUTPUT', &
623 filstat_opt=
'REPLACE')
624 write (
iout, fmtcsvout) trim(fname), this%icsvinnerout
626 write (errmsg,
'(a)')
'Optional CSV_INNER_OUTPUT '// &
627 'keyword must be followed by FILEOUT'
631 call this%parser%GetStringCaps(keyword)
632 select case (keyword)
643 this%iallowptc = ival
644 write (
iout,
'(1x,A)')
'PSEUDO-TRANSIENT CONTINUATION DISABLED FOR'// &
645 ' '//trim(adjustl(msg))//
' STRESS-PERIOD(S)'
646 case (
'ATS_OUTER_MAXIMUM_FRACTION')
647 rval = this%parser%GetDouble()
649 write (errmsg,
'(a,G0)')
'Value for ATS_OUTER_MAXIMUM_FRAC must be &
650 &between 0 and 0.5. Found ', rval
654 write (
iout,
'(1x,A,G0)')
'ADAPTIVE TIME STEP SETTING FOUND. FRACTION &
655 &OF OUTER MAXIMUM USED TO INCREASE OR DECREASE TIME STEP SIZE IS ',&
660 call this%parser%GetStringCaps(keyword)
661 if (keyword ==
'FILEOUT')
then
662 call this%parser%GetString(fname)
664 call openfile(this%icsvouterout,
iout, fname,
'CSV_OUTPUT', &
665 filstat_opt=
'REPLACE')
666 write (
iout, fmtcsvout) trim(fname), this%icsvouterout
669 write (warnmsg,
'(a)') &
670 'OUTER ITERATION INFORMATION WILL BE SAVED TO '//trim(fname)
674 warnmsg, this%parser%GetUnit())
676 write (errmsg,
'(a)')
'Optional CSV_OUTPUT '// &
677 'keyword must be followed by FILEOUT'
686 call this%parser%DevOpt()
688 write (
iout,
'(1x,A)')
'PSEUDO-TRANSIENT CONTINUATION ENABLED'
689 case (
'DEV_PTC_OUTPUT')
690 call this%parser%DevOpt()
692 call this%parser%GetStringCaps(keyword)
693 if (keyword ==
'FILEOUT')
then
694 call this%parser%GetString(fname)
697 filstat_opt=
'REPLACE')
698 write (
iout, fmtptcout) trim(fname), this%iptcout
700 write (errmsg,
'(a)') &
701 'Optional PTC_OUTPUT keyword must be followed by FILEOUT'
704 case (
'DEV_PTC_OPTION')
705 call this%parser%DevOpt()
708 write (
iout,
'(1x,A)') &
709 'PSEUDO-TRANSIENT CONTINUATION USES BNORM AND L2NORM TO '// &
711 case (
'DEV_PTC_EXPONENT')
712 call this%parser%DevOpt()
713 rval = this%parser%GetDouble()
714 if (rval <
dzero)
then
715 write (errmsg,
'(a)')
'PTC_EXPONENT must be > 0.'
720 write (
iout,
'(1x,A,1x,g15.7)') &
721 'PSEUDO-TRANSIENT CONTINUATION EXPONENT', this%ptcexp
723 case (
'DEV_PTC_DEL0')
724 call this%parser%DevOpt()
725 rval = this%parser%GetDouble()
726 if (rval <
dzero)
then
727 write (errmsg,
'(a)')
'IMS sln_ar: PTC_DEL0 must be > 0.'
732 write (
iout,
'(1x,A,1x,g15.7)') &
733 'PSEUDO-TRANSIENT CONTINUATION INITIAL TIMESTEP', this%ptcdel0
736 write (errmsg,
'(a,2(1x,a))') &
737 'Unknown IMS option (', trim(keyword),
').'
741 write (
iout,
'(1x,a)')
'END OF IMS OPTIONS'
743 write (
iout,
'(1x,a)')
'NO IMS OPTION BLOCK DETECTED.'
746 00021
FORMAT(1x,
'SIMPLE OPTION:', /, &
747 1x,
'DEFAULT SOLVER INPUT VALUES FOR FAST SOLUTIONS')
748 00023
FORMAT(1x,
'MODERATE OPTION:', /, 1x,
'DEFAULT SOLVER', &
749 ' INPUT VALUES REFLECT MODERETELY NONLINEAR MODEL')
750 00025
FORMAT(1x,
'COMPLEX OPTION:', /, 1x,
'DEFAULT SOLVER', &
751 ' INPUT VALUES REFLECT STRONGLY NONLINEAR MODEL')
755 call this%sln_setouter(ifdparam)
758 call this%parser%GetBlock(
'NONLINEAR', isfound, ierr, &
759 supportopenclose=.true., blockrequired=.false.)
763 write (
iout,
'(/1x,a)')
'PROCESSING IMS NONLINEAR'
765 call this%parser%GetNextLine(endofblock)
767 call this%parser%GetStringCaps(keyword)
769 select case (keyword)
770 case (
'OUTER_DVCLOSE')
771 this%dvclose = this%parser%GetDouble()
772 case (
'OUTER_MAXIMUM')
773 this%mxiter = this%parser%GetInteger()
774 case (
'UNDER_RELAXATION')
775 call this%parser%GetStringCaps(keyword)
777 if (keyword ==
'NONE')
then
779 else if (keyword ==
'SIMPLE')
then
781 else if (keyword ==
'COOLEY')
then
783 else if (keyword ==
'DBD')
then
786 write (errmsg,
'(3a)') &
787 'Unknown UNDER_RELAXATION specified (', trim(keyword),
').'
791 case (
'LINEAR_SOLVER')
792 call this%parser%GetStringCaps(keyword)
794 if (keyword .eq.
'DEFAULT' .or. &
795 keyword .eq.
'LINEAR')
then
798 write (errmsg,
'(3a)') &
799 'Unknown LINEAR_SOLVER specified (', trim(keyword),
').'
802 this%linsolver = ival
803 case (
'UNDER_RELAXATION_THETA')
804 this%theta = this%parser%GetDouble()
805 case (
'UNDER_RELAXATION_KAPPA')
806 this%akappa = this%parser%GetDouble()
807 case (
'UNDER_RELAXATION_GAMMA')
808 this%gamma = this%parser%GetDouble()
809 case (
'UNDER_RELAXATION_MOMENTUM')
810 this%amomentum = this%parser%GetDouble()
811 case (
'BACKTRACKING_NUMBER')
812 this%numtrack = this%parser%GetInteger()
813 IF (this%numtrack > 0) this%ibflag = 1
814 case (
'BACKTRACKING_TOLERANCE')
815 this%btol = this%parser%GetDouble()
816 case (
'BACKTRACKING_REDUCTION_FACTOR')
817 this%breduc = this%parser%GetDouble()
818 case (
'BACKTRACKING_RESIDUAL_LIMIT')
819 this%res_lim = this%parser%GetDouble()
822 case (
'OUTER_HCLOSE')
823 this%dvclose = this%parser%GetDouble()
826 write (warnmsg,
'(a)') &
827 'SETTING OUTER_DVCLOSE TO OUTER_HCLOSE VALUE'
831 warnmsg, this%parser%GetUnit())
832 case (
'OUTER_RCLOSEBND')
835 write (warnmsg,
'(a)') &
836 'OUTER_DVCLOSE IS USED TO EVALUATE PACKAGE CONVERGENCE'
840 warnmsg, this%parser%GetUnit())
842 write (errmsg,
'(3a)') &
843 'Unknown IMS NONLINEAR keyword (', trim(keyword),
').'
847 write (
iout,
'(1x,a)')
'END OF IMS NONLINEAR DATA'
849 if (ifdparam .EQ. 0)
then
850 write (errmsg,
'(a)')
'NO IMS NONLINEAR block detected.'
855 if (this%theta <
dem3)
then
860 if (this%nonmeth < 1)
then
865 if (this%mxiter <= 0)
then
866 write (errmsg,
'(a)')
'Outer iteration number must be > 0.'
871 if (this%nonmeth > 0)
then
872 WRITE (
iout, *)
'**UNDER-RELAXATION WILL BE USED***'
874 elseif (this%nonmeth == 0)
then
875 WRITE (
iout, *)
'***UNDER-RELAXATION WILL NOT BE USED***'
878 WRITE (errmsg,
'(a)') &
879 'Incorrect value for variable NONMETH was specified.'
884 if (this%nonmeth == 1)
then
885 if (this%gamma == 0)
then
886 WRITE (errmsg,
'(a)') &
887 'GAMMA must be greater than zero if SIMPLE under-relaxation is used.'
892 if (this%solver_mode ==
'PETSC')
then
897 call this%linear_settings%init(this%memory_path)
898 call this%linear_settings%preset_config(ifdparam)
899 call this%linear_settings%read_from_file(this%parser,
iout)
901 if (this%linear_settings%ilinmeth ==
cg_method)
then
907 if (this%solver_mode ==
"IMS")
then
908 allocate (this%imslinear)
909 WRITE (
iout, *)
'***IMS LINEAR SOLVER WILL BE USED***'
910 call this%imslinear%imslinear_allocate(this%name,
iout, this%iprims, &
911 this%mxiter, this%neq, &
912 this%system_matrix, this%rhs, &
913 this%x, this%linear_settings)
916 else if (this%solver_mode ==
"PETSC")
then
917 call this%linear_solver%initialize(this%system_matrix, &
918 this%linear_settings, &
923 write (errmsg,
'(a)') &
924 'Incorrect value for linear solution method specified.'
929 if (this%isymmetric == 1)
then
930 write (
iout,
'(1x,a,/)')
'A symmetric matrix will be solved'
932 write (
iout,
'(1x,a,/)')
'An asymmetric matrix will be solved'
937 if (this%isymmetric == 1)
then
940 do i = 1, this%modellist%Count()
942 if (mp%get_iasym() /= 0)
then
943 write (errmsg, fmterrasym)
'MODEL', trim(adjustl(mp%name))
949 do i = 1, this%exchangelist%Count()
951 if (cp%get_iasym() /= 0)
then
952 write (errmsg, fmterrasym)
'EXCHANGE', trim(adjustl(cp%name))
962 WRITE (
iout, 9002) this%dvclose, this%mxiter, &
963 this%iprims, this%nonmeth, this%linsolver
966 9002
FORMAT(1x,
'OUTER ITERATION CONVERGENCE CRITERION (DVCLOSE) = ', e15.6, &
967 /1x,
'MAXIMUM NUMBER OF OUTER ITERATIONS (MXITER) = ', i0, &
968 /1x,
'SOLVER PRINTOUT INDEX (IPRIMS) = ', i0, &
969 /1x,
'NONLINEAR ITERATION METHOD (NONLINMETH) = ', i0, &
970 /1x,
'LINEAR SOLUTION METHOD (LINMETH) = ', i0)
972 if (this%nonmeth == 1)
then
973 write (
iout, 9003) this%gamma
974 else if (this%nonmeth == 2)
then
975 write (
iout, 9004) this%gamma
976 else if (this%nonmeth == 3)
then
977 write (
iout, 9005) this%theta, this%akappa, this%gamma, this%amomentum
981 if (this%numtrack /= 0)
write (
iout, 9006) this%numtrack, this%btol, &
982 this%breduc, this%res_lim
985 9003
FORMAT(1x,
'UNDER-RELAXATION FACTOR (GAMMA) = ', e15.6)
986 9004
FORMAT(1x,
'UNDER-RELAXATION PREVIOUS HISTORY FACTOR (GAMMA) = ', e15.6)
987 9005
FORMAT(1x,
'UNDER-RELAXATION WEIGHT REDUCTION FACTOR (THETA) = ', e15.6, &
988 /1x,
'UNDER-RELAXATION WEIGHT INCREASE INCREMENT (KAPPA) = ', e15.6, &
989 /1x,
'UNDER-RELAXATION PREVIOUS HISTORY FACTOR (GAMMA) = ', e15.6, &
990 /1x,
'UNDER-RELAXATION MOMENTUM TERM (AMOMENTUM) = ', e15.6)
993 9006
FORMAT(1x,
'MAXIMUM NUMBER OF BACKTRACKS (NUMTRACK) = ', i0, &
994 /1x,
'BACKTRACKING TOLERANCE FACTOR (BTOL) = ', e15.6, &
995 /1x,
'BACKTRACKING REDUCTION FACTOR (BREDUC) = ', e15.6, &
996 /1x,
'BACKTRACKING RESIDUAL LIMIT (RES_LIM) = ', e15.6)
1000 call this%imslinear%imslinear_summary(this%mxiter)
1002 call this%linear_solver%print_summary()
1008 call this%parser%StoreErrorUnit()
1013 call mem_reallocate(this%lrch, 3, this%mxiter,
'LRCH', this%name)
1016 if (this%nonmeth == 3)
then
1021 this%wsave(i) =
dzero
1022 this%hchold(i) =
dzero
1023 this%deold(i) =
dzero
1030 if (this%iprims == 2 .or. this%icsvinnerout > 0)
then
1031 this%nitermax = this%linear_settings%iter1 * this%mxiter
1036 allocate (this%caccel(this%nitermax))
1040 call this%cnvg_summary%reinit(this%nitermax)
1045 call this%parser%StoreErrorUnit()
1049 call this%parser%Clear()
1068 integer(I4B) :: idir
1069 real(DP) :: delt_temp
1070 real(DP) :: fact_lower
1071 real(DP) :: fact_upper
1077 if (this%atsfrac > dzero)
then
1079 fact_lower = this%mxiter * this%atsfrac
1080 fact_upper = this%mxiter - fact_lower
1081 if (this%iouttot_timestep < int(fact_lower))
then
1084 else if (this%iouttot_timestep > int(fact_upper))
then
1111 if (
kper == 1 .and.
kstp == 1)
then
1112 call this%writeCSVHeader()
1116 call this%writePTCInfoToFile(
kper)
1120 this%itertot_timestep = 0
1121 this%iouttot_timestep = 0
1153 write (
iout,
'(//1x,a,1x,a,1x,a)') &
1154 'Solution', trim(adjustl(this%name)),
'summary'
1155 write (
iout,
"(1x,70('-'))")
1156 write (
iout,
'(1x,a,1x,g0,1x,a)') &
1157 'Total formulate time: ', this%ttform,
'seconds'
1158 write (
iout,
'(1x,a,1x,g0,1x,a,/)') &
1159 'Total solution time: ', this%ttsoln,
'seconds'
1179 call this%imslinear%imslinear_da()
1180 deallocate (this%imslinear)
1184 call this%modellist%Clear()
1185 call this%exchangelist%Clear()
1186 deallocate (this%modellist)
1187 deallocate (this%exchangelist)
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)
1198 deallocate (this%caccel)
1201 if (
associated(this%innertab))
then
1202 call this%innertab%table_da()
1203 deallocate (this%innertab)
1204 nullify (this%innertab)
1208 if (
associated(this%outertab))
then
1209 call this%outertab%table_da()
1210 deallocate (this%outertab)
1211 nullify (this%outertab)
1226 call this%cnvg_summary%destroy()
1227 deallocate (this%cnvg_summary)
1230 call this%linear_solver%destroy()
1231 deallocate (this%linear_solver)
1234 call this%linear_settings%destroy()
1235 deallocate (this%linear_settings)
1290 subroutine sln_ca(this, isgcnvg, isuppress_output)
1293 integer(I4B),
intent(inout) :: isgcnvg
1294 integer(I4B),
intent(in) :: isuppress_output
1297 character(len=LINELENGTH) :: line
1298 character(len=LINELENGTH) :: fmt
1300 integer(I4B) :: kiter
1304 call this%prepareSolve()
1308 line =
'mode="validation" -- Skipping matrix assembly and solution.'
1310 do im = 1, this%modellist%Count()
1312 call mp%model_message(line, fmt=fmt)
1316 outerloop:
do kiter = 1, this%mxiter
1319 call this%solve(kiter)
1322 if (this%icnvg == 1)
then
1329 call this%finalizeSolve(kiter, isgcnvg, isuppress_output)
1349 if (this%icsvouterout > 0)
then
1350 write (this%icsvouterout,
'(*(G0,:,","))') &
1351 'total_inner_iterations',
'totim',
'kper',
'kstp',
'nouter', &
1352 'inner_iterations',
'solution_outer_dvmax', &
1353 'solution_outer_dvmax_model',
'solution_outer_dvmax_package', &
1354 'solution_outer_dvmax_node'
1358 if (this%icsvinnerout > 0)
then
1359 write (this%icsvinnerout,
'(*(G0,:,","))', advance=
'NO') &
1360 'total_inner_iterations',
'totim',
'kper',
'kstp',
'nouter', &
1361 'ninner',
'solution_inner_dvmax',
'solution_inner_dvmax_model', &
1362 'solution_inner_dvmax_node'
1363 write (this%icsvinnerout,
'(*(G0,:,","))', advance=
'NO') &
1364 '',
'solution_inner_rmax',
'solution_inner_rmax_model', &
1365 'solution_inner_rmax_node'
1368 write (this%icsvinnerout,
'(*(G0,:,","))', advance=
'NO') &
1369 '',
'solution_inner_alpha'
1370 if (this%imslinear%ilinmeth == 2)
then
1371 write (this%icsvinnerout,
'(*(G0,:,","))', advance=
'NO') &
1372 '',
'solution_inner_omega'
1376 if (this%convnmod > 1)
then
1377 do im = 1, this%modellist%Count()
1379 write (this%icsvinnerout,
'(*(G0,:,","))', advance=
'NO') &
1380 '', trim(adjustl(mp%name))//
'_inner_dvmax', &
1381 trim(adjustl(mp%name))//
'_inner_dvmax_node', &
1382 trim(adjustl(mp%name))//
'_inner_rmax', &
1383 trim(adjustl(mp%name))//
'_inner_rmax_node'
1386 write (this%icsvinnerout,
'(a)')
''
1401 integer(I4B),
intent(in) :: kper
1403 integer(I4B) :: n, im, iallowptc, iptc
1408 do im = 1, this%modellist%Count()
1412 if (this%iallowptc < 0)
then
1420 iallowptc = this%iallowptc
1423 if (iallowptc > 0)
then
1425 call mp%model_ptcchk(iptc)
1432 write (
iout,
'(//)')
1435 write (
iout,
'(1x,a,1x,i0,1x,3a)') &
1436 'PSEUDO-TRANSIENT CONTINUATION WILL BE APPLIED TO MODEL', im,
'("', &
1437 trim(adjustl(mp%name)),
'") DURING THIS TIME STEP'
1461 do ic = 1, this%exchangelist%Count()
1467 do im = 1, this%modellist%Count()
1492 integer(I4B),
intent(in) :: kiter
1496 character(len=LINELENGTH) :: title
1497 character(len=LINELENGTH) :: tag
1498 character(len=LENPAKLOC) :: cmod
1499 character(len=LENPAKLOC) :: cpak
1500 character(len=LENPAKLOC) :: cpakout
1501 character(len=LENPAKLOC) :: strh
1502 character(len=25) :: cval
1503 character(len=7) :: cmsg
1505 integer(I4B) :: im, m_idx, model_id
1506 integer(I4B) :: icsv0
1507 integer(I4B) :: kcsv0
1508 integer(I4B) :: ntabrows
1509 integer(I4B) :: ntabcols
1510 integer(I4B) :: i0, i1
1511 integer(I4B) :: itestmat, n
1512 integer(I4B) :: iter
1513 integer(I4B) :: inewtonur
1514 integer(I4B) :: locmax_nur
1515 integer(I4B) :: iend
1516 integer(I4B) :: icnvgmod
1517 integer(I4B) :: iptc
1518 integer(I4B) :: node_user
1519 integer(I4B) :: ipak
1520 integer(I4B) :: ipos0
1521 integer(I4B) :: ipos1
1522 real(DP) :: dxmax_nur
1523 real(DP) :: dxold_max
1528 real(DP) :: outer_hncg
1531 icsv0 = max(1, this%itertot_sim + 1)
1532 kcsv0 = max(1, this%itertot_timestep + 1)
1535 if (this%iprims > 0)
then
1536 if (.not.
associated(this%outertab))
then
1542 if (this%numtrack > 0)
then
1543 ntabcols = ntabcols + 4
1547 title = trim(this%memory_path)//
' OUTER ITERATION SUMMARY'
1548 call table_cr(this%outertab, this%name, title)
1549 call this%outertab%table_df(ntabrows, ntabcols,
iout, &
1551 tag =
'OUTER ITERATION STEP'
1552 call this%outertab%initialize_column(tag, 25, alignment=
tableft)
1553 tag =
'OUTER ITERATION'
1554 call this%outertab%initialize_column(tag, 10, alignment=
tabright)
1555 tag =
'INNER ITERATION'
1556 call this%outertab%initialize_column(tag, 10, alignment=
tabright)
1557 if (this%numtrack > 0)
then
1558 tag =
'BACKTRACK FLAG'
1559 call this%outertab%initialize_column(tag, 10, alignment=
tabright)
1560 tag =
'BACKTRACK ITERATIONS'
1561 call this%outertab%initialize_column(tag, 10, alignment=
tabright)
1562 tag =
'INCOMING RESIDUAL'
1563 call this%outertab%initialize_column(tag, 15, alignment=
tabright)
1564 tag =
'OUTGOING RESIDUAL'
1565 call this%outertab%initialize_column(tag, 15, alignment=
tabright)
1567 tag =
'MAXIMUM CHANGE'
1568 call this%outertab%initialize_column(tag, 15, alignment=
tabright)
1569 tag =
'STEP SUCCESS'
1570 call this%outertab%initialize_column(tag, 7, alignment=
tabright)
1571 tag =
'MAXIMUM CHANGE MODEL-(CELLID) OR MODEL-PACKAGE-(NUMBER)'
1572 call this%outertab%initialize_column(tag, 34, alignment=
tabright)
1577 if (this%numtrack > 0)
then
1578 call this%sln_backtracking(mp, cp, kiter)
1584 call this%sln_buildsystem(kiter, inewton=1)
1587 call this%sln_calc_ptc(iptc, ptcf)
1590 do im = 1, this%modellist%Count()
1592 call mp%model_nr(kiter, this%system_matrix, 1)
1598 call this%sln_ls(kiter,
kstp,
kper, iter, iptc, ptcf)
1604 this%itertot_timestep = this%itertot_timestep + iter
1605 this%iouttot_timestep = this%iouttot_timestep + 1
1606 this%itertot_sim = this%itertot_sim + iter
1612 if (itestmat /= 0)
then
1613 open (99, file=
'sol_MF6.TXT')
1614 WRITE (99, *)
'MATRIX SOLUTION FOLLOWS'
1615 WRITE (99,
'(10(I8,G15.4))') (n, this%x(n), n=1, this%NEQ)
1622 call this%sln_get_dxmax(this%hncg(kiter), this%lrch(1, kiter))
1624 if (this%sln_has_converged(this%hncg(kiter)))
then
1629 if (this%icnvg == 0)
then
1637 if (kiter == this%mxiter)
then
1642 if (this%iprims > 0)
then
1644 call this%sln_get_loc(this%lrch(1, kiter), strh)
1647 call this%outertab%add_term(cval)
1648 call this%outertab%add_term(kiter)
1649 call this%outertab%add_term(iter)
1650 if (this%numtrack > 0)
then
1651 call this%outertab%add_term(
' ')
1652 call this%outertab%add_term(
' ')
1653 call this%outertab%add_term(
' ')
1654 call this%outertab%add_term(
' ')
1656 call this%outertab%add_term(this%hncg(kiter))
1657 call this%outertab%add_term(cmsg)
1658 call this%outertab%add_term(trim(strh))
1662 do ic = 1, this%exchangelist%Count()
1664 call cp%exg_cc(this%icnvg)
1668 icnvgmod = this%icnvg
1672 do im = 1, this%modellist%Count()
1674 call mp%get_mcellid(0, cmod)
1675 call mp%model_cc(this%itertot_sim, kiter, iend, icnvgmod, &
1678 ipos0 = index(cpak,
'-', back=.true.)
1679 ipos1 = len_trim(cpak)
1680 write (cpakout,
'(a,a,"-(",i0,")",a)') &
1681 trim(cmod), cpak(1:ipos0 - 1), ipak, cpak(ipos0:ipos1)
1688 if (this%icnvg == 1)
then
1689 this%icnvg = this%sln_package_convergence(dpak, cpakout, iend)
1692 if (this%iprims > 0)
then
1694 if (this%icnvg /= 1)
then
1699 if (len_trim(cpakout) > 0)
then
1702 call this%outertab%add_term(cval)
1703 call this%outertab%add_term(kiter)
1704 call this%outertab%add_term(
' ')
1705 if (this%numtrack > 0)
then
1706 call this%outertab%add_term(
' ')
1707 call this%outertab%add_term(
' ')
1708 call this%outertab%add_term(
' ')
1709 call this%outertab%add_term(
' ')
1711 call this%outertab%add_term(dpak)
1712 call this%outertab%add_term(cmsg)
1713 call this%outertab%add_term(cpakout)
1719 if (this%icnvg /= 1)
then
1720 if (this%nonmeth > 0)
then
1721 call this%sln_underrelax(kiter, this%hncg(kiter), this%neq, &
1722 this%active, this%x, this%xtemp)
1724 call this%sln_calcdx(this%neq, this%active, &
1725 this%x, this%xtemp, this%dxold)
1732 do im = 1, this%modellist%Count()
1734 i0 = mp%moffset + 1 - this%matrix_offset
1735 i1 = i0 + mp%neq - 1
1736 call mp%model_nur(mp%neq, this%x(i0:i1), this%xtemp(i0:i1), &
1737 this%dxold(i0:i1), inewtonur, dxmax_nur, locmax_nur)
1741 inewtonur = this%sln_sync_newtonur_flag(inewtonur)
1744 if (inewtonur /= 0)
then
1748 call this%sln_maxval(this%neq, this%dxold, dxold_max)
1751 if (this%sln_nur_has_converged(dxold_max, this%hncg(kiter)))
then
1757 call this%sln_get_dxmax(this%hncg(kiter), this%lrch(1, kiter))
1761 if (this%iprims > 0)
then
1762 cval =
'Newton under-relaxation'
1764 call this%sln_get_loc(this%lrch(1, kiter), strh)
1767 call this%outertab%add_term(cval)
1768 call this%outertab%add_term(kiter)
1769 call this%outertab%add_term(iter)
1770 if (this%numtrack > 0)
then
1771 call this%outertab%add_term(
' ')
1772 call this%outertab%add_term(
' ')
1773 call this%outertab%add_term(
' ')
1774 call this%outertab%add_term(
' ')
1776 call this%outertab%add_term(this%hncg(kiter))
1777 call this%outertab%add_term(cmsg)
1778 call this%outertab%add_term(trim(strh))
1785 if (this%icsvouterout > 0)
then
1788 outer_hncg = this%hncg(kiter)
1791 if (abs(outer_hncg) > abs(dpak))
then
1794 call this%sln_get_nodeu(this%lrch(1, kiter), m_idx, node_user)
1796 else if (outer_hncg ==
dzero .and. dpak ==
dzero)
then
1806 ipos0 = index(cmod,
'_')
1807 read (cmod(1:ipos0 - 1), *) m_idx
1809 ipos0 = index(cpak,
'-', back=.true.)
1810 cpakout = cpak(1:ipos0 - 1)
1820 write (this%icsvouterout,
'(*(G0,:,","))') &
1822 outer_hncg, model_id, trim(cpakout), node_user
1826 if (this%icsvinnerout > 0)
then
1827 call this%csv_convergence_summary(this%icsvinnerout,
totim,
kper,
kstp, &
1828 kiter, iter, icsv0, kcsv0)
1831 end subroutine solve
1843 integer(I4B),
intent(in) :: kiter
1844 integer(I4B),
intent(inout) :: isgcnvg
1845 integer(I4B),
intent(in) :: isuppress_output
1847 integer(I4B) :: ic, im
1851 character(len=*),
parameter :: fmtnocnvg = &
1852 "(1X,'Solution ', i0, ' did not converge for stress period ', i0, &
1853 &' and time step ', i0)"
1854 character(len=*),
parameter :: fmtcnvg = &
1855 "(1X, I0, ' CALLS TO NUMERICAL SOLUTION ', 'IN TIME STEP ', I0, &
1856 &' STRESS PERIOD ',I0,/1X,I0,' TOTAL ITERATIONS')"
1859 if (this%iprims > 0)
then
1860 call this%outertab%finalize_table()
1866 if (this%icnvg /= 0)
then
1867 if (this%iprims > 0)
then
1868 write (
iout, fmtcnvg) kiter,
kstp,
kper, this%itertot_timestep
1877 if (this%iprims == 2)
then
1880 do im = 1, this%modellist%Count()
1882 call this%convergence_summary(mp%iout, im, this%itertot_timestep)
1886 call this%convergence_summary(
iout, this%convnmod + 1, &
1887 this%itertot_timestep)
1891 if (this%icnvg == 0) isgcnvg = 0
1894 do im = 1, this%modellist%Count()
1896 call mp%model_cq(this%icnvg, isuppress_output)
1900 do ic = 1, this%exchangelist%Count()
1902 call cp%exg_cq(isgcnvg, isuppress_output, this%id)
1906 do im = 1, this%modellist%Count()
1908 call mp%model_bd(this%icnvg, isuppress_output)
1912 do ic = 1, this%exchangelist%Count()
1914 call cp%exg_bd(isgcnvg, isuppress_output, this%id)
1922 integer(I4B),
intent(in) :: kiter
1923 integer(I4B),
intent(in) :: inewton
1925 integer(I4B) :: im, ic
1930 call this%sln_reset()
1933 do im = 1, this%modellist%Count()
1935 call mp%model_reset()
1943 do ic = 1, this%exchangelist%Count()
1945 call cp%exg_cf(kiter)
1949 do im = 1, this%modellist%Count()
1951 call mp%model_cf(kiter)
1959 do ic = 1, this%exchangelist%Count()
1961 call cp%exg_fc(kiter, this%system_matrix, this%rhs, inewton)
1965 do im = 1, this%modellist%Count()
1967 call mp%model_fc(kiter, this%system_matrix, inewton)
1982 integer(I4B),
intent(in) :: iu
1983 integer(I4B),
intent(in) :: im
1984 integer(I4B),
intent(in) :: itertot_timestep
1986 character(len=LINELENGTH) :: title
1987 character(len=LINELENGTH) :: tag
1988 character(len=LENPAKLOC) :: loc_dvmax_str
1989 character(len=LENPAKLOC) :: loc_rmax_str
1990 integer(I4B) :: ntabrows
1991 integer(I4B) :: ntabcols
1992 integer(I4B) :: iinner
1994 integer(I4B) :: iouter
1997 integer(I4B) :: locdv
1998 integer(I4B) :: locdr
2010 if (.not.
associated(this%innertab))
then
2014 ntabrows = itertot_timestep
2018 title = trim(this%memory_path)//
' INNER ITERATION SUMMARY'
2019 call table_cr(this%innertab, this%name, title)
2020 call this%innertab%table_df(ntabrows, ntabcols, iu)
2021 tag =
'TOTAL ITERATION'
2022 call this%innertab%initialize_column(tag, 10, alignment=
tabright)
2023 tag =
'OUTER ITERATION'
2024 call this%innertab%initialize_column(tag, 10, alignment=
tabright)
2025 tag =
'INNER ITERATION'
2026 call this%innertab%initialize_column(tag, 10, alignment=
tabright)
2027 tag =
'MAXIMUM CHANGE'
2028 call this%innertab%initialize_column(tag, 15, alignment=
tabright)
2029 tag =
'MAXIMUM CHANGE MODEL-(CELLID)'
2031 tag =
'MAXIMUM RESIDUAL'
2032 call this%innertab%initialize_column(tag, 15, alignment=
tabright)
2033 tag =
'MAXIMUM RESIDUAL MODEL-(CELLID)'
2038 call this%innertab%set_maxbound(itertot_timestep)
2039 call this%innertab%set_iout(iu)
2044 do k = 1, itertot_timestep
2045 iinner = this%cnvg_summary%itinner(k)
2046 if (iinner <= i0)
then
2049 if (im > this%convnmod)
then
2052 do j = 1, this%convnmod
2053 if (abs(this%cnvg_summary%convdvmax(j, k)) > abs(dv))
then
2054 locdv = this%cnvg_summary%convlocdv(j, k)
2055 dv = this%cnvg_summary%convdvmax(j, k)
2057 if (abs(this%cnvg_summary%convrmax(j, k)) > abs(res))
then
2058 locdr = this%cnvg_summary%convlocr(j, k)
2059 res = this%cnvg_summary%convrmax(j, k)
2063 locdv = this%cnvg_summary%convlocdv(im, k)
2064 locdr = this%cnvg_summary%convlocr(im, k)
2065 dv = this%cnvg_summary%convdvmax(im, k)
2066 res = this%cnvg_summary%convrmax(im, k)
2068 call this%sln_get_loc(locdv, loc_dvmax_str)
2069 call this%sln_get_loc(locdr, loc_rmax_str)
2072 call this%innertab%add_term(k)
2073 call this%innertab%add_term(iouter)
2074 call this%innertab%add_term(iinner)
2075 call this%innertab%add_term(dv)
2076 call this%innertab%add_term(adjustr(trim(loc_dvmax_str)))
2077 call this%innertab%add_term(res)
2078 call this%innertab%add_term(adjustr(trim(loc_rmax_str)))
2094 niter, istart, kstart)
2099 integer(I4B),
intent(in) :: iu
2100 real(DP),
intent(in) :: totim
2101 integer(I4B),
intent(in) :: kper
2102 integer(I4B),
intent(in) :: kstp
2103 integer(I4B),
intent(in) :: kouter
2104 integer(I4B),
intent(in) :: niter
2105 integer(I4B),
intent(in) :: istart
2106 integer(I4B),
intent(in) :: kstart
2108 integer(I4B) :: itot
2109 integer(I4B) :: m_idx, j, k
2110 integer(I4B) :: kpos
2111 integer(I4B) :: loc_dvmax
2112 integer(I4B) :: loc_rmax
2113 integer(I4B) :: model_id, node_user
2124 kpos = kstart + k - 1
2125 write (iu,
'(*(G0,:,","))', advance=
'NO') &
2126 itot, totim, kper, kstp, kouter, k
2131 do j = 1, this%convnmod
2132 if (abs(this%cnvg_summary%convdvmax(j, kpos)) > abs(dvmax))
then
2133 loc_dvmax = this%cnvg_summary%convlocdv(j, kpos)
2134 dvmax = this%cnvg_summary%convdvmax(j, kpos)
2136 if (abs(this%cnvg_summary%convrmax(j, kpos)) > abs(rmax))
then
2137 loc_rmax = this%cnvg_summary%convlocr(j, kpos)
2138 rmax = this%cnvg_summary%convrmax(j, kpos)
2143 if (dvmax ==
dzero) loc_dvmax = 0
2144 if (rmax ==
dzero) loc_rmax = 0
2147 if (loc_dvmax > 0)
then
2148 call this%sln_get_nodeu(loc_dvmax, m_idx, node_user)
2150 model_id = num_mod%id
2155 write (iu,
'(*(G0,:,","))', advance=
'NO')
'', dvmax, model_id, node_user
2158 if (loc_rmax > 0)
then
2159 call this%sln_get_nodeu(loc_rmax, m_idx, node_user)
2161 model_id = num_mod%id
2166 write (iu,
'(*(G0,:,","))', advance=
'NO')
'', rmax, model_id, node_user
2170 write (iu,
'(*(G0,:,","))', advance=
'NO') &
2171 '', trim(adjustl(this%caccel(kpos)))
2175 if (this%convnmod > 1)
then
2176 do j = 1, this%cnvg_summary%convnmod
2177 loc_dvmax = this%cnvg_summary%convlocdv(j, kpos)
2178 dvmax = this%cnvg_summary%convdvmax(j, kpos)
2179 loc_rmax = this%cnvg_summary%convlocr(j, kpos)
2180 rmax = this%cnvg_summary%convrmax(j, kpos)
2183 if (loc_dvmax > 0)
then
2184 call this%sln_get_nodeu(loc_dvmax, m_idx, node_user)
2188 write (iu,
'(*(G0,:,","))', advance=
'NO')
'', dvmax, node_user
2191 if (loc_rmax > 0)
then
2192 call this%sln_get_nodeu(loc_rmax, m_idx, node_user)
2196 write (iu,
'(*(G0,:,","))', advance=
'NO')
'', rmax, node_user
2201 write (iu,
'(a)')
''
2226 character(len=*),
intent(in) :: filename
2228 integer(I4B) :: inunit
2231 select type (spm => this%system_matrix)
2234 open (unit=inunit, file=filename, status=
'unknown')
2235 write (inunit, *)
'ia'
2236 write (inunit, *) spm%ia
2237 write (inunit, *)
'ja'
2238 write (inunit, *) spm%ja
2239 write (inunit, *)
'amat'
2240 write (inunit, *) spm%amat
2241 write (inunit, *)
'rhs'
2242 write (inunit, *) this%rhs
2243 write (inunit, *)
'x'
2244 write (inunit, *) this%x
2286 models => this%modellist
2303 select type (exchange)
2317 type(
listtype),
pointer :: exchanges
2319 exchanges => this%exchangelist
2344 do im = 1, this%modellist%Count()
2346 call mp%model_ac(this%sparse)
2353 do ic = 1, this%exchangelist%Count()
2355 call cp%exg_ac(this%sparse)
2360 call this%sparse%sort()
2361 call this%system_matrix%init(this%sparse, this%name)
2362 call this%sparse%destroy()
2367 do im = 1, this%modellist%Count()
2369 call mp%model_mc(this%system_matrix)
2373 do ic = 1, this%exchangelist%Count()
2375 call cp%exg_mc(this%system_matrix)
2393 call this%system_matrix%zero_entries()
2394 call this%vec_rhs%zero_entries()
2406 subroutine sln_ls(this, kiter, kstp, kper, in_iter, iptc, ptcf)
2409 integer(I4B),
intent(in) :: kiter
2410 integer(I4B),
intent(in) :: kstp
2411 integer(I4B),
intent(in) :: kper
2412 integer(I4B),
intent(inout) :: in_iter
2413 integer(I4B),
intent(inout) :: iptc
2414 real(DP),
intent(in) :: ptcf
2416 logical(LGP) :: lsame
2418 integer(I4B) :: irow_glo
2419 integer(I4B) :: itestmat
2420 integer(I4B) :: ipos
2421 integer(I4B) :: icol_s
2422 integer(I4B) :: icol_e
2423 integer(I4B) :: jcol
2424 integer(I4B) :: iptct
2425 integer(I4B) :: iallowptc
2431 character(len=50) :: fname
2432 character(len=*),
parameter :: fmtfname =
"('mf6mat_', i0, '_', i0, &
2433 &'_', i0, '_', i0, '.txt')"
2436 do ieq = 1, this%neq
2439 irow_glo = ieq + this%matrix_offset
2442 this%xtemp(ieq) = this%x(ieq)
2446 if (this%active(ieq) > 0)
then
2448 adiag = abs(this%system_matrix%get_diag_value(irow_glo))
2449 if (adiag <
dem15)
then
2450 call this%system_matrix%set_diag_value(irow_glo, diagval)
2451 this%rhs(ieq) = this%rhs(ieq) + diagval * this%x(ieq)
2455 call this%system_matrix%set_diag_value(irow_glo,
done)
2456 call this%system_matrix%zero_row_offdiag(irow_glo)
2457 this%rhs(ieq) = this%x(ieq)
2463 do ieq = 1, this%neq
2464 if (this%active(ieq) > 0)
then
2465 icol_s = this%system_matrix%get_first_col_pos(ieq)
2466 icol_e = this%system_matrix%get_last_col_pos(ieq)
2467 do ipos = icol_s, icol_e
2468 jcol = this%system_matrix%get_column(ipos)
2469 if (jcol == ieq) cycle
2470 if (this%active(jcol) < 0)
then
2471 this%rhs(ieq) = this%rhs(ieq) - &
2472 (this%system_matrix%get_value_pos(ipos) * &
2474 call this%system_matrix%set_value_pos(ipos,
dzero)
2486 if (this%iallowptc < 0)
then
2495 iallowptc = this%iallowptc
2499 iptct = iptc * iallowptc
2503 if (iptct /= 0)
then
2504 call this%sln_l2norm(l2norm)
2507 if (kiter == 1)
then
2508 if (kper > 1 .or. kstp > 1)
then
2509 if (l2norm <= this%l2norm0)
then
2514 lsame =
is_close(l2norm, this%l2norm0)
2520 iptct = iptc * iallowptc
2521 if (iptct /= 0)
then
2522 if (kiter == 1)
then
2523 if (this%iptcout > 0)
then
2524 write (this%iptcout,
'(A10,6(1x,A15))')
'OUTER ITER', &
2525 ' PTCDEL',
' L2NORM0',
' L2NORM', &
2526 ' RHSNORM',
' 1/PTCDEL',
' RHSNORM/L2NORM'
2528 if (this%ptcdel0 >
dzero)
then
2529 this%ptcdel = this%ptcdel0
2531 if (this%iptcopt == 0)
then
2534 this%ptcdel =
done / ptcf
2537 do ieq = 1, this%neq
2538 if (this%active(ieq) .gt. 0)
then
2539 bnorm = bnorm + this%rhs(ieq) * this%rhs(ieq)
2543 this%ptcdel = bnorm / l2norm
2547 if (l2norm >
dzero)
then
2548 this%ptcdel = this%ptcdel * (this%l2norm0 / l2norm)**this%ptcexp
2553 if (this%ptcdel >
dzero)
then
2554 ptcval =
done / this%ptcdel
2559 do ieq = 1, this%neq
2560 irow_glo = ieq + this%matrix_offset
2561 if (this%active(ieq) > 0)
then
2562 diagval = abs(this%system_matrix%get_diag_value(irow_glo))
2563 bnorm = bnorm + this%rhs(ieq) * this%rhs(ieq)
2564 call this%system_matrix%add_diag_value(irow_glo, -ptcval)
2565 this%rhs(ieq) = this%rhs(ieq) - ptcval * this%x(ieq)
2569 if (this%iptcout > 0)
then
2570 write (this%iptcout,
'(i10,5(1x,e15.7),1(1x,f15.6))') &
2571 kiter, this%ptcdel, this%l2norm0, l2norm, bnorm, &
2572 ptcval, bnorm / l2norm
2574 this%l2norm0 = l2norm
2581 if (itestmat == 1)
then
2582 write (fname, fmtfname) this%id, kper, kstp, kiter
2583 print *,
'Saving amat to: ', trim(adjustl(fname))
2586 open (itestmat, file=trim(adjustl(fname)))
2587 write (itestmat, *)
'NODE, RHS, AMAT FOLLOW'
2588 do ieq = 1, this%neq
2589 irow_glo = ieq + this%matrix_offset
2590 icol_s = this%system_matrix%get_first_col_pos(irow_glo)
2591 icol_e = this%system_matrix%get_last_col_pos(irow_glo)
2592 write (itestmat,
'(*(G0,:,","))') &
2595 (this%system_matrix%get_column(ipos), ipos=icol_s, icol_e), &
2596 (this%system_matrix%get_value_pos(ipos), ipos=icol_s, icol_e)
2607 call this%imslinear%imslinear_apply(this%icnvg, kstp, kiter, in_iter, &
2608 this%nitermax, this%convnmod, &
2609 this%convmodstart, this%caccel, &
2612 call this%linear_solver%solve(kiter, this%vec_rhs, &
2613 this%vec_x, this%cnvg_summary)
2614 in_iter = this%linear_solver%iteration_number
2615 this%icnvg = this%linear_solver%is_converged
2631 integer(I4B),
intent(in) :: ifdparam
2634 select case (ifdparam)
2642 this%amomentum =
dzero
2646 this%res_lim =
dzero
2654 this%akappa = 0.0001d0
2656 this%amomentum =
dzero
2660 this%res_lim =
dzero
2668 this%akappa = 0.0001d0
2670 this%amomentum =
dzero
2674 this%res_lim = 0.002d0
2693 integer(I4B),
intent(in) :: kiter
2695 character(len=7) :: cmsg
2697 integer(I4B) :: btflag
2698 integer(I4B) :: ibflag
2699 integer(I4B) :: ibtcnt
2707 call this%sln_buildsystem(kiter, inewton=0)
2711 if (kiter == 1)
then
2712 call this%sln_l2norm(this%res_prev)
2713 resin = this%res_prev
2716 call this%sln_l2norm(this%res_new)
2717 resin = this%res_new
2721 if (this%res_new > this%res_prev * this%btol)
then
2724 btloop:
do nb = 1, this%numtrack
2727 call this%sln_backtracking_xupdate(btflag)
2730 if (btflag == 0)
then
2738 call this%sln_buildsystem(kiter, inewton=0)
2742 call this%sln_l2norm(this%res_new)
2745 if (nb == this%numtrack)
then
2749 if (this%res_new < this%res_prev * this%btol)
then
2753 if (this%res_new < this%res_lim)
then
2759 this%res_prev = this%res_new
2763 if (this%iprims > 0)
then
2764 if (ibtcnt > 0)
then
2771 call this%outertab%add_term(
'Backtracking')
2772 call this%outertab%add_term(kiter)
2773 call this%outertab%add_term(
' ')
2774 if (this%numtrack > 0)
then
2775 call this%outertab%add_term(ibflag)
2776 call this%outertab%add_term(ibtcnt)
2777 call this%outertab%add_term(resin)
2778 call this%outertab%add_term(this%res_prev)
2780 call this%outertab%add_term(
' ')
2781 call this%outertab%add_term(cmsg)
2782 call this%outertab%add_term(
' ')
2798 integer(I4B),
intent(inout) :: bt_flag
2800 bt_flag = this%get_backtracking_flag()
2803 if (bt_flag > 0)
then
2804 call this%apply_backtracking()
2813 integer(I4B) :: bt_flag
2818 real(dp) :: dx_abs_max
2826 if (this%active(n) < 1) cycle
2827 dx = this%x(n) - this%xtemp(n)
2829 if (dx_abs > dx_abs_max) dx_abs_max = dx_abs
2833 if (this%breduc * dx_abs_max >= this%dvclose)
then
2848 if (this%active(n) < 1) cycle
2849 delx = this%breduc * (this%x(n) - this%xtemp(n))
2850 this%x(n) = this%xtemp(n) + delx
2873 vec_resid => this%system_matrix%create_vec(this%neq)
2874 call this%sln_calc_residual(vec_resid)
2877 l2norm = vec_resid%norm2()
2880 call vec_resid%destroy()
2881 deallocate (vec_resid)
2894 integer(I4B),
intent(in) :: nsize
2895 real(DP),
dimension(nsize),
intent(in) :: v
2896 real(DP),
intent(inout) :: vmax
2908 if (denom ==
dzero)
then
2913 dnorm = abs(d) / denom
2914 if (dnorm >
done)
then
2931 integer(I4B),
intent(in) :: neq
2932 integer(I4B),
dimension(neq),
intent(in) :: active
2933 real(DP),
dimension(neq),
intent(in) :: x
2934 real(DP),
dimension(neq),
intent(in) :: xtemp
2935 real(DP),
dimension(neq),
intent(inout) :: dx
2942 if (active(n) < 1)
then
2945 dx(n) = x(n) - xtemp(n)
2957 integer(I4B) :: iptc
2968 vec_resid => this%system_matrix%create_vec(this%neq)
2969 call this%sln_calc_residual(vec_resid)
2972 do im = 1, this%modellist%Count()
2974 call mp%model_ptc(vec_resid, iptc, ptcf)
2978 call vec_resid%destroy()
2979 deallocate (vec_resid)
2991 call this%system_matrix%multiply(this%vec_x, vec_resid)
2993 call vec_resid%axpy(-1.0_dp, this%vec_rhs)
2996 if (this%active(n) < 1)
then
2997 call vec_resid%set_value_local(n, 0.0_dp)
3011 integer(I4B),
intent(in) :: kiter
3012 real(DP),
intent(in) :: bigch
3013 integer(I4B),
intent(in) :: neq
3014 integer(I4B),
dimension(neq),
intent(in) :: active
3015 real(DP),
dimension(neq),
intent(inout) :: x
3016 real(DP),
dimension(neq),
intent(in) :: xtemp
3028 if (this%nonmeth == 1)
then
3032 if (active(n) < 1) cycle
3035 delx = x(n) - xtemp(n)
3036 this%dxold(n) = delx
3039 x(n) = xtemp(n) + this%gamma * delx
3043 else if (this%nonmeth == 2)
then
3049 if (kiter == 1)
then
3051 this%relaxold =
done
3052 this%bigchold = bigch
3056 es = this%bigch / (this%bigchold * this%relaxold)
3058 if (es < -
done)
then
3064 this%relaxold = relax
3067 this%bigchold = (
done - this%gamma) * this%bigch + this%gamma * &
3071 if (relax <
done)
then
3075 if (active(n) < 1) cycle
3078 delx = x(n) - xtemp(n)
3079 this%dxold(n) = delx
3080 x(n) = xtemp(n) + relax * delx
3085 else if (this%nonmeth == 3)
then
3089 if (active(n) < 1) cycle
3092 delx = x(n) - xtemp(n)
3095 if (kiter == 1)
then
3096 this%wsave(n) =
done
3097 this%hchold(n) =
dem20
3098 this%deold(n) =
dzero
3105 if (this%deold(n) * delx <
dzero)
then
3106 ww = this%theta * this%wsave(n)
3109 ww = this%wsave(n) + this%akappa
3115 if (kiter == 1)
then
3116 this%hchold(n) = delx
3118 this%hchold(n) = (
done - this%gamma) * delx + &
3119 this%gamma * this%hchold(n)
3123 this%deold(n) = delx
3124 this%dxold(n) = delx
3128 if (kiter > 4) amom = this%amomentum
3129 delx = delx * ww + amom * this%hchold(n)
3130 x(n) = xtemp(n) + delx
3148 real(DP),
intent(inout) :: hncg
3149 integer(I4B),
intent(inout) :: lrch
3163 if (this%active(n) < 1) cycle
3164 hdif = this%x(n) - this%xtemp(n)
3166 if (ahdif > abigch)
then
3184 logical(LGP) :: has_converged
3186 has_converged = .false.
3187 if (abs(max_dvc) <= this%dvclose)
then
3188 has_converged = .true.
3198 real(dp),
intent(in) :: dpak
3199 character(len=LENPAKLOC),
intent(in) :: cpakout
3200 integer(I4B),
intent(in) :: iend
3202 integer(I4B) :: ivalue
3204 if (abs(dpak) > this%dvclose)
then
3209 'PACKAGE (', trim(cpakout),
') CAUSED CONVERGENCE FAILURE'
3221 integer(I4B),
intent(in) :: inewtonur
3223 integer(I4B) :: ivalue
3232 result(has_converged)
3234 real(dp),
intent(in) :: dxold_max
3235 real(dp),
intent(in) :: hncg
3236 logical(LGP) :: has_converged
3238 has_converged = .false.
3239 if (abs(dxold_max) <= this%dvclose .and. &
3240 abs(hncg) <= this%dvclose)
then
3241 has_converged = .true.
3253 integer(I4B),
intent(in) :: nodesln
3254 character(len=*),
intent(inout) :: str
3258 integer(I4B) :: istart
3259 integer(I4B) :: iend
3260 integer(I4B) :: noder
3261 integer(I4B) :: nglo
3270 nglo = nodesln + this%matrix_offset
3273 do i = 1, this%modellist%Count()
3277 call mp%get_mrange(istart, iend)
3278 if (nglo >= istart .and. nglo <= iend)
then
3279 noder = nglo - istart + 1
3280 call mp%get_mcellid(noder, str)
3297 integer(I4B),
intent(in) :: nodesln
3298 integer(I4B),
intent(inout) :: im
3299 integer(I4B),
intent(inout) :: nodeu
3303 integer(I4B) :: istart
3304 integer(I4B) :: iend
3305 integer(I4B) :: noder, nglo
3311 nglo = nodesln + this%matrix_offset
3314 do i = 1, this%modellist%Count()
3318 call mp%get_mrange(istart, iend)
3319 if (nglo >= istart .and. nglo <= iend)
then
3320 noder = nglo - istart + 1
3321 call mp%get_mnodeu(noder, nodeu)
3338 class(*),
pointer,
intent(inout) :: obj
3346 if (.not.
associated(obj))
return
3365 type(
listtype),
intent(inout) :: list
3366 integer(I4B),
intent(in) :: idx
3370 class(*),
pointer :: obj
3372 obj => list%GetItem(idx)
subroutine, public ats_submit_delt(kstp, kper, dt, sloc, idir)
@ brief Allow and external caller to submit preferred time step
subroutine, public addbasesolutiontolist(list, solution)
This module contains block parser methods.
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
@ tabright
right justified table column
@ tableft
left justified table column
@ mvalidate
validation mode - do not run time steps
@ mnormal
normal output mode
real(dp), parameter dem20
real constant 1e-20
real(dp), parameter dep3
real constant 1000
integer(i4b), parameter lensolutionname
maximum length of the solution name
real(dp), parameter dep6
real constant 1000000
real(dp), parameter donethird
real constant 1/3
real(dp), parameter dnodata
real no data constant
integer(i4b), parameter lenpakloc
maximum length of a package location
real(dp), parameter dep20
real constant 1e20
real(dp), parameter dem1
real constant 1e-1
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dem3
real constant 1e-3
real(dp), parameter dem4
real constant 1e-4
real(dp), parameter dem6
real constant 1e-6
real(dp), parameter dzero
real constant zero
real(dp), parameter dprec
real constant machine precision
real(dp), parameter dem15
real constant 1e-15
real(dp), parameter dem2
real constant 1e-2
real(dp), parameter dtwo
real constant 2
integer(i4b), parameter lenmempath
maximum length of the memory path
real(dp), parameter dthree
real constant 3
real(dp), parameter done
real constant 1
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
subroutine allocate_scalars(this)
@ brief Allocate and initialize scalars
integer(i4b), parameter, public cg_method
This module defines variable data types.
class(linearsolverbasetype) function, pointer, public create_linear_solver(solver_mode, sln_name)
Factory method to create the linear solver object.
type(listtype), public basesolutionlist
pure logical function, public is_close(a, b, rtol, atol, symmetric)
Check if a real value is approximately equal to another.
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.
subroutine, public write_message(text, iunit, fmt, skipbefore, skipafter, advance)
Write a message to an output unit.
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 solve(this, kiter)
@ brief Build and solve the simulation
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
This module contains simulation methods.
subroutine, public store_error(msg, terminate)
Store an error message.
integer(i4b) function, public count_errors()
Return number of errors.
subroutine, public deprecation_warning(cblock, cvar, cver, endmsg, iunit)
Store deprecation warning message.
integer(i4b), parameter, public stg_bfr_exg_ad
before exchange advance (per solution)
integer(i4b), parameter, public stg_bfr_exg_cf
before exchange calculate (per solution)
integer(i4b), parameter, public stg_bfr_exg_fc
before exchange formulate (per solution)
integer(i4b), parameter, public stg_bfr_exg_ac
before exchange add connections (per solution)
This module contains simulation variables.
character(len=maxcharlen) errmsg
error message string
character(len=linelength) simulation_mode
integer(i4b) iout
file unit number for simulation output
integer(i4b) isim_mode
simulation mode
subroutine, public table_cr(this, name, title)
real(dp), pointer, public totim
time relative to start of simulation
integer(i4b), pointer, public kstp
current time step number
integer(i4b), pointer, public kper
current stress period number
real(dp), pointer, public delt
length of the current time step
subroutine, public code_timer(it, t1, ts)
Get end time and calculate elapsed time.
This module contains version information.
integer(i4b), parameter idevelopmode
Highest level model type. All models extend this parent type.
This structure stores the generic convergence info for a solution.
Abstract type for linear solver.
A generic heterogeneous doubly-linked list.