LCOV - code coverage report
Current view: | top level - utils - gmean_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 38 | 110 | 34.5 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 3 | 6 | 50.0 %

              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     2233344 :    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     4466688 :       real(r8)                   :: rel_diff(2, nflds)
          99             :       integer                    :: ifld ! field index
         100             :       integer                    :: num_err
         101             :       logical                    :: write_warning
         102             :       !
         103             :       !-----------------------------------------------------------------------
         104             :       !
         105     2233344 :       call t_startf('gmean_arr')
         106     2233344 :       call t_startf ('gmean_fixed_repro')
         107     2233344 :       call gmean_fixed_repro(arr, arr_gmean, rel_diff, nflds)
         108     2233344 :       call t_stopf ('gmean_fixed_repro')
         109             :
         110             :       ! check that "fast" reproducible sum is accurate enough. If not, calculate
         111             :       ! using old method
         112     2233344 :       write_warning = masterproc
         113     2233344 :       num_err = 0
         114     2233344 :       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     2233344 :       call t_stopf('gmean_arr')
         126     2233344 :       if (num_err > 0) then
         127           0 :          call endrun('gmean: '//toString(num_err)//' reprosum errors found')
         128             :       end if
         129             :
         130     2233344 :    end subroutine gmean_arr
         131             :
         132             :    !
         133             :    !========================================================================
         134             :    !
         135             :
         136     3724800 :    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     3724800 :       real(r8)           :: array(pcols, begchunk:endchunk, nflds)
         158             :       integer            :: ncols, lchnk
         159             :
         160     9370200 :       do lchnk = begchunk, endchunk
         161     7507800 :          ncols = get_ncols_p(lchnk)
         162   127225200 :          array(:ncols, lchnk, 1) = arr(:ncols, lchnk)
         163             :       end do
         164     1862400 :       call gmean_arr(array, gmean_array, nflds)
         165     1862400 :       gmean = gmean_array(1)
         166             :
         167     1862400 :    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     1862400 :       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     4466688 :    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     2233344 :       real(r8), allocatable :: xfld(:,:)         ! weighted summands
         258             :       integer               :: nlcols
         259             :       !
         260             :       !-----------------------------------------------------------------------
         261             :       !
         262     2233344 :       nlcols = get_nlcols_p()
         263     8933376 :       allocate(xfld(nlcols, nflds))
         264             :
         265             :       ! pre-weight summands
         266     5579520 :       do ifld = 1, nflds
         267     3346176 :          count = 0
         268    19068792 :          do lchnk = begchunk, endchunk
         269    13489272 :             ncols = get_ncols_p(lchnk)
         270    13489272 :             call get_wght_all_p(lchnk, ncols, wght)
         271   228585648 :             do icol = 1, ncols
         272   211750200 :                count = count + 1
         273   225239472 :                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     2233344 :            gbl_count=ngcols_p, commid=mpicom, rel_diff=rel_diff)
         281             :
         282     2233344 :       deallocate(xfld)
         283             :       ! final normalization
         284     5579520 :       arr_gmean(:) = arr_gmean(:) / (4.0_r8 * pi)
         285             :
         286     2233344 :    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     2233344 :       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

Generated by: LCOV version 1.14