Test_RegridderManager.pf Source File

#define _VERIFY(status) \
   if(status /= 0) then; \
      call assert_that(status, is(0), location=SourceLocation(__FILE__,__LINE__)); \
      if (anyExceptions()) return; \
   endif
#define _RC rc=status); _VERIFY(status

! Helper procedures
#define _SUCCESS 0
#define _RC2 rc=status); _VERIFY2(status
#define _VERIFY2(status)  if (status /= 0) then; if (present(rc)) rc=status; return; endif
#define _RETURN(status) if (present(rc)) rc=status; return

module Test_RegridderManager
   use pfunit
   use mapl3g_regridder_mgr
   use mapl3g_geom_mgr
   use mapl_BaseMod, only: MAPL_UNDEF
   use esmf_TestMethod_mod ! mapl
   use esmf
   implicit none

contains

   ! Helper procedures
   ! TODO add error handling to helper procedures
   
   function make_geom(geom_mgr, hconfig, rc) result(geom)
      type(ESMF_Geom) :: geom
      type(GeomManager), intent(inout) :: geom_mgr
      type(ESMF_HConfig), optional, intent(in) :: hconfig
      integer, optional, intent(out) :: rc

      type(MaplGeom), pointer :: mapl_geom
      integer :: status
      type(ESMF_HConfig) :: hconfig_

      hconfig_ = ESMF_HConfigCreate(content="{class: latlon, im_world: 12, jm_world: 13, pole: PC, dateline: DC, nx: 1, ny: 1}", &
           _RC2)
      if (present(hconfig)) hconfig_ = hconfig

      mapl_geom => geom_mgr%get_mapl_geom(hconfig_, _RC2)
      geom = mapl_geom%get_geom()

      _RETURN(_SUCCESS)
   end function make_geom

   function make_field(geom, name, value, lm, rc) result(field)
      type(ESMF_Field) :: field
      type(ESMF_Geom), intent(in) :: geom
      character(*), intent(in) :: name
      real(kind=ESMF_KIND_R4), intent(in) :: value
      integer, optional, intent(in) :: lm
      integer, optional, intent(out) :: rc

      real(kind=ESMF_KIND_R4), pointer :: x(:,:)
      real(kind=ESMF_KIND_R4), pointer :: x_3d(:,:,:)
      integer :: status

      field = ESMF_FieldEmptyCreate(name=name, _RC2)
      call ESMF_FieldEmptySet(field, geom, _RC2)
      if (present(lm)) then
         call ESMF_FieldEmptyComplete(field, typekind=ESMF_TYPEKIND_R4, ungriddedLBound=[1], ungriddedUBound=[lm], _RC2)
         call ESMF_FieldGet(field, farrayptr=x_3d,_RC2)
         x_3d = value
      else
         call ESMF_FieldEmptyComplete(field, typekind=ESMF_TYPEKIND_R4, _RC2)
         call ESMF_FieldGet(field, farrayptr=x, _RC2)
         x = value
      end if

      _RETURN(_SUCCESS)
   end function make_field

   @test(type=ESMF_TestMethod, npes=[1])
   ! Just execute a series of plausible commands and ensure that no
   ! failures are indicated Regrid a constant field onto identical
   ! geometry should not change any values.
   subroutine test_basic(this)
      class(ESMF_TestMethod), intent(inout) :: this
      type(GeomManager) :: geom_mgr
      type(RegridderManager) :: regridder_mgr
      type(RegridderSpec) :: regridder_spec
      integer :: status
      class(Regridder), pointer :: my_regridder
      type(ESMF_Geom) :: geom

      type(ESMF_Field) :: f1, f2
      real(kind=ESMF_KIND_R4), pointer :: x(:,:)


      geom_mgr = GeomManager()
      regridder_mgr = RegridderManager()

      geom = make_geom(geom_mgr, _RC)

      ! use default esmf regrid parameters: method, zero region, etc
      regridder_spec = RegridderSpec(EsmfRegridderParam(), geom, geom)

      my_regridder => regridder_mgr%get_regridder(regridder_spec, _RC)

      f1 = make_field(geom, 'f1', value=3._ESMF_KIND_R4, _RC)
      f2 = make_field(geom, 'f2', value=0._ESMF_KIND_R4, _RC)

      call my_regridder%regrid(f1, f2, _RC)
      call ESMF_FieldGet(f2, farrayptr=x, _RC)

      @assert_that(x, every_item(is(equal_to(3._ESMF_KIND_R4))))

   end subroutine test_basic

   @test(type=ESMF_TestMethod, npes=[1])
   ! Test that identical spec returns same regridder object.  I.e.,
   ! that the manager is properly caching.
   subroutine test_reuse_regridder(this)
      class(ESMF_TestMethod), intent(inout) :: this
      type(GeomManager) :: geom_mgr
      type(RegridderManager), target :: regridder_mgr
      type(RegridderSpec) :: regridder_spec
      integer :: status
      class(Regridder), pointer :: regridder_1, regridder_2
      type(ESMF_Geom) :: geom

      geom_mgr = GeomManager()
      regridder_mgr = RegridderManager()

      geom = make_geom(geom_mgr, _RC)

      regridder_spec = RegridderSpec(EsmfRegridderParam(), geom, geom)

      regridder_1 => regridder_mgr%get_regridder(regridder_spec, _RC)

      regridder_2 => regridder_mgr%get_regridder(regridder_spec, _RC)

      @assertTrue(associated(regridder_2, regridder_1))
   end subroutine test_reuse_regridder

   @test(type=ESMF_TestMethod, npes=[1])
   ! Test that different spec returns different regridder object.  I.e.,
   ! that the manager is properly caching.
   subroutine test_do_not_reuse_regridder(this)
      class(ESMF_TestMethod), intent(inout) :: this
      type(GeomManager) :: geom_mgr
      type(RegridderManager), target :: regridder_mgr
      type(RegridderSpec) :: spec_1, spec_2
      integer :: status
      class(Regridder), pointer :: regridder_1, regridder_2
      type(ESMF_Geom) :: geom_1, geom_2
      type(ESMF_HConfig) :: hconfig

      geom_mgr = GeomManager()
      regridder_mgr = RegridderManager()

      geom_1 = make_geom(geom_mgr, _RC)

      hconfig = ESMF_HConfigCreate(content="{class: latlon, im_world: 12, jm_world: 11, pole: PC, dateline: DC, nx: 1, ny: 1}", _RC)
      geom_2 = make_geom(geom_mgr, hconfig, _RC) ! variant of geom_1

      
      spec_1 = RegridderSpec(EsmfRegridderParam(), geom_1, geom_1)
      regridder_1 => regridder_mgr%get_regridder(spec_1, _RC)

      spec_2 = RegridderSpec(EsmfRegridderParam(), geom_1, geom_2)
      regridder_2 => regridder_mgr%get_regridder(spec_2, _RC)

      @assertFalse(associated(regridder_1, regridder_2))
   end subroutine test_do_not_reuse_regridder

   @test(type=ESMF_TestMethod, npes=[1])
   ! Test realistic regridding.  A checkerboard input field (in
   ! longitude) with constant spacing should produce a constant output
   ! grid with default bilinear regrid method.
   subroutine test_regrid_values(this)
      class(ESMF_TestMethod), intent(inout) :: this
      type(GeomManager) :: geom_mgr
      type(RegridderManager), target :: regridder_mgr
      type(RegridderSpec) :: spec
      integer :: status
      class(Regridder), pointer :: my_regridder
      type(ESMF_Geom) :: geom_1, geom_2
      type(ESMF_HConfig) :: hconfig
      type(ESMF_Field) :: f1, f2
      real(kind=ESMF_KIND_R4), pointer :: x1(:,:)
      real(kind=ESMF_KIND_R4), pointer :: x2(:,:)

      geom_mgr = GeomManager()
      regridder_mgr = RegridderManager()

      hconfig = ESMF_HConfigCreate(content="{class: latlon, im_world: 12, jm_world: 11, pole: PC, dateline: DE, nx: 1, ny: 1}", _RC)
      geom_1 = make_geom(geom_mgr, hconfig, _RC)

      hconfig = ESMF_HConfigCreate(content="{class: latlon, im_world: 6, jm_world: 5, pole: PC, dateline: DE, nx: 1, ny: 1}", _RC)
      geom_2 = make_geom(geom_mgr, hconfig, _RC) ! variant of geom_1

      spec = RegridderSpec(EsmfRegridderParam(regridmethod=ESMF_REGRIDMETHOD_CONSERVE), geom_1, geom_2)
      my_regridder => regridder_mgr%get_regridder(spec, _RC)

      f1 = make_field(geom_1, 'f1', value=2._ESMF_KIND_R4, _RC)
      call ESMF_FieldGet(f1, farrayptr=x1, _RC)
      x1(2::2,:) = 0 ! checkerboard
      
      f2 = make_field(geom_2, 'f2', value=0._ESMF_KIND_R4, _RC)

      ! (0 + 2)/2 == 1
      call my_regridder%regrid(f1, f2, _RC)
      call ESMF_FieldGet(f2, farrayptr=x2, _RC)
 
      @assert_that(x2, every_item(is(equal_to(1._ESMF_KIND_R4))))

   end subroutine test_regrid_values


   @test(type=ESMF_TestMethod, npes=[1])
   ! Test regridding on fields with ungridded dimensions.  ESMF does
   ! not directly support this case, and this test is to drive the
   ! creation of a wrapper layer in MAPL.
   subroutine test_regrid_3d(this)
      class(ESMF_TestMethod), intent(inout) :: this
      type(GeomManager) :: geom_mgr
      type(RegridderManager), target :: regridder_mgr
      type(RegridderSpec) :: spec
      integer :: status
      class(Regridder), pointer :: my_regridder
      type(ESMF_Geom) :: geom_1, geom_2
      type(ESMF_HConfig) :: hconfig
      type(ESMF_Field) :: f1, f2
      real(kind=ESMF_KIND_R4), pointer :: x1(:,:,:)
      real(kind=ESMF_KIND_R4), pointer :: x2(:,:,:)

      type(DynamicMask) :: dyn_mask
      
      geom_mgr = GeomManager()
      regridder_mgr = RegridderManager()

      hconfig = ESMF_HConfigCreate(content="{class: latlon, im_world: 12, jm_world: 11, pole: PC, dateline: DE, nx: 1, ny: 1}", _RC)
      geom_1 = make_geom(geom_mgr, hconfig, _RC)

      hconfig = ESMF_HConfigCreate(content="{class: latlon, im_world: 6, jm_world: 5, pole: PC, dateline: DE, nx: 1, ny: 1}", _RC)
      geom_2 = make_geom(geom_mgr, hconfig, _RC) ! variant of geom_1

      dyn_mask = DynamicMask(mask_type='missing_value', src_mask_value=real(MAPL_UNDEF,kind=ESMF_KIND_R8), handleAllElements=.true.,_RC)
      
      spec = RegridderSpec(EsmfRegridderParam(regridmethod=ESMF_REGRIDMETHOD_CONSERVE, dyn_mask=dyn_mask), geom_1, geom_2)
      my_regridder => regridder_mgr%get_regridder(spec, _RC)

      f1 = make_field(geom_1, 'f1', value=2._ESMF_KIND_R4, lm=2, _RC)
      call ESMF_FieldGet(f1, farrayptr=x1)
      x1(::4,6,1) = MAPL_UNDEF ! missing bits in level 1
      x1(1::2,:,2) = 0 ! checkerboard on level 2
      
      f2 = make_field(geom_2, 'f2', value=0._ESMF_KIND_R4, lm=2, _RC)

      call my_regridder%regrid(f1, f2, _RC)

      call ESMF_FieldGet(f2, farrayptr=x2, _RC)


      ! Missing elements case
      @assert_that(x2(1:2,:,1), every_item(is(equal_to(2._ESMF_KIND_R4))))
      ! Non missing elements case

      ! Weirdly this introduces roundoff that was not present in the
      ! previous test.  This has been reported to the ESMF core team.
      @assert_that(x2(:,:,2), every_item(is(near(1._ESMF_KIND_R4, 1.e-6))))

   end subroutine test_regrid_3d


end module Test_RegridderManager