diff --git a/config_src/drivers/unit_tests/test_MOM_array_transform.F90 b/config_src/drivers/unit_tests/test_MOM_array_transform.F90 new file mode 100644 index 0000000000..e700791ce3 --- /dev/null +++ b/config_src/drivers/unit_tests/test_MOM_array_transform.F90 @@ -0,0 +1,21 @@ +! This file is part of MOM6, the Modular Ocean Model version 6. +! See the LICENSE file for licensing information. +! SPDX-License-Identifier: Apache-2.0 + +program test_MOM_array_transform + +use MOM_array_transform, only : symmetric_sum_unit_tests +use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, set_skip_mpi + +character(len=120) :: fail_message + +call set_skip_mpi(.true.) ! This unit tests is not expecting MPI to be used + +if ( symmetric_sum_unit_tests(fail_message) ) then + call MOM_error(FATAL, "symmetric_sum_unit_tests FAILED: "//trim(fail_message)) + stop 1 +else + call MOM_mesg("symmetric_sum_unit_tests PASSED.") +endif + +end program test_MOM_array_transform diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90 index 33218f4e34..cae23fbb91 100644 --- a/src/core/MOM_unit_tests.F90 +++ b/src/core/MOM_unit_tests.F90 @@ -7,6 +7,7 @@ module MOM_unit_tests ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_array_transform, only : symmetric_sum_unit_tests use MOM_diag_buffers, only : diag_buffer_unit_tests_2d, diag_buffer_unit_tests_3d use MOM_error_handler, only : MOM_error, FATAL, is_root_pe use MOM_hor_bnd_diffusion, only : near_boundary_unit_tests @@ -31,6 +32,7 @@ subroutine unit_tests(verbosity) ! Arguments integer, intent(in) :: verbosity !< The verbosity level ! Local variables + character(len=120) :: fail_message logical :: verbose verbose = verbosity>=5 @@ -38,6 +40,8 @@ subroutine unit_tests(verbosity) if (is_root_pe()) then ! The following need only be tested on 1 PE if (string_functions_unit_tests(verbose)) call MOM_error(FATAL, & "MOM_unit_tests: string_functions_unit_tests FAILED") + if (symmetric_sum_unit_tests(fail_message)) call MOM_error(FATAL, & + "MOM_unit_tests: symmetric_sum_unit_tests FAILED: "//trim(fail_message)) if (EOS_unit_tests(verbose)) call MOM_error(FATAL, & "MOM_unit_tests: EOS_unit_tests FAILED") if (remapping_unit_tests(verbose)) call MOM_error(FATAL, & diff --git a/src/framework/MOM_array_transform.F90 b/src/framework/MOM_array_transform.F90 index 30aeadcdf1..8e8418257c 100644 --- a/src/framework/MOM_array_transform.F90 +++ b/src/framework/MOM_array_transform.F90 @@ -14,15 +14,22 @@ !! !! 90 degree rotations change the shape of the field, and are handled !! separately from 180 degree rotations. +!! +!! It also provides the symmetric_sum functions to do a rotationally invariant +!! sum of the contents of a 1d or 2d array. module MOM_array_transform +use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit + implicit none ; private public rotate_array public rotate_array_pair public rotate_vector public allocate_rotated_array +public symmetric_sum +public symmetric_sum_unit_tests !> Rotate the elements of an array to the rotated set of indices. @@ -62,8 +69,7 @@ module MOM_array_transform end interface rotate_vector -!> Allocate an array based on the rotated index map of an unrotated reference -!! array. +!> Allocate an array based on the rotated index map of an unrotated reference array. interface allocate_rotated_array module procedure allocate_rotated_array_real_2d module procedure allocate_rotated_array_real_3d @@ -71,6 +77,13 @@ module MOM_array_transform module procedure allocate_rotated_array_integer end interface allocate_rotated_array + +!> Return a rotationally symmetric sum of the elements an array. +interface symmetric_sum + module procedure symmetric_sum_1d, symmetric_sum_2d +end interface symmetric_sum + + contains !> Rotate the elements of a 2d real array along first and second axes. @@ -358,4 +371,197 @@ subroutine allocate_rotated_array_integer(A_in, lb, turns, A) endif end subroutine allocate_rotated_array_integer + +!> Do a rotationally symmetric sum of a 1-d array +function symmetric_sum_1d(field) result(sum) + real, dimension(1:), intent(in) :: field !< The field to sum in arbitrary units [A ~> a] + real :: sum !< The rotationally symmetric sum of the entries in field [A ~> a] + + ! Local variables + integer :: i, szi, szi_2 + + szi = size(field, 1) + szi_2 = szi / 2 ! Note that for an odd number sz_2 is rounded down. + sum = 0.0 + if (2*szi_2 < szi) sum = field(szi_2+1) + ! Add pairs of values, working from the inside out. + do i=szi_2,1,-1 + sum = sum + (field(i) + field(szi+1-i)) + enddo +end function symmetric_sum_1d + + +!> Do a rotationally symmetric sum of a 2-d array using a recursive "Union-Jack" pattern of addition. +recursive function symmetric_sum_2d(field) result(sum) + real, dimension(1:,1:), intent(in) :: field !< The field to sum in arbitrary units [A ~> a] + real :: sum !< The rotationally symmetric sum of the entries in field [A ~> a] + + ! Local variables + real :: quad_sum(2,2) ! The sums in each of the quadrants [A ~> a] + logical :: odd_i, odd_j + integer :: ij, szi, szj, szi_2, szj_2, ic, jc + + szi = size(field, 1) ; szj = size(field, 2) + ! These 5 special cases are equivalent to the general case, but they reduce the use + ! of complicated logic for common simple cases. + if ((szi == 1) .and. (szj == 1)) then + sum = field(1,1) + elseif ((szi == 2) .and. (szj == 2)) then + sum = (field(1,1) + field(2,2)) + (field(2,1) + field(1,2)) + elseif ((szi == 3) .and. (szj == 3)) then + sum = (field(2,2) + ((field(1,2) + field(3,2)) + (field(2,1) + field(2,3)))) + & + ((field(1,1) + field(3,3)) + (field(3,1) + field(1,3))) + elseif (szi == 1) then + sum = symmetric_sum_1d(field(1,:)) + elseif (szj == 1) then + sum = symmetric_sum_1d(field(:,1)) + else + ! This is the general case. + ! Note that for an odd number szi_2 is rounded down. + szi_2 = szi / 2 + szj_2 = szj / 2 + + odd_i = (2*szi_2 < szi) ! This could be (modulo(szi,2) == 1) + odd_j = (2*szj_2 < szj) + ! Start by finding the sums along the central axes if there are an odd number of points. + if (odd_i .and. odd_j) then + ic = szi_2+1 ; jc = szj_2+1 ! The index of the central point + sum = field(ic,jc) + ! Add pairs of pairs of values, working from the inside out. + do ij=1,min(szi_2,szj_2) + sum = sum + ((field(ic-ij,jc) + field(ic+ij,jc)) + (field(ic,jc-ij) + field(ic,jc+ij))) + enddo + ! Add extra pairs of values, working from the inside out. + if (szi_2 > szj_2) then + do ij=szj_2+1,szi_2 + sum = sum + (field(ic-ij,jc) + field(ic+ij,jc)) + enddo + elseif (szj_2 > szi_2) then + do ij=szi_2+1,szj_2 + sum = sum + (field(ic,jc-ij) + field(ic,jc+ij)) + enddo + endif + elseif (odd_i) then + sum = symmetric_sum_1d(field(szi_2+1,1:szj)) + elseif (odd_j) then + sum = symmetric_sum_1d(field(1:szi,szj_2+1)) + else + sum = 0.0 + endif + + ! Find the sums in the four quadrants of the array. + if ((szi_2 > 1) .and. (szj_2 > 1)) then + ! Use a recursive call to symmetric_sum_2d to determine the sums in the corner quadrants. + quad_sum(1,1) = symmetric_sum_2d(field(1:szi_2,1:szj_2)) + quad_sum(2,1) = symmetric_sum_2d(field(szi+1-szi_2:szi,1:szj_2)) + quad_sum(1,2) = symmetric_sum_2d(field(1:szi_2,szj+1-szj_2:szj)) + quad_sum(2,2) = symmetric_sum_2d(field(szi+1-szi_2:szi,szj+1-szj_2:szj)) + elseif (szi_2 > 1) then + quad_sum(1,1) = symmetric_sum_1d(field(1:szi_2,1)) + quad_sum(2,1) = symmetric_sum_1d(field(szi+1-szi_2:szi,1)) + quad_sum(1,2) = symmetric_sum_1d(field(1:szi_2,szj)) + quad_sum(2,2) = symmetric_sum_1d(field(szi+1-szi_2:szi,szj)) + elseif (szj_2 > 1) then + quad_sum(1,1) = symmetric_sum_1d(field(1,1:szj_2)) + quad_sum(2,1) = symmetric_sum_1d(field(szi,1:szj_2)) + quad_sum(1,2) = symmetric_sum_1d(field(1,szj+1-szj_2:szj)) + quad_sum(2,2) = symmetric_sum_1d(field(szi,szj+1-szj_2:szj)) + else + quad_sum(1,1) = field(1,1) + quad_sum(2,1) = field(szi,1) + quad_sum(1,2) = field(1,szj) + quad_sum(2,2) = field(szi,szj) + endif + + sum = sum + ((quad_sum(1,1) + quad_sum(2,2)) + (quad_sum(2,1) + quad_sum(1,2))) + endif +end function symmetric_sum_2d + + +!> Do a naive non-rotationally symmetric sum of a 2-d array. This function is only here for testing. +function naive_sum_2d(field, abs_val) result(sum) + real, dimension(1:,1:), intent(in) :: field !< The field to sum in arbitrary units [A ~> a] + logical, optional, intent(in) :: abs_val !< If present and true, sum the absolute values + real :: sum !< The rotation dependent sum of the entries in field [A ~> a] + + ! Local variables + logical :: sum_abs_val + integer :: i, j, szi, szj + + szi = size(field, 1) ; szj = size(field, 2) + sum_abs_val = .false. ; if (present(abs_val)) sum_abs_val = abs_val + sum = 0.0 + if (sum_abs_val) then + do j=1,szj ; do i=1,szi + sum = sum + abs(field(i,j)) + enddo ; enddo + else + do j=1,szj ; do i=1,szi + sum = sum + field(i,j) + enddo ; enddo + endif +end function naive_sum_2d + + +!> Returns true if a unit test of the symmetric sums fails. +logical function symmetric_sum_unit_tests(fail_message) + ! Arguments + character(len=120), intent(out) :: fail_message !< Blank or a description of the first failed test. + ! Local variables + integer, parameter :: sz=13 ! The maximum size of the test arrays + real :: array(sz,sz) ! An array of inexact real values for testing in arbitrary units [A] + real :: ar_90(sz,sz) ! Array rotated by 90 degrees in arbitrary units [A] + real :: ar_180(sz,sz) ! Array rotated by 180 degrees in arbitrary units [A] + real :: ar_270(sz,sz) ! Array rotated by 270 degrees in arbitrary units [A] + real :: sum(5) ! Different versions of sums over a sub-array [A] + real :: abs_sum ! The sum of the absolute values of the array [A] + real :: tol ! The tolerance for an inexact test [A] + + character(len=120) :: mesg + integer :: i, j, n, m, r + logical :: fail + + fail = .false. + fail_message = "" + + ! Fill the array with real numbers that can not be represented exactly. + do j=1,sz ; do i=1,sz + array(i,j) = 1.0 / (2.0*(j*sz + i) + 1.0) + ! Combining positive and negative numbers amplifies differences from the order of arithmetic. + if (modulo(i+j, 2) == 0) array(i,j) = -array(i,j) + enddo ; enddo + call rotate_array_real_2d(array, 1, ar_90) + call rotate_array_real_2d(array, 2, ar_180) + call rotate_array_real_2d(array, 3, ar_270) + + do n = 1, sz ; do m = 1, sz + sum(1) = symmetric_sum(array(1:n,1:m)) + sum(2) = symmetric_sum(ar_90(sz+1-m:sz,1:n)) + sum(3) = symmetric_sum(ar_180(sz+1-n:sz,sz+1-m:sz)) + sum(4) = symmetric_sum(ar_270(1:m,sz+1-n:sz)) + sum(5) = naive_sum_2d(array(1:n,1:m)) + abs_sum = naive_sum_2d(array(1:n,1:m), abs_val=.true.) + tol = 2.0 * abs_sum * epsilon(abs_sum) + if (abs(sum(1) - sum(5)) > tol) then + write(mesg,'(i0," x ",i0," symmetric vs naive sum, sum=",ES13.5," diff=",ES13.5)') & + n, m, sum(1), sum(5) - sum(1) + write(stdout,*) "Symmetric_sum_failure: "//trim(mesg) + write(stderr,*) "Symmetric_sum_failure: "//trim(mesg) + if (.not.fail) fail_message = mesg ! This is the first failed test. + fail = .true. + endif + do r = 2, 4 ; if (abs(sum(1) - sum(r)) > 0.0) then + write(mesg,'(i0," x ",i0," with ",i0," degree rotation, sum=",ES13.5," diff=",ES13.5)') & + n, m, 90*(r-1), sum(1), sum(r) - sum(1) + write(stdout,*) "Symmetric_sum_failure: "//trim(mesg) + write(stderr,*) "Symmetric_sum_failure: "//trim(mesg) + if (.not.fail) fail_message = mesg ! This is the first failed test. + fail = .true. + endif ; enddo + enddo ; enddo + + symmetric_sum_unit_tests = fail + +end function symmetric_sum_unit_tests + end module MOM_array_transform