Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 6 commits into from
Nov 14, 2023
Merged

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

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
206a093
Homme: Modify dcmip2016_test1_pg_forcing to handle planar case.
ambrad Oct 31, 2023
4b54afa
Homme: Add planar capability to gllfvremap_mod.
ambrad Nov 1, 2023
804b3b9
Homme(xx): Set up C++-vs-F90 standalone planar pg2 test.
ambrad Nov 2, 2023
e1a1356
Homme: Add gllfvremap_planar_ut unit test.
ambrad Nov 2, 2023
50d641d
Homme: Fix latitude in dcmip2016_test1_forcing.
ambrad Nov 6, 2023
8324c66
Homme/SL: Address a warning.
ambrad Nov 13, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions components/homme/src/share/control_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down
2 changes: 2 additions & 0 deletions components/homme/src/share/cxx/GllFvRemapImpl.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand why we need this barrier (and the one below)

Copy link
Member Author

@ambrad ambrad Nov 6, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Neither is clearly needed. However, because there is still a potential source of very rare nondeterminism lurking in our code, I've been inspecting old code as I work on new related code. All the data I've collected about this rare nondeterminism can't rule out its occurring in physics-dynamics grid remap, so I read through this file again while working on the planar physgrid case. I decided to put these two barriers in based on this fresh read-through.

The one below, in my analysis, is needed if uv_dim != 3, but I decided to put it in w/o a check on that condition for the reason I'll state in a minute. (In practice, uv_dim == 3.)

I put this one in the for the same reason, although here again the loop patterns in the two adjacent loops arguably don't require it, though in a much subtler way than the usual case of identical loop patterns.

The principle I decided to follow here was strict adherence to protecting work arrays, disregarding implicit protection from loop patterns, just to make reasoning about determinism here simpler. I believe with these two additional team_barriers, that principle is then followed throughout the whole file.

}

// (u,v)
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion components/homme/src/share/cxx/mpi/Connectivity.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
2 changes: 1 addition & 1 deletion components/homme/src/share/cxx/mpi/Connectivity.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
112 changes: 86 additions & 26 deletions components/homme/src/share/gllfvremap_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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(:)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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(:)
Expand All @@ -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)

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i assume this is needed for SL?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is for the metric terms in dynamics-physics grid remap. This PR is to support pg2 in planar mode.

! 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.

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
Loading