Merge branch 'ambrad/homme/dp-pg' (PR #6057)

Homme(xx)/pgN: Add pgN feature to doubly periodic planar mode.

* Add pgN feature to planar doubly periodic mode.
* Add some standalone-hommexx infrastructure to support using the same unit-test
  exe in multiple tests.
* Add an init-planar capability to Hommexx unit tests.
* Add a unit test: gllfvremap_planar_ut, planar version of old gllfvremap_ut
  unit test.
* Add a BFB test for F90/C++ dycore consistency in planar mode with pg2 physics
  grid.

e3sm_developer passes on Chrysalis. Additional tests for GPU/EAMxx pass on
various machines.

[BFB] except for new tests in homme_integration

components/homme/src/share/compose/compose_cedr.cpp

@@ -239,7 +239,7 @@ void renumber_leaves (const tree::Node::Ptr& node, const Int horiz_nleaf,
 void attach_and_renumber_horizontal_trees (const tree::Node::Ptr& supnode,
                                            const tree::Node::Ptr& htree,
                                            const Int horiz_nleaf) {
-  Int level = -1, rank;
+  Int level = -1, rank = -1;
   for (Int k = 0; k < supnode->nkids; ++k) {
     auto& kid = supnode->kids[k];
     if (kid->nkids) {

components/homme/src/share/control_mod.F90

@@ -679,7 +679,7 @@ subroutine set_planar_defaults()
        Ly = 51.2D0 * 1000.0D0
        Sx = -25.6D0 * 1000.0D0
        Sy = -25.6D0 * 1000.0D0
-    else if (test_case == "planar_rising_bubble" ) then
+    else if (test_case == "planar_rising_bubble" .or. test_case == "planar_rising_bubble_pg2") then
        Lx = 2.0D0 * 10000.0D0
        Ly = 2.0D0 * 10000.0D0
        Sx = -10000.0D0
@@ -704,7 +704,7 @@ subroutine set_planar_defaults()
     endif
     endif !if lx,ly,sx,sy are not set in nl
 
-    if (test_case == "planar_rising_bubble" ) then
+    if (test_case == "planar_rising_bubble" .or. test_case == "planar_rising_bubble_pg2") then
        case_planar_bubble = .TRUE.
     end if
 

components/homme/src/share/cxx/GllFvRemapImpl.cpp

@@ -514,6 +514,7 @@ ::run_dyn_to_fv_phys (const int timeidx, const Phys1T& ps, const Phys1T& phis, c
       kv.team_barrier(); // w2, w3, w4 in use
       // T_f
       loop_ik(ttrf, tvr, [&] (int i, int k) { T(ie,i,k) = th_f(i,k)*exner_f(i,k); });
+      kv.team_barrier(); // w2, w3 in use
     }
 
     // (u,v)
@@ -647,6 +648,7 @@ run_fv_phys_to_dyn (const int timeidx, const CPhys2T& Ts, const CPhys3T& uvs,
         // modify omega. Thus, zero FM so that the third component is 0.
         loop_ik(ttrg, tvr, [&] (int i, int k) { fm_ie(2,i,k) = 0; });
       }
+      kv.team_barrier();
     }
 
     { // T

components/homme/src/share/cxx/mpi/Connectivity.cpp

@@ -240,7 +240,7 @@ void Connectivity::setup_ucon () {
   {
     /* Check that the connection counts are OK.
          Let max_elements_attached_to_node be the number of elements attached to
-       an element corner. For example, in a regular planar mesh, it's 4. for RRM
+       an element corner. For example, in a regular planar mesh, it's 4. For RRM
        it's <= 7, as established in src/share/dimensions_mod.F90.
          Let max_corner_elem be the max number of elements attached to a corner
        that are not accounted for by the two edges and the parent element. Thus,

components/homme/src/share/cxx/mpi/Connectivity.hpp

@@ -156,7 +156,7 @@ class Connectivity
   ExecViewManaged<int*>::HostMirror h_ucon_ptr;
   ExecViewManaged<int*>             d_ucon_dir_ptr;
   ExecViewManaged<int*>::HostMirror h_ucon_dir_ptr;
-  // Helper used to accumulated connections during add_connection phase. Emptied
+  // Helper used to accumulate connections during add_connection phase. Emptied
   // in finalize. l_ is local; r_ is remote.
   struct UConInfo {
     int l_lid, l_gid, r_lid, r_gid, r_pid;

components/homme/src/share/gllfvremap_mod.F90

@@ -35,6 +35,7 @@ module gllfvremap_mod
   ! nearly identical out to day 30, whereas pg1 fields differ visibly.
   !
   ! AMB 2019/07-2020/06 Initial
+  ! AMB 2023/11         Handle planar geometry
 
   use hybrid_mod, only: hybrid_t
   use kinds, only: real_kind
@@ -84,7 +85,7 @@ module gllfvremap_mod
   ! GLL remap.
   type, private :: GllFvRemap_t
      integer :: nphys, npi, check
-     logical :: have_fv_topo_file_phis, boost_pg1, check_ok
+     logical :: have_fv_topo_file_phis, boost_pg1, check_ok, is_planar
      real(kind=real_kind) :: tolfac ! for checking
      real(kind=real_kind) :: &
           ! Node or cell weights
@@ -182,7 +183,8 @@ subroutine gfr_init(par, elem, nphys, check, boost_pg1)
     use kinds, only: iulog
     use dimensions_mod, only: nlev
     use parallel_mod, only: parallel_t, abortmp
-    use quadrature_mod, only : gausslobatto, quadrature_t
+    use quadrature_mod, only: gausslobatto, quadrature_t
+    use control_mod, only: geometry, cubed_sphere_map
 
     type (parallel_t), intent(in) :: par
     type (element_t), intent(in) :: elem(:)
@@ -236,6 +238,8 @@ subroutine gfr_init(par, elem, nphys, check, boost_pg1)
     gfr%npi = max(2, nphys)
     nphys2 = nphys*nphys
 
+    gfr%is_planar = trim(geometry) == 'plane'
+
     call gfr_init_w_gg(np, gfr%w_gg)
     call gfr_init_w_gg(gfr%npi, gfr%w_sgsg)
     call gfr_init_w_ff(nphys, gfr%w_ff)
@@ -1118,6 +1122,7 @@ subroutine gfr_init_geometry(elem, gfr)
     use coordinate_systems_mod, only: cartesian3D_t, spherical_polar_t, &
          sphere_tri_area, change_coordinates
     use cube_mod, only: ref2sphere
+    use physical_constants, only: dx, dy
 
     type (element_t), intent(in) :: elem(:)
     type (GllFvRemap_t), intent(inout) :: gfr
@@ -1154,14 +1159,19 @@ subroutine gfr_init_geometry(elem, gfr)
                       end if
                       fv_corners_xyz(ai,bi) = elem(ie)%corners3D(idx)
                    else
-                      ! p_sphere is unused. fv_corners_xyz(ai,bi) contains
-                      ! the cartesian point before it's converted to lat-lon.
-                      p_sphere = ref2sphere(ae(ai), be(bi), elem(ie)%corners3D, &
-                           cubed_sphere_map, elem(ie)%corners, elem(ie)%facenum, &
-                           fv_corners_xyz(ai,bi))
-                      if (cubed_sphere_map == 0) then
-                         ! In this case, fv_corners_xyz above is not set.
-                         fv_corners_xyz(ai,bi) = change_coordinates(p_sphere)
+                      if (gfr%is_planar) then
+                         call ref2plane(elem(ie)%corners3D, ae(ai), be(bi), &
+                              fv_corners_xyz(ai,bi))
+                      else
+                         ! fv_corners_xyz(ai,bi) contains the cartesian point
+                         ! before it's converted to lat-lon.
+                         p_sphere = ref2sphere(ae(ai), be(bi), elem(ie)%corners3D, &
+                              cubed_sphere_map, elem(ie)%corners, elem(ie)%facenum, &
+                              fv_corners_xyz(ai,bi))
+                         if (cubed_sphere_map == 0) then
+                            ! In this case, fv_corners_xyz above is not set.
+                            fv_corners_xyz(ai,bi) = change_coordinates(p_sphere)
+                         end if
                       end if
                    end if
                    ctr%x = ctr%x + fv_corners_xyz(ai,bi)%x
@@ -1170,16 +1180,29 @@ subroutine gfr_init_geometry(elem, gfr)
                 end do
              end do
 
-             if (cubed_sphere_map == 2) then
+             if (gfr%is_planar) then
+                ! The cell area is (dx*dy)/nf^2. This must then be divided by
+                ! w_ff = 4/nf^2, leaving (dx*dy)/4.
+                gfr%fv_metdet(k,ie) = (dx*dy)/four
+             elseif (cubed_sphere_map == 2) then
                 call sphere_tri_area(fv_corners_xyz(1,1), fv_corners_xyz(2,1), &
                      fv_corners_xyz(2,2), spherical_area)
                 call sphere_tri_area(fv_corners_xyz(1,1), fv_corners_xyz(2,2), &
                      fv_corners_xyz(1,2), tmp)
                 spherical_area = spherical_area + tmp
                 gfr%fv_metdet(k,ie) = spherical_area/gfr%w_ff(k)
-
-                ! Center is average of 4 corner points projected to sphere.
+             end if
+             if (gfr%is_planar .or. cubed_sphere_map == 2) then
                 ctr%x = ctr%x/four; ctr%y = ctr%y/four; ctr%z = ctr%z/four
+                if (cubed_sphere_map == 2) then
+                   ! [Projection bug 2023/11]: In the initial implementation, I
+                   ! forgot to project to the sphere. Mathematically, it doesn't
+                   ! matter: change_coordinates is invariant to the norm, and
+                   ! sphere2ref does the right thing since it solves a
+                   ! least-squares problem. In finite precision, fixing this bug
+                   ! would cause eps-level changes in ref2sphere and thus cause
+                   ! all pg2 tests to be non-BFB. I'm leaving it as is.
+                end if
                 gfr%center_f(i,j,ie) = ctr
              end if
 
@@ -1199,7 +1222,7 @@ subroutine gfr_init_geometry(elem, gfr)
        end do
     end do
 
-    if (cubed_sphere_map == 0) then
+    if (cubed_sphere_map == 0 .and. .not. gfr%is_planar) then
        ! For cubed_sphere_map == 0, we set the center so that it maps to the ref
        ! element center and set fv_metdet so that it corresponds to the integral
        ! of metdet over the FV subcell. TempestRemap establishes the
@@ -1256,6 +1279,7 @@ end subroutine gfr_f_ref_edges
   subroutine gfr_init_Dmap(elem, gfr)
     use control_mod, only: cubed_sphere_map
     use cube_mod, only: Dmap, ref2sphere
+    use planar_mod, only: plane_Dmap
     use coordinate_systems_mod, only: cartesian3D_t, change_coordinates
 
     type (element_t), intent(in) :: elem(:)
@@ -1272,16 +1296,21 @@ subroutine gfr_init_Dmap(elem, gfr)
     do ie = 1,nelemd
        do j = 1,nf
           do i = 1,nf
-             if (cubed_sphere_map == 2) then
+             if (gfr%is_planar .or. cubed_sphere_map == 0) then
+                call gfr_f_ref_center(nf, i, a)
+                call gfr_f_ref_center(nf, j, b)                
+             else
                 call gfr_f_get_cartesian3d(ie, i, j, sphere)
                 call sphere2ref(elem(ie)%corners3D, sphere, a, b)
-             else
-                call gfr_f_ref_center(nf, i, a)
-                call gfr_f_ref_center(nf, j, b)
              end if
 
-             call Dmap(wrk, a, b, elem(ie)%corners3D, cubed_sphere_map, elem(ie)%cartp, &
-                  elem(ie)%facenum)
+             if (gfr%is_planar) then
+                call plane_Dmap(wrk, a, b, elem(ie)%corners3D, cubed_sphere_map, &
+                     elem(ie)%cartp, elem(ie)%facenum)
+             else
+                call       Dmap(wrk, a, b, elem(ie)%corners3D, cubed_sphere_map, &
+                     elem(ie)%cartp, elem(ie)%facenum)
+             end if
 
              det = wrk(1,1)*wrk(2,2) - wrk(1,2)*wrk(2,1)
 
@@ -2081,10 +2110,14 @@ end subroutine gfr_f_get_corner_latlon
 
   subroutine gfr_f_get_cartesian3d(ie, i, j, p)
     ! Get (x,y,z) of FV point i,j.
+    !   [Projection bug] In the case of sphere goemetry, output is not
+    ! normalized to the sphere.
 
     integer, intent(in) :: ie, i, j
     type (cartesian3D_t), intent(out) :: p
 
+    real(kind=real_kind) :: r
+
     p = gfr%center_f(i,j,ie)
   end subroutine gfr_f_get_cartesian3d
 
@@ -2175,6 +2208,29 @@ subroutine limiter_clip_and_sum(spheremp, qmin, qmax, dp, q)
     call limiter1_clip_and_sum(np*np, spheremp, qmin, qmax, dp, q)
   end subroutine limiter_clip_and_sum
 
+  subroutine ref2plane(corners, a, b, plane)
+    ! For planar geometry.
+
+    use coordinate_systems_mod, only: cartesian3D_t
+
+    type (cartesian3D_t), intent(in) :: corners(4)
+    real(real_kind), intent(in) :: a, b
+    type (cartesian3D_t), intent(out) :: plane
+
+    real(real_kind) :: q(4), c(4,3), p(3)
+    integer :: i
+
+    do i = 1,4
+       c(i,1) = corners(i)%x; c(i,2) = corners(i)%y; c(i,3) = corners(i)%z
+    end do
+    q(1) = (1-a)*(1-b); q(2) = (1+a)*(1-b); q(3) = (1+a)*(1+b); q(4) = (1-a)*(1+b)
+    q = q/four
+    do i = 1,3
+       p(i) = sum(c(:,i)*q)
+    end do
+    plane%x = p(1); plane%y = p(2); plane%z = p(3)
+  end subroutine ref2plane
+
   subroutine ref2spherea_deriv(c, a, b, s_ab, s)
     ! For cubed_sphere_map = 2.
 
@@ -2479,7 +2535,7 @@ subroutine check_areas(par, gfr, elem, nets, nete)
 
     use kinds, only: iulog
     use parallel_mod, only: parallel_t, global_shared_buf, global_shared_sum
-    use physical_constants, only: dd_pi
+    use physical_constants, only: dd_pi, Lx, Ly
     ! Can't use wrap_repro_sum because this routine needs to support
     ! unit testing in an already threaded region.
 #ifdef CAM
@@ -2494,7 +2550,7 @@ subroutine check_areas(par, gfr, elem, nets, nete)
     integer, intent(in) :: nets, nete
 
     integer :: ie, i, j, nf
-    real(kind=real_kind) :: area, sphere_area, re
+    real(kind=real_kind) :: area, area_true, re
 
     nf = gfr%nphys
     do ie = nets,nete
@@ -2509,20 +2565,24 @@ subroutine check_areas(par, gfr, elem, nets, nete)
        global_shared_buf(ie,3) = sum(elem(ie)%spheremp)
     end do
     call repro_sum(global_shared_buf, global_shared_sum, nelemd, nelemd, 3, commid=par%comm)
-    sphere_area = 4*dd_pi
+    if (gfr%is_planar) then
+       area_true = Lx*Ly
+    else
+       area_true = 4*dd_pi
+    end if
     if (par%masterproc) then
        write(iulog,*) 'gfr> area fv raw', global_shared_sum(1), &
-            abs(global_shared_sum(1) - sphere_area)/sphere_area
+            abs(global_shared_sum(1) - area_true)/area_true
        ! fv vs gll
        re = abs(global_shared_sum(2) - global_shared_sum(3))/global_shared_sum(3)
        write(iulog,*) 'gfr> area fv adj', &
-            abs(global_shared_sum(2) - sphere_area)/sphere_area, re
+            abs(global_shared_sum(2) - area_true)/area_true, re
        if (re > 2*eps) then
           write(iulog,*) 'gfr> check_areas ERROR'
           gfr%check_ok = .false.
        end if
        write(iulog,*) 'gfr> area gll   ', &
-            abs(global_shared_sum(3) - sphere_area)/sphere_area
+            abs(global_shared_sum(3) - area_true)/area_true
     end if
     deallocate(gfr%check_areas)
   end subroutine check_areas