Line data Source code
1 : module gmean_mod
2 : !-----------------------------------------------------------------------
3 : !
4 : ! Purpose:
5 : ! Perform global mean calculations for energy conservation and other checks.
6 : !
7 : ! Method:
8 : ! Reproducible (scalable):
9 : ! Convert to fixed point (integer representation) to enable
10 : ! reproducibility when using MPI collectives.
11 : ! If error checking is on (via setting reprosum_diffmax > 0 and
12 : ! reprosum_recompute = .true. in user_nl_cpl), shr_reprosum_calc will
13 : ! check the accuracy of its computation with a fast but
14 : ! non-reproducible algorithm. If any error is reported, report
15 : ! the difference and the expected sum and abort run (call endrun)
16 : !
17 : !
18 : !-----------------------------------------------------------------------
19 : use shr_kind_mod, only: r8 => shr_kind_r8
20 : use ppgrid, only: pcols, begchunk, endchunk
21 : use shr_reprosum_mod, only: shr_reprosum_calc, shr_reprosum_tolExceeded
22 : use shr_reprosum_mod, only: shr_reprosum_reldiffmax, shr_reprosum_recompute
23 : use perf_mod, only: t_startf, t_stopf
24 : use cam_logfile, only: iulog
25 :
26 : implicit none
27 : private
28 :
29 : public :: gmean ! compute global mean of 2D fields on physics decomposition
30 : public :: gmean_init ! Initialize gmean (maybe run tests)
31 : public :: test_gmean ! test accuracy of gmean
32 :
33 : interface gmean
34 : module procedure gmean_arr
35 : module procedure gmean_scl
36 : end interface gmean
37 :
38 : private :: gmean_fixed_repro
39 : private :: gmean_float_norepro
40 :
41 : ! Set do_gmean_tests to .true. to run a gmean challenge test
42 : logical, private :: do_gmean_tests = .false.
43 :
44 : CONTAINS
45 :
46 : !
47 : !========================================================================
48 : !
49 :
50 0 : subroutine gmean_init(do_test)
51 : !-----------------------------------------------------------------------
52 : !
53 : ! Purpose: Possibly run a test
54 : !
55 : !-----------------------------------------------------------------------
56 : !
57 : logical, optional, intent(in) :: do_test
58 :
59 : logical :: do_test_use
60 :
61 0 : if (present(do_test)) then
62 0 : do_test_use = do_test
63 : else
64 0 : do_test_use = do_gmean_tests
65 : end if
66 :
67 0 : if (do_test_use) then
68 0 : call test_gmean()
69 : end if
70 :
71 0 : end subroutine gmean_init
72 :
73 : !
74 : !========================================================================
75 : !
76 :
77 370944 : subroutine gmean_arr (arr, arr_gmean, nflds)
78 : use shr_strconvert_mod, only: toString
79 : use spmd_utils, only: masterproc
80 : use cam_abortutils, only: endrun
81 : !-----------------------------------------------------------------------
82 : !
83 : ! Purpose:
84 : ! Compute the global mean of each field in "arr" in the physics
85 : ! chunked decomposition
86 : !
87 : ! Method is to call shr_reprosum_calc (called from gmean_fixed_repro)
88 : !-----------------------------------------------------------------------
89 : !
90 : ! Arguments
91 : !
92 : integer, intent(in) :: nflds ! number of fields
93 : real(r8), intent(in) :: arr(pcols, begchunk:endchunk, nflds)
94 : real(r8), intent(out) :: arr_gmean(nflds) ! global means
95 : !
96 : ! Local workspace
97 : !
98 741888 : real(r8) :: rel_diff(2, nflds)
99 : integer :: ifld ! field index
100 : integer :: num_err
101 : logical :: write_warning
102 : !
103 : !-----------------------------------------------------------------------
104 : !
105 370944 : call t_startf('gmean_arr')
106 370944 : call t_startf ('gmean_fixed_repro')
107 370944 : call gmean_fixed_repro(arr, arr_gmean, rel_diff, nflds)
108 370944 : call t_stopf ('gmean_fixed_repro')
109 :
110 : ! check that "fast" reproducible sum is accurate enough. If not, calculate
111 : ! using old method
112 370944 : write_warning = masterproc
113 370944 : num_err = 0
114 370944 : if (shr_reprosum_tolExceeded('gmean', nflds, write_warning, &
115 : iulog, rel_diff)) then
116 0 : if (shr_reprosum_recompute) then
117 0 : do ifld = 1, nflds
118 0 : if (rel_diff(1, ifld) > shr_reprosum_reldiffmax) then
119 0 : call gmean_float_norepro(arr(:,:,ifld), arr_gmean(ifld), ifld)
120 0 : num_err = num_err + 1
121 : end if
122 : end do
123 : end if
124 : end if
125 370944 : call t_stopf('gmean_arr')
126 370944 : if (num_err > 0) then
127 0 : call endrun('gmean: '//toString(num_err)//' reprosum errors found')
128 : end if
129 :
130 370944 : end subroutine gmean_arr
131 :
132 : !
133 : !========================================================================
134 : !
135 :
136 0 : subroutine gmean_scl (arr, gmean)
137 : use phys_grid, only: get_ncols_p
138 :
139 : !-----------------------------------------------------------------------
140 : !
141 : ! Purpose:
142 : ! Compute the global mean of each field in "arr" in the physics
143 : ! chunked decomposition
144 : !
145 : !-----------------------------------------------------------------------
146 : !
147 : ! Arguments
148 : !
149 : real(r8), intent(in) :: arr(pcols, begchunk:endchunk)
150 : ! Input array, chunked
151 : real(r8), intent(out):: gmean ! global means
152 : !
153 : ! Local workspace
154 : !
155 : integer, parameter :: nflds = 1
156 : real(r8) :: gmean_array(nflds)
157 0 : real(r8) :: array(pcols, begchunk:endchunk, nflds)
158 : integer :: ncols, lchnk
159 :
160 0 : do lchnk = begchunk, endchunk
161 0 : ncols = get_ncols_p(lchnk)
162 0 : array(:ncols, lchnk, 1) = arr(:ncols, lchnk)
163 : end do
164 0 : call gmean_arr(array, gmean_array, nflds)
165 0 : gmean = gmean_array(1)
166 :
167 0 : end subroutine gmean_scl
168 :
169 : !
170 : !========================================================================
171 : !
172 :
173 0 : subroutine gmean_float_norepro(arr, repro_sum, index)
174 : !-----------------------------------------------------------------------
175 : !
176 : ! Purpose:
177 : ! Compute the global mean of <arr> in the physics chunked
178 : ! decomposition using a fast but non-reproducible algorithm.
179 : ! Log that value along with the value computed by
180 : ! shr_reprosum_calc (<repro_sum>)
181 : !
182 : !-----------------------------------------------------------------------
183 :
184 0 : use physconst, only: pi
185 : use spmd_utils, only: masterproc, masterprocid, MPI_REAL8, MPI_SUM, mpicom
186 : use phys_grid, only: get_ncols_p, get_wght_p
187 : !
188 : ! Arguments
189 : !
190 : real(r8), intent(in) :: arr(pcols, begchunk:endchunk)
191 : real(r8), intent(in) :: repro_sum ! Value computed by reprosum
192 : integer, intent(in) :: index ! Index of field in original call
193 : !
194 : ! Local workspace
195 : !
196 : integer :: lchnk, ncols, icol
197 : integer :: ierr
198 : real(r8) :: wght
199 : real(r8) :: check
200 : real(r8) :: check_sum
201 : real(r8), parameter :: pi4 = 4.0_r8 * pi
202 :
203 : !
204 : !-----------------------------------------------------------------------
205 : !
206 : ! Calculate and print out non-reproducible value
207 0 : check = 0.0_r8
208 0 : do lchnk = begchunk, endchunk
209 0 : ncols = get_ncols_p(lchnk)
210 0 : do icol = 1, ncols
211 0 : wght = get_wght_p(lchnk, icol)
212 0 : check = check + arr(icol, lchnk) * wght
213 : end do
214 : end do
215 : call MPI_reduce(check, check_sum, 1, MPI_REAL8, check_sum, MPI_SUM, &
216 0 : masterprocid, mpicom, ierr)
217 0 : if (masterproc) then
218 0 : write(iulog, '(a,i0,2(a,e20.13e2))') 'gmean(', index, ') = ', &
219 0 : check_sum / pi4, ', reprosum reported ', repro_sum
220 : end if
221 :
222 0 : end subroutine gmean_float_norepro
223 :
224 : !
225 : !========================================================================
226 : !
227 741888 : subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds)
228 : !-----------------------------------------------------------------------
229 : !
230 : ! Purpose:
231 : ! Compute the global mean of each field in "arr" in the physics
232 : ! chunked decomposition with a reproducible yet scalable implementation
233 : ! based on a fixed-point algorithm.
234 : !
235 : !-----------------------------------------------------------------------
236 0 : use spmd_utils, only: mpicom
237 : use phys_grid, only: get_ncols_p, get_wght_all_p, get_nlcols_p
238 : use phys_grid, only: ngcols_p => num_global_phys_cols
239 : use physconst, only: pi
240 : !
241 : ! Arguments
242 : !
243 : integer, intent(in) :: nflds ! number of fields
244 : real(r8), intent(in) :: arr(pcols,begchunk:endchunk,nflds)
245 : ! arr_gmean: output global sums
246 : real(r8), intent(out) :: arr_gmean(nflds)
247 : ! rel_diff: relative and absolute differences from shr_reprosum_calc
248 : real(r8), intent(out) :: rel_diff(2, nflds)
249 : !
250 : ! Local workspace
251 : !
252 : integer :: lchnk, icol, ifld ! chunk, column, field indices
253 : integer :: ncols ! # columns in current chunk
254 : integer :: count ! summand count
255 :
256 : real(r8) :: wght(pcols) ! integration weights
257 370944 : real(r8), allocatable :: xfld(:,:) ! weighted summands
258 : integer :: nlcols
259 : !
260 : !-----------------------------------------------------------------------
261 : !
262 370944 : nlcols = get_nlcols_p()
263 1483776 : allocate(xfld(nlcols, nflds))
264 :
265 : ! pre-weight summands
266 1854720 : do ifld = 1, nflds
267 1483776 : count = 0
268 7836192 : do lchnk = begchunk, endchunk
269 5981472 : ncols = get_ncols_p(lchnk)
270 5981472 : call get_wght_all_p(lchnk, ncols, wght)
271 101360448 : do icol = 1, ncols
272 93895200 : count = count + 1
273 99876672 : xfld(count, ifld) = arr(icol, lchnk, ifld) * wght(icol)
274 : end do
275 : end do
276 : end do
277 :
278 : ! call fixed-point algorithm
279 : call shr_reprosum_calc (xfld, arr_gmean, count, nlcols, nflds, &
280 370944 : gbl_count=ngcols_p, commid=mpicom, rel_diff=rel_diff)
281 :
282 370944 : deallocate(xfld)
283 : ! final normalization
284 1854720 : arr_gmean(:) = arr_gmean(:) / (4.0_r8 * pi)
285 :
286 370944 : end subroutine gmean_fixed_repro
287 :
288 0 : subroutine test_gmean(max_diff)
289 : ! Test gmean on some different field patterns
290 : ! Test 1: Just 1, easy peasy
291 : ! Test 2: Positive definite, moderate dynamic range
292 : ! Test 3: Positive definite, large dynamic range (pattern 1)
293 : ! Test 4: Positive definite, large dynamic range (pattern 2)
294 : ! Test 5: Large dynamic range (pattern 1)
295 : ! Test 6: Large dynamic range (pattern 2)
296 370944 : use shr_kind_mod, only: SHR_KIND_CL, INT64 => SHR_KIND_I8
297 : use physconst, only: pi
298 : use spmd_utils, only: iam, masterproc
299 : use cam_abortutils, only: endrun
300 : use cam_logfile, only: iulog
301 : use phys_grid, only: get_ncols_p, get_gcol_p, get_wght_p
302 : use phys_grid, only: ngcols_p => num_global_phys_cols
303 :
304 : ! Dummy argument
305 : real(r8), optional, intent(in) :: max_diff
306 : ! Local variables
307 : integer, parameter :: num_tests = 6
308 : integer :: lchnk, ncols, icol, gcol, findex
309 : integer(INT64) :: test_val
310 0 : real(r8) :: test_arr(pcols,begchunk:endchunk,num_tests)
311 : real(r8) :: test_mean(num_tests)
312 : real(r8) :: expect(num_tests)
313 : real(r8) :: diff, wght
314 : real(r8) :: max_diff_use
315 : real(r8), parameter :: fact2 = 1.0e-8_r8
316 : real(r8), parameter :: ifact = 1.0e6_r8
317 : real(r8), parameter :: pi4 = 4.0_r8 * pi
318 : real(r8), parameter :: max_diff_def = 1.0e-14_r8
319 : character(len=SHR_KIND_CL) :: errmsg(num_tests)
320 : character(len=*), parameter :: subname = 'test_gmean: '
321 :
322 0 : if (present(max_diff)) then
323 0 : max_diff_use = max_diff
324 : else
325 : max_diff_use = max_diff_def
326 : end if
327 0 : test_arr = 0.0_r8
328 0 : do lchnk = begchunk, endchunk
329 0 : ncols = get_ncols_p(lchnk)
330 0 : do icol = 1, ncols
331 0 : gcol = get_gcol_p(lchnk, icol)
332 0 : test_arr(icol, lchnk, 1) = 1.0_r8
333 0 : wght = get_wght_p(lchnk, icol)
334 0 : test_arr(icol, lchnk, 2) = real(gcol, r8) * pi4 / wght
335 0 : test_arr(icol, lchnk, 3) = test_arr(icol, lchnk, 2) * fact2
336 0 : if (mod(gcol, 2) == 1) then
337 0 : test_arr(icol, lchnk, 4) = test_arr(icol, lchnk, 3) + ifact
338 0 : test_arr(icol, lchnk, 6) = test_arr(icol, lchnk, 3) + ifact
339 : else
340 0 : test_arr(icol, lchnk, 4) = test_arr(icol, lchnk, 3)
341 0 : test_arr(icol, lchnk, 6) = test_arr(icol, lchnk, 3) - ifact
342 : end if
343 0 : if (gcol > (ngcols_p / 2)) then
344 0 : test_arr(icol, lchnk, 5) = test_arr(icol, lchnk, 3) + ifact
345 0 : test_arr(icol, lchnk, 3) = test_arr(icol, lchnk, 3) + ifact
346 : else
347 : ! test_arr 3 already has correct value
348 0 : test_arr(icol, lchnk, 5) = test_arr(icol, lchnk, 3) - ifact
349 : end if
350 : end do
351 : end do
352 0 : test_mean(:) = -2.71828_r8 * pi
353 0 : expect(1) = 1.0_r8
354 0 : test_val = int(ngcols_p, INT64)
355 0 : test_val = (test_val + 1) * test_val / 2_INT64
356 0 : expect(2) = real(test_val, r8)
357 0 : expect(3) = (expect(2) * fact2) + (ifact / 2.0_r8)
358 0 : expect(4) = expect(3)
359 0 : expect(5) = expect(2) * fact2
360 0 : expect(6) = expect(5)
361 0 : call gmean(test_arr, test_mean, num_tests)
362 0 : errmsg = ''
363 0 : do findex = 1, num_tests
364 0 : diff = abs(test_mean(findex) - expect(findex)) / expect(findex)
365 0 : if (diff > max_diff_use) then
366 0 : write(errmsg(findex), '(i0,a,i0,3(a,e20.13e2))') iam, &
367 0 : ': test_mean(', findex, ') FAIL: ', test_mean(findex), &
368 0 : ' /= ', expect(findex), ', diff = ', diff
369 : end if
370 : end do
371 0 : if (ANY(len_trim(errmsg) > 0)) then
372 : call endrun(subname//trim(errmsg(1))//'\n'//trim(errmsg(2))//'\n'// &
373 : trim(errmsg(3))//'\n'//trim(errmsg(4))//'\n'// &
374 0 : trim(errmsg(5))//'\n'//trim(errmsg(6)))
375 : end if
376 0 : if (masterproc) then
377 0 : do findex = 1, num_tests
378 0 : write(iulog, '(2a,i0,a,e20.13e2)') subname, 'test_mean(', findex, &
379 0 : ') = ', test_mean(findex)
380 : end do
381 : end if
382 0 : end subroutine test_gmean
383 :
384 : !
385 : !========================================================================
386 : !
387 :
388 : end module gmean_mod
|
