Skip to content

Homme: Fix planar-mode pgN get_latlon routines. #2982

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

Merged
merged 2 commits into from
Sep 9, 2024
Merged

Homme: Fix planar-mode pgN get_latlon routines. #2982

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ac1f0ed
Homme: Fix planar-mode pgN get_latlon routines.
ambrad Aug 30, 2024
cbf149b
Hommexx: Use portable PRIx64 instead of lx in printf.
ambrad Sep 3, 2024
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
6 changes: 4 additions & 2 deletions components/homme/src/share/cxx/utilities/InternalDiagnostics.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@
#include "TimeLevel.hpp"
#include "mpi/Comm.hpp"

#include <cinttypes>

namespace Homme {

namespace {
Expand Down Expand Up @@ -50,10 +52,10 @@ void print_global_state_hash (const std::string& label) {
all_reduce_HashType(comm.mpi_comm(), accum, gaccum, 5);
if (comm.root()) {
for (int i = 0; i < NUM_TIME_LEVELS; ++i)
fprintf(stderr, "hxxhash> %14d %1d %16lx (E %s)\n",
fprintf(stderr, "hxxhash> %14d %1d %16" PRIx64 " (E %s)\n",
tl.nstep, i, gaccum[i], label.c_str());
for (int i = 0; i < Q_NUM_TIME_LEVELS; ++i)
fprintf(stderr, "hxxhash> %14d %1d %16lx (T %s)\n",
fprintf(stderr, "hxxhash> %14d %1d %16" PRIx64 " (T %s)\n",
tl.nstep, i, gaccum[NUM_TIME_LEVELS+i], label.c_str());
}
}
Expand Down
97 changes: 75 additions & 22 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 @@ -2088,9 +2088,14 @@ subroutine gfr_f_get_latlon(ie, i, j, lat, lon)

type (spherical_polar_t) :: p

p = change_coordinates(gfr%center_f(i,j,ie))
lat = p%lat
lon = p%lon
if (gfr%is_planar) then
lon = gfr%center_f(i,j,ie)%x
lat = gfr%center_f(i,j,ie)%y
else
p = change_coordinates(gfr%center_f(i,j,ie))
lat = p%lat
lon = p%lon
end if
end subroutine gfr_f_get_latlon

subroutine gfr_f_get_corner_latlon(ie, i, j, c, lat, lon)
Expand All @@ -2103,9 +2108,14 @@ subroutine gfr_f_get_corner_latlon(ie, i, j, c, lat, lon)

type (spherical_polar_t) :: p

p = change_coordinates(gfr%corners_f(c,i,j,ie))
lat = p%lat
lon = p%lon
if (gfr%is_planar) then
lon = gfr%corners_f(c,i,j,ie)%x
lat = gfr%corners_f(c,i,j,ie)%y
else
p = change_coordinates(gfr%corners_f(c,i,j,ie))
lat = p%lat
lon = p%lon
end if
end subroutine gfr_f_get_corner_latlon

subroutine gfr_f_get_cartesian3d(ie, i, j, p)
Expand Down Expand Up @@ -2318,6 +2328,7 @@ subroutine set_ps_Q(elem, nets, nete, timeidx, qidx, nlev)
! Make up a test function for use in unit tests.

use coordinate_systems_mod, only: cartesian3D_t, change_coordinates
use physical_constants, only: dd_pi, Lx, Ly

type (element_t), intent(inout) :: elem(:)
integer, intent(in) :: nets, nete, timeidx, qidx, nlev
Expand All @@ -2329,7 +2340,13 @@ subroutine set_ps_Q(elem, nets, nete, timeidx, qidx, nlev)
do ie = nets, nete
do j = 1,np
do i = 1,np
p = change_coordinates(elem(ie)%spherep(i,j))
if (gfr%is_planar) then
p%x = elem(ie)%spherep(i,j)%lon*2*dd_pi/Lx
p%y = elem(ie)%spherep(i,j)%lat*2*dd_pi/Ly
p%z = 0
else
p = change_coordinates(elem(ie)%spherep(i,j))
end if
elem(ie)%state%ps_v(i,j,timeidx) = &
1.0d3*(1 + 0.05*sin(2*p%x+0.5)*sin(p%y+1.5)*sin(3*p%z+2.5))
q = 0.5*(1 + sin(3*p%x)*sin(3*p%y)*sin(4*p%z))
Expand Down Expand Up @@ -2609,7 +2626,8 @@ function check(par, dom_mt, gfr, elem, verbose) result(nerr)
logical, intent(in) :: verbose

real(kind=real_kind) :: a, b, rd, x, y, f0(np*np), f1(np*np), g(np,np), &
wf(np*np), wg(np,np), qmin, qmax, qmin1, qmax1
wf(np*np), wg(np,np), qmin, qmax, qmin1, qmax1, lat, lon, latext(2), &
lonext(2), tol
integer :: nf, nf2, ie, i, j, k, iremap, info, ilimit, it
real(kind=real_kind), allocatable :: Qdp_fv(:,:), ps_v_fv(:,:), &
qmins(:,:,:), qmaxs(:,:,:)
Expand Down Expand Up @@ -2669,6 +2687,31 @@ function check(par, dom_mt, gfr, elem, verbose) result(nerr)
nerr = nerr + 1
write(iulog,*) 'gfr> Dinv', ie, rd
end if
if (gfr%is_planar) then
lonext(1) = minval(elem(ie)%spherep(:,:)%lon)
lonext(2) = maxval(elem(ie)%spherep(:,:)%lon)
latext(1) = minval(elem(ie)%spherep(:,:)%lat)
latext(2) = maxval(elem(ie)%spherep(:,:)%lat)
tol = max(lonext(2) - lonext(1), latext(2) - latext(1))*1.e4*eps
do i = 1, nf
do j = 1, nf
call gfr_f_get_latlon(ie, i, j, lat, lon)
if ( lon < lonext(1) .or. lon > lonext(2) .or. &
lat < latext(1) .or. lat > latext(2)) then
write(iulog,*) 'gfr> planar latlon ctr', ie, latext, lonext, lat, lon
nerr = nerr + 1
end if
do k = 1, 4
call gfr_f_get_corner_latlon(ie, i, j, k, lat, lon)
if ( lon < lonext(1) - tol .or. lon > lonext(2) + tol .or. &
lat < latext(1) - tol .or. lat > latext(2) + tol) then
write(iulog,*) 'gfr> planar latlon crn', ie, latext, lonext, lat, lon
nerr = nerr + 1
end if
end do
end do
end do
end if

! Check that FV -> GLL -> FV recovers the original FV values exactly
! (with no DSS and no limiter).
Expand Down Expand Up @@ -2839,8 +2882,8 @@ function test_sphere2ref() result(nerr)
use kinds, only: iulog

type (cartesian3D_t) :: corners(4), sphere
real (real_kind) :: refin(2), refout(2), err
integer :: i, j, n, nerr
real (real_kind) :: refin(2), refout(2), err, tol
integer :: i, j, n, trial, nerr

nerr = 0

Expand All @@ -2854,18 +2897,28 @@ function test_sphere2ref() result(nerr)
do i = 1,n
refin(1) = -1 + (1.0_real_kind/(n-1))*(i-1)
do j = 1,n
refin(2) = -1 + (1.0_real_kind/(n-1))*(j-1)
call ref2sphere(corners, refin(1), refin(2), sphere)
call sphere2ref(corners, sphere, refout(1), refout(2))
err = abs(refin(1) - refout(1)) + abs(refin(2) - refout(2))
if (err > 15*eps .or. &
maxval(abs(refout)) > 1 + 5*eps .or. &
any(refout /= refout)) then
write(iulog,*) refin(1), refin(2)
write(iulog,*) refout(1), refout(2)
write(iulog,*) err
nerr = nerr + 1
end if
do trial = 1, 2
refin(2) = -1 + (1.0_real_kind/(n-1))*(j-1)
call ref2sphere(corners, refin(1), refin(2), sphere)
if (trial == 2) then
sphere%x = 0.9d0*sphere%x
sphere%y = 0.9d0*sphere%y
sphere%z = 0.9d0*sphere%z
end if
call sphere2ref(corners, sphere, refout(1), refout(2))
err = abs(refin(1) - refout(1)) + abs(refin(2) - refout(2))
tol = 15*eps
if (trial == 2) tol = 1e-6
if (err > tol .or. &
maxval(abs(refout)) > 1 + 5*eps .or. &
any(refout /= refout)) then
write(iulog,*) 'gfr> test_sphere2ref trial', trial
write(iulog,*) refin(1), refin(2)
write(iulog,*) refout(1), refout(2)
write(iulog,*) err
nerr = nerr + 1
end if
end do
end do
end do
if (nerr /= 0) write(iulog,*) 'test_sphere2ref FAILED'
Expand Down