MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
ParallelSolution.f90
Go to the documentation of this file.
2  use kindmodule, only: dp, lgp, i4b
6  use mpi
8  implicit none
9  private
10 
11  public :: parallelsolutiontype
12 
14  integer(I4B) :: tmr_convergence = -1 !< timer for convergence check
15  integer(I4B) :: tmr_pkg_cnvg = -1 !< timer for package convergence check
16  integer(I4B) :: tmr_sync_nur = -1 !< timer for NUR synchronization
17  integer(I4B) :: tmr_nur_cnvg = -1 !< timer for NUR convergence check
18  integer(I4B) :: tmr_calcptc = -1 !< timer for PTC calculation
19  integer(I4B) :: tmr_underrelax = -1 !< timer for underrelaxation
20  integer(I4B) :: tmr_backtracking = -1 !< timer for backtracking
21  contains
22  ! override
23  procedure :: sln_has_converged => par_has_converged
24  procedure :: sln_package_convergence => par_package_convergence
25  procedure :: sln_sync_newtonur_flag => par_sync_newtonur_flag
26  procedure :: sln_nur_has_converged => par_nur_has_converged
27  procedure :: sln_calc_ptc => par_calc_ptc
28  procedure :: sln_underrelax => par_underrelax
29  procedure :: sln_backtracking_xupdate => par_backtracking_xupdate
30 
31  end type parallelsolutiontype
32 
33 contains
34 
35  !> @brief Check global convergence. The local maximum dependent
36  !! variable change is reduced over MPI with all other processes
37  !< that are running this parallel numerical solution.
38  function par_has_converged(this, max_dvc) result(has_converged)
39  class(parallelsolutiontype) :: this !< parallel solution
40  real(dp) :: max_dvc !< the LOCAL maximum dependent variable change
41  logical(LGP) :: has_converged !< True, when GLOBALLY converged
42  ! local
43  real(dp) :: global_max_dvc
44  real(dp) :: abs_max_dvc
45  integer :: ierr
46  type(mpiworldtype), pointer :: mpi_world
47 
48  call g_prof%start("Parallel Solution (cnvg check)", this%tmr_convergence)
49 
50  mpi_world => get_mpi_world()
51 
52  has_converged = .false.
53  abs_max_dvc = abs(max_dvc)
54  call mpi_allreduce(abs_max_dvc, global_max_dvc, 1, mpi_double_precision, &
55  mpi_max, mpi_world%comm, ierr)
56  call check_mpi(ierr)
57  if (global_max_dvc <= this%dvclose) then
58  has_converged = .true.
59  end if
60 
61  call g_prof%stop(this%tmr_convergence)
62 
63  end function par_has_converged
64 
65  function par_package_convergence(this, dpak, cpakout, iend) &
66  result(icnvg_global)
67  class(parallelsolutiontype) :: this !< parallel solution instance
68  real(dp), intent(in) :: dpak !< Newton Under-relaxation flag
69  character(len=LENPAKLOC), intent(in) :: cpakout
70  integer(I4B), intent(in) :: iend
71  ! local
72  integer(I4B) :: icnvg_global
73  integer(I4B) :: icnvg_local
74  integer :: ierr
75  type(mpiworldtype), pointer :: mpi_world
76 
77  call g_prof%start("Parallel Solution (package cnvg)", this%tmr_pkg_cnvg)
78 
79  mpi_world => get_mpi_world()
80 
81  icnvg_local = &
82  this%NumericalSolutionType%sln_package_convergence(dpak, cpakout, iend)
83 
84  call mpi_allreduce(icnvg_local, icnvg_global, 1, mpi_integer, &
85  mpi_min, mpi_world%comm, ierr)
86  call check_mpi(ierr)
87 
88  call g_prof%stop(this%tmr_pkg_cnvg)
89 
90  end function par_package_convergence
91 
92  function par_sync_newtonur_flag(this, inewtonur) result(ivalue)
93  class(parallelsolutiontype) :: this !< parallel solution instance
94  integer(I4B), intent(in) :: inewtonur !< local Newton Under-relaxation flag
95  ! local
96  integer(I4B) :: ivalue !< Maximum of all local values (1 = under-relaxation applied)
97  integer :: ierr
98  type(mpiworldtype), pointer :: mpi_world
99 
100  call g_prof%start("Parallel Solution (NUR)", this%tmr_sync_nur)
101 
102  mpi_world => get_mpi_world()
103  call mpi_allreduce(inewtonur, ivalue, 1, mpi_integer, &
104  mpi_max, mpi_world%comm, ierr)
105  call check_mpi(ierr)
106 
107  call g_prof%stop(this%tmr_sync_nur)
108 
109  end function par_sync_newtonur_flag
110 
111  function par_nur_has_converged(this, dxold_max, hncg) &
112  result(has_converged)
113  class(parallelsolutiontype) :: this !< parallel solution instance
114  real(dp), intent(in) :: dxold_max !< the maximum dependent variable change for cells not adjusted by NUR
115  real(dp), intent(in) :: hncg !< largest dep. var. change at end of Picard iter.
116  logical(LGP) :: has_converged !< True, when converged
117  ! local
118  integer(I4B) :: icnvg_local
119  integer(I4B) :: icnvg_global
120  integer :: ierr
121  type(mpiworldtype), pointer :: mpi_world
122 
123  call g_prof%start("Parallel Solution (NUR cnvg)", this%tmr_nur_cnvg)
124 
125  mpi_world => get_mpi_world()
126 
127  has_converged = .false.
128  icnvg_local = 0
129  if (this%NumericalSolutionType%sln_nur_has_converged( &
130  dxold_max, hncg)) then
131  icnvg_local = 1
132  end if
133 
134  call mpi_allreduce(icnvg_local, icnvg_global, 1, mpi_integer, &
135  mpi_min, mpi_world%comm, ierr)
136  call check_mpi(ierr)
137  if (icnvg_global == 1) has_converged = .true.
138 
139  call g_prof%stop(this%tmr_nur_cnvg)
140 
141  end function par_nur_has_converged
142 
143  !> @brief Calculate pseudo-transient continuation factor
144  !< for the parallel case
145  subroutine par_calc_ptc(this, iptc, ptcf)
146  class(parallelsolutiontype) :: this !< parallel solution
147  integer(I4B) :: iptc !< PTC (1) or not (0)
148  real(DP) :: ptcf !< the (global) PTC factor calculated
149  ! local
150  integer(I4B) :: iptc_loc
151  real(DP) :: ptcf_loc, ptcf_glo_max
152  integer :: ierr
153  type(mpiworldtype), pointer :: mpi_world
154 
155  call g_prof%start("Parallel Solution (PTC calc)", this%tmr_calcptc)
156 
157  mpi_world => get_mpi_world()
158  call this%NumericalSolutionType%sln_calc_ptc(iptc_loc, ptcf_loc)
159  if (iptc_loc == 0) ptcf_loc = dzero
160 
161  ! now reduce
162  call mpi_allreduce(ptcf_loc, ptcf_glo_max, 1, mpi_double_precision, &
163  mpi_max, mpi_world%comm, ierr)
164  call check_mpi(ierr)
165 
166  iptc = 0
167  ptcf = dzero
168  if (ptcf_glo_max > dzero) then
169  iptc = 1
170  ptcf = ptcf_glo_max
171  end if
172 
173  call g_prof%stop(this%tmr_calcptc)
174 
175  end subroutine par_calc_ptc
176 
177  !> @brief apply under-relaxation in sync over all processes
178  !<
179  subroutine par_underrelax(this, kiter, bigch, neq, active, x, xtemp)
180  class(parallelsolutiontype) :: this !< parallel instance
181  integer(I4B), intent(in) :: kiter !< Picard iteration number
182  real(DP), intent(in) :: bigch !< maximum dependent-variable change
183  integer(I4B), intent(in) :: neq !< number of equations
184  integer(I4B), dimension(neq), intent(in) :: active !< active cell flag (1)
185  real(DP), dimension(neq), intent(inout) :: x !< current dependent-variable
186  real(DP), dimension(neq), intent(in) :: xtemp !< previous dependent-variable
187  ! local
188  real(DP) :: dvc_global_max, dvc_global_min
189  integer :: ierr
190  type(mpiworldtype), pointer :: mpi_world
191 
192  call g_prof%start("Parallel Solution (underrelax)", this%tmr_underrelax)
193 
194  mpi_world => get_mpi_world()
195 
196  ! first reduce largest change over all processes
197  call mpi_allreduce(bigch, dvc_global_max, 1, mpi_double_precision, &
198  mpi_max, mpi_world%comm, ierr)
199  call check_mpi(ierr)
200  call mpi_allreduce(bigch, dvc_global_min, 1, mpi_double_precision, &
201  mpi_min, mpi_world%comm, ierr)
202  call check_mpi(ierr)
203 
204  if (abs(dvc_global_min) > abs(dvc_global_max)) then
205  dvc_global_max = dvc_global_min
206  end if
207 
208  ! call local underrelax routine
209  call this%NumericalSolutionType%sln_underrelax(kiter, dvc_global_max, &
210  neq, active, x, xtemp)
211 
212  call g_prof%stop(this%tmr_underrelax)
213 
214  end subroutine par_underrelax
215 
216  !> @brief synchronize backtracking flag over processes
217  !<
218  subroutine par_backtracking_xupdate(this, bt_flag)
219  ! -- dummy variables
220  class(parallelsolutiontype), intent(inout) :: this !< ParallelSolutionType instance
221  integer(I4B), intent(inout) :: bt_flag !< global backtracking flag (1) backtracking performed (0) backtracking not performed
222  ! -- local variables
223  integer(I4B) :: btflag_local
224  type(mpiworldtype), pointer :: mpi_world
225  integer :: ierr
226 
227  call g_prof%start("Parallel Solution (backtrack)", this%tmr_backtracking)
228 
229  mpi_world => get_mpi_world()
230 
231  ! get local bt flag
232  btflag_local = this%NumericalSolutionType%get_backtracking_flag()
233 
234  ! reduce into global decision (if any, then all)
235  call mpi_allreduce(btflag_local, bt_flag, 1, mpi_integer, &
236  mpi_max, mpi_world%comm, ierr)
237  call check_mpi(ierr)
238 
239  ! perform backtracking if ...
240  if (bt_flag > 0) then
241  call this%NumericalSolutionType%apply_backtracking()
242  end if
243 
244  call g_prof%stop(this%tmr_backtracking)
245 
246  end subroutine par_backtracking_xupdate
247 
248 end module parallelsolutionmodule
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter lenpakloc
maximum length of a package location
Definition: Constants.f90:50
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
This module defines variable data types.
Definition: kind.f90:8
type(mpiworldtype) function, pointer, public get_mpi_world()
Definition: MpiWorld.f90:32
subroutine, public check_mpi(mpi_error_code)
Check the MPI error code, report, and.
Definition: MpiWorld.f90:116
logical(lgp) function par_has_converged(this, max_dvc)
Check global convergence. The local maximum dependent variable change is reduced over MPI with all ot...
integer(i4b) function par_package_convergence(this, dpak, cpakout, iend)
logical(lgp) function par_nur_has_converged(this, dxold_max, hncg)
integer(i4b) function par_sync_newtonur_flag(this, inewtonur)
subroutine par_underrelax(this, kiter, bigch, neq, active, x, xtemp)
apply under-relaxation in sync over all processes
subroutine par_backtracking_xupdate(this, bt_flag)
synchronize backtracking flag over processes
subroutine par_calc_ptc(this, iptc, ptcf)
Calculate pseudo-transient continuation factor.
type(profilertype), public g_prof
the global timer object (to reduce trivial lines of code)
Definition: Profiler.f90:65