From 206a093d1416684686de2cd15af3cead6f0340c1 Mon Sep 17 00:00:00 2001 From: "Andrew M. Bradley" Date: Tue, 31 Oct 2023 18:42:31 -0500 Subject: [PATCH 1/6] Homme: Modify dcmip2016_test1_pg_forcing to handle planar case. --- components/homme/src/share/control_mod.F90 | 4 +- components/homme/src/share/namelist_mod.F90 | 1 + components/homme/src/test_mod.F90 | 13 ++- .../homme/src/test_src/dcmip16_wrapper.F90 | 79 +++++++++++++------ 4 files changed, 66 insertions(+), 31 deletions(-) diff --git a/components/homme/src/share/control_mod.F90 b/components/homme/src/share/control_mod.F90 index 7a2aa3c7221b..0e9494f5a6cd 100644 --- a/components/homme/src/share/control_mod.F90 +++ b/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 diff --git a/components/homme/src/share/namelist_mod.F90 b/components/homme/src/share/namelist_mod.F90 index f0ba74ab2926..1d47090182ba 100644 --- a/components/homme/src/share/namelist_mod.F90 +++ b/components/homme/src/share/namelist_mod.F90 @@ -565,6 +565,7 @@ subroutine readnl(par) test_case == "planar_nonhydro_mtn_wave" .or. & test_case == "planar_schar_mtn_wave" .or. & test_case == "planar_rising_bubble" .or. & + test_case == "planar_rising_bubble_pg2" .or. & test_case == "planar_density_current" .or. & test_case == "planar_baroclinic_instab" .or. & test_case == "planar_moist_rising_bubble" .or. & diff --git a/components/homme/src/test_mod.F90 b/components/homme/src/test_mod.F90 index 5b60e4baf81e..361465218e5d 100644 --- a/components/homme/src/test_mod.F90 +++ b/components/homme/src/test_mod.F90 @@ -24,7 +24,7 @@ module test_mod dcmip2012_test4_init, mtest_init, dcmip2012_test1_1_conv use dcmip16_wrapper, only: dcmip2016_test1, dcmip2016_test2, dcmip2016_test3, & dcmip2016_test1_forcing, dcmip2016_test2_forcing, dcmip2016_test3_forcing, & - dcmip2016_test1_pg, dcmip2016_test1_pg_forcing, dcmip2016_init + dcmip2016_pg_init, dcmip2016_test1_pg, dcmip2016_test1_pg_forcing, dcmip2016_init use held_suarez_mod, only: hs0_init_state use dry_planar_tests, only: planar_hydro_gravity_wave_init, planar_nonhydro_gravity_wave_init @@ -76,8 +76,8 @@ subroutine set_test_initial_conditions(elem, deriv, hybrid, hvcoord, tl, nets, n case('dcmip2012_test2_2'); test_with_forcing = .true. ; case('dcmip2012_test3'); case('dcmip2012_test4'); - case('dcmip2016_test1'); call dcmip2016_init(); test_with_forcing = .true. ; - case('dcmip2016_test1_pg1', 'dcmip2016_test1_pg2', 'dcmip2016_test1_pg3', 'dcmip2016_test1_pg4') + case('dcmip2016_test1', & + 'dcmip2016_test1_pg1', 'dcmip2016_test1_pg2', 'dcmip2016_test1_pg3', 'dcmip2016_test1_pg4') call dcmip2016_init(); test_with_forcing = .true. ; case('dcmip2016_test2'); call dcmip2016_init(); test_with_forcing = .true. ; case('dcmip2016_test3'); call dcmip2016_init(); test_with_forcing = .true. ; @@ -91,7 +91,7 @@ subroutine set_test_initial_conditions(elem, deriv, hybrid, hvcoord, tl, nets, n case('planar_hydro_mtn_wave'); case('planar_nonhydro_mtn_wave'); case('planar_schar_mtn_wave'); - case('planar_rising_bubble'); + case('planar_rising_bubble', 'planar_rising_bubble_pg2') if (bubble_moist) then call dcmip2016_init(); test_with_forcing = .true. ; @@ -146,6 +146,9 @@ subroutine set_test_initial_conditions(elem, deriv, hybrid, hvcoord, tl, nets, n case('planar_nonhydro_mtn_wave'); call planar_nonhydro_mountain_wave_init(elem,hybrid,hvcoord,nets,nete) case('planar_schar_mtn_wave'); call planar_schar_mountain_wave_init(elem,hybrid,hvcoord,nets,nete) case('planar_rising_bubble'); call planar_rising_bubble_init(elem,hybrid,hvcoord,nets,nete) + case('planar_rising_bubble_pg2') + call dcmip2016_pg_init(elem,hybrid,hvcoord,nets,nete,2) + call planar_rising_bubble_init(elem,hybrid,hvcoord,nets,nete) case('planar_density_current'); call planar_density_current_init(elem,hybrid,hvcoord,nets,nete) case('planar_baroclinic_instab'); call planar_baroclinic_instab_init(elem,hybrid,hvcoord,nets,nete) case('planar_moist_rising_bubble'); call planar_moist_rising_bubble_init(elem,hybrid,hvcoord,nets,nete) @@ -251,6 +254,8 @@ subroutine compute_test_forcing(elem,hybrid,hvcoord,nt,ntQ,dt,nets,nete,tl) case('planar_rising_bubble'); if (bubble_moist) call dcmip2016_test1_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl) + case('planar_rising_bubble_pg2'); + if (bubble_moist) call dcmip2016_test1_pg_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl) case('held_suarez0'); do ie=nets,nete diff --git a/components/homme/src/test_src/dcmip16_wrapper.F90 b/components/homme/src/test_src/dcmip16_wrapper.F90 index 719710e0bb27..0806bfe744a5 100644 --- a/components/homme/src/test_src/dcmip16_wrapper.F90 +++ b/components/homme/src/test_src/dcmip16_wrapper.F90 @@ -153,7 +153,7 @@ subroutine dcmip2016_test1(elem,hybrid,hvcoord,nets,nete) end subroutine -subroutine dcmip2016_test1_pg(elem,hybrid,hvcoord,nets,nete,nphys) +subroutine dcmip2016_pg_init(elem,hybrid,hvcoord,nets,nete,nphys) use gllfvremap_mod, only: gfr_init type(element_t), intent(inout), target :: elem(:) ! element array @@ -173,6 +173,18 @@ subroutine dcmip2016_test1_pg(elem,hybrid,hvcoord,nets,nete,nphys) pg_data%q(ncol,nlev,qsize,nelemd)) end if !$omp barrier +end subroutine dcmip2016_pg_init + +subroutine dcmip2016_test1_pg(elem,hybrid,hvcoord,nets,nete,nphys) + use gllfvremap_mod, only: gfr_init + + type(element_t), intent(inout), target :: elem(:) ! element array + type(hybrid_t), intent(in) :: hybrid ! hybrid parallel structure + type(hvcoord_t), intent(inout) :: hvcoord ! hybrid vertical coordinates + integer, intent(in) :: nets,nete ! start, end element index + integer, intent(in) :: nphys ! pgN, N parameter, for physgrid + + call dcmip2016_pg_init(elem,hybrid,hvcoord,nets,nete,nphys) call dcmip2016_test1(elem,hybrid,hvcoord,nets,nete) sample_period = 3600*24 end subroutine dcmip2016_test1_pg @@ -634,15 +646,13 @@ subroutine dcmip2016_test1_pg_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl type(TimeLevel_t), intent(in) :: tl ! time level structure integer, parameter :: iqv = 1 - real(rl), parameter :: one = 1.0_rl + integer, parameter :: test = 1 integer :: i,j,k,ie,qi - real(rl), dimension(np,np,nlev) :: u,v,w,T,p,dp,rho,rho_dry,z,exner_kess,theta_kess - real(rl), dimension(np,np,nlev) :: rho_new,p_pk - real(rl), dimension(nlev) :: u_c,v_c,p_c,qv_c,qc_c,qr_c,rho_c,z_c, th_c + real(rl), dimension(np,np,nlev) :: w,dp + real(rl), dimension(nlev) :: u_c,v_c,p_c,qv_c,qc_c,qr_c,rho_c,z_c,th_c real(rl) :: max_w, max_precl, min_ps - real(rl) :: lat, lon, dz_top(np,np),zi(np,np,nlevp),zi_c(nlevp), ps(np,np), & - wrk(np,np), rd, wrk3(np,np,nlev), wrk4(np,np,nlev,2), wf(np*np,1) + real(rl) :: lat, lon, zi_c(nlevp), wrk3(np,np,nlev), wrk4(np,np,nlev,2), wf(np*np,1) integer :: nf, ncol real(rl), dimension(np,np,nlev) :: dp_fv, p_fv, u_fv, v_fv, T_fv, exner_kess_fv, & @@ -652,14 +662,21 @@ subroutine dcmip2016_test1_pg_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl real(rl), dimension(np,np) :: zs_fv, ps_fv, delta_ps real(rl) :: precl_fv(np,np,1), rcd(6) real(rl), allocatable :: qmin(:,:,:), qmax(:,:,:) - integer :: pbl_type, prec_type - integer, parameter :: test = 1 + logical :: toy_chemistry_on nf = pg_data%nphys ncol = nf*nf - prec_type = dcmip16_prec_type + if (case_planar_bubble) then + toy_chemistry_on = .false. + prec_type = bubble_prec_type + if (qsize .ne. 3) call abortmp('ERROR: moist bubble test requires qsize=3') + else + toy_chemistry_on = .true. + prec_type = dcmip16_prec_type + endif + pbl_type = dcmip16_pbl_type max_w = -huge(rl) @@ -735,9 +752,12 @@ subroutine dcmip2016_test1_pg_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl zi_c = zi_fv(i,j,nlevp:1:-1) th_c = theta_kess_fv(i,j,nlev:1:-1) - call gfr_f_get_latlon(ie, i, j, lat, lon) - ! Get forced versions of u,v,p,qv,qc,qr. rho is constant. + if (case_planar_bubble) then + lat = 0 + else + call gfr_f_get_latlon(ie, i, j, lat, lon) + end if call DCMIP2016_PHYSICS(test, u_c, v_c, p_c, th_c, qv_c, qc_c, qr_c, rho_c, dt, & z_c, zi_c, lat, nlev, precl_fv(i,j,1), pbl_type, prec_type) @@ -748,18 +768,25 @@ subroutine dcmip2016_test1_pg_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl Q_fv(i,j,:,3) = qr_c(nlev:1:-1) theta_kess_fv(i,j,:) = th_c(nlev:1:-1) - do k=1,nlev - call tendency_terminator(lat*rad2dg, lon*rad2dg, Q_fv(i,j,k,4), Q_fv(i,j,k,5), & - dt, ddt_cl(i,j,k), ddt_cl2(i,j,k)) - enddo + if (toy_chemistry_on) then + call gfr_f_get_latlon(ie, i, j, lat, lon) + qi = 4 + do k=1,nlev + call tendency_terminator(lat*rad2dg, lon*rad2dg, Q_fv(i,j,k,qi), Q_fv(i,j,k,qi+1), & + dt, ddt_cl(i,j,k), ddt_cl2(i,j,k)) + end do + end if enddo enddo do i = 1,3 Q_fv(:nf,:nf,:,i) = (rho_dry_fv(:nf,:nf,:)/rho_fv(:nf,:nf,:))*Q_fv(:nf,:nf,:,i) end do - Q_fv(:nf,:nf,:,4) = Q_fv(:nf,:nf,:,4) + dt*ddt_cl(:nf,:nf,:) - Q_fv(:nf,:nf,:,5) = Q_fv(:nf,:nf,:,5) + dt*ddt_cl2(:nf,:nf,:) + if (toy_chemistry_on) then + qi = 4 + Q_fv(:nf,:nf,:,qi ) = Q_fv(:nf,:nf,:,qi ) + dt*ddt_cl (:nf,:nf,:) + Q_fv(:nf,:nf,:,qi+1) = Q_fv(:nf,:nf,:,qi+1) + dt*ddt_cl2(:nf,:nf,:) + end if ! Convert from theta to T w.r.t. new model state. ! Assume hydrostatic pressure pi changed by qv forcing. @@ -803,14 +830,16 @@ subroutine dcmip2016_test1_pg_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl call gfr_f2g_dss(hybrid, elem, nets, nete) call gfr_pg1_reconstruct(hybrid, nt, hvcoord, elem, nets, nete) - call toy_init(rcd) - do ie = nets,nete - do i = 1,2 - wrk4(:,:,:,i) = elem(ie)%state%Q(:,:,:,i+3) + if (toy_chemistry_on) then + call toy_init(rcd) + do ie = nets,nete + do i = 1,2 + wrk4(:,:,:,i) = elem(ie)%state%Q(:,:,:,i+3) + end do + call toy_rcd(wrk4, rcd) end do - call toy_rcd(wrk4, rcd) - end do - call toy_print(hybrid, tl%nstep, rcd) + call toy_print(hybrid, tl%nstep, rcd) + end if if (ftype == 0) then ! Convert FQ from state to Qdp tendency. From 4b54afa19289dff9829aef1aaa64f08ee40f2687 Mon Sep 17 00:00:00 2001 From: "Andrew M. Bradley" Date: Wed, 1 Nov 2023 15:44:19 -0500 Subject: [PATCH 2/6] Homme: Add planar capability to gllfvremap_mod. --- components/homme/src/share/gllfvremap_mod.F90 | 112 ++++++++++++++---- components/homme/src/share/planar_mod.F90 | 8 +- .../homme/src/test_src/dcmip16_wrapper.F90 | 3 +- 3 files changed, 92 insertions(+), 31 deletions(-) diff --git a/components/homme/src/share/gllfvremap_mod.F90 b/components/homme/src/share/gllfvremap_mod.F90 index 1e0ee9b81842..e0e0fa6c4daa 100644 --- a/components/homme/src/share/gllfvremap_mod.F90 +++ b/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 diff --git a/components/homme/src/share/planar_mod.F90 b/components/homme/src/share/planar_mod.F90 index 64bf02a7c0fa..2af18ddcece3 100644 --- a/components/homme/src/share/planar_mod.F90 +++ b/components/homme/src/share/planar_mod.F90 @@ -43,6 +43,7 @@ module planar_mod public :: PlaneEdgeCount public :: PlaneElemCount public :: convert_gbl_index_plane + public :: plane_Dmap ! =============================== ! Private methods @@ -50,7 +51,6 @@ module planar_mod private :: coordinates_atomic private :: metric_atomic private :: coriolis_init_atomic - private :: Dmap contains @@ -310,7 +310,7 @@ subroutine coordinates_atomic(elem,gll_points) end subroutine coordinates_atomic - subroutine Dmap(D, a,b, corners3D, ref_map, cartp, facenum) + subroutine plane_Dmap(D, a,b, corners3D, ref_map, cartp, facenum) real (kind=real_kind), intent(out) :: D(2,2) real (kind=real_kind), intent(in) :: a,b type (cartesian3D_t) :: corners3D(4) !x,y,z coords of element corners @@ -328,7 +328,7 @@ subroutine Dmap(D, a,b, corners3D, ref_map, cartp, facenum) D(2,1) = 0.0D0 D(2,2) = dy/2.0d0 - end subroutine Dmap + end subroutine plane_Dmap ! ========================================= ! metric_atomic: @@ -395,7 +395,7 @@ subroutine metric_atomic(elem,gll_points,alpha) do j=1,np do i=1,np - call Dmap(elem%D(i,j,:,:),1.0D0,1.0D0,elem%corners3D,cubed_sphere_map) + call plane_Dmap(elem%D(i,j,:,:),1.0D0,1.0D0,elem%corners3D,cubed_sphere_map) ! Numerical metric tensor based on analytic D: met = D^T times D ! (D maps between physical plane and reference element) diff --git a/components/homme/src/test_src/dcmip16_wrapper.F90 b/components/homme/src/test_src/dcmip16_wrapper.F90 index 0806bfe744a5..a71a71b2090a 100644 --- a/components/homme/src/test_src/dcmip16_wrapper.F90 +++ b/components/homme/src/test_src/dcmip16_wrapper.F90 @@ -841,7 +841,8 @@ subroutine dcmip2016_test1_pg_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl call toy_print(hybrid, tl%nstep, rcd) end if - if (ftype == 0) then + ! In standalone Homme, for all ftype values, FQ is tendency. + if (.true.) then !(ftype == 0) then ! Convert FQ from state to Qdp tendency. do ie = nets,nete do k = 1,nlev From 804b3b95367a214b6f206f508597149a5a3cf6bd Mon Sep 17 00:00:00 2001 From: "Andrew M. Bradley" Date: Wed, 1 Nov 2023 20:30:07 -0500 Subject: [PATCH 3/6] Homme(xx): Set up C++-vs-F90 standalone planar pg2 test. --- .../namelists/thetanh-moist-bubble-sl-pg2.nl | 72 +++++++++++++++++++ .../namelists/thetanh-moist-bubble-sl.nl | 1 - .../test/reg_test/run_tests/test-list.cmake | 7 +- .../thetanh-moist-bubble-sl-pg2-kokkos.cmake | 6 ++ .../thetanh-moist-bubble-sl-pg2.cmake | 6 ++ 5 files changed, 89 insertions(+), 3 deletions(-) create mode 100644 components/homme/test/reg_test/namelists/thetanh-moist-bubble-sl-pg2.nl create mode 100644 components/homme/test/reg_test/run_tests/thetanh-moist-bubble-sl-pg2-kokkos.cmake create mode 100644 components/homme/test/reg_test/run_tests/thetanh-moist-bubble-sl-pg2.cmake diff --git a/components/homme/test/reg_test/namelists/thetanh-moist-bubble-sl-pg2.nl b/components/homme/test/reg_test/namelists/thetanh-moist-bubble-sl-pg2.nl new file mode 100644 index 000000000000..7cc866df46e9 --- /dev/null +++ b/components/homme/test/reg_test/namelists/thetanh-moist-bubble-sl-pg2.nl @@ -0,0 +1,72 @@ +&ctl_nl + nthreads = 1 + partmethod = 4 + topology = "plane" + geometry = "plane" + test_case = "planar_rising_bubble_pg2" + theta_hydrostatic_mode = .false. + theta_advect_form = 1 + vert_remap_q_alg = 10 + moisture = 'moist' + qsize = 3 + statefreq = 1000 + restartfreq = -1 + runtype = 0 + qsplit = -1 + rsplit = -1 + integration = 'explicit' + tstep_type = 9 + ne_x = 32 + ne_y = 5 + nmax = 300 ! 1800 is good for an actual simulation + tstep = 0.8 + hypervis_order = 2 + hypervis_scaling = 3.2 + nu = 0.01 + nu_top = 0.0 + hypervis_order = 2 ! 2 = hyperviscosity + hypervis_subcycle = 1 ! 1 = no hypervis subcycling + hypervis_subcycle_tom = 1 + limiter_option = 9 + planar_slice=.true. + lx = 20000.0 + ly = 1000.0 + sx = -10000.0 + sy = -500.0 + bubble_zcenter=2000.0 + bubble_ztop=20000.0 + bubble_xyradius=5000.0 + bubble_zradius =1000.0 + bubble_cosine=.true. + bubble_moist=.true. + bubble_T0=300.0 + bubble_dT=4.0 + bubble_moist_drh=1.0 + bubble_prec_type=1 + dt_remap_factor = 5 + transport_alg = 12 + dt_tracer_factor = 10 + hypervis_subcycle_q = 10 + semi_lagrange_hv_q = 1 + semi_lagrange_nearest_point_lev = 256 +/ +&vert_nl + vanalytic = 1 ! set vcoords in initialization routine + vtop = 2.73919e-1 ! vertical coordinate at top of atm (z=10000m) +/ +&analysis_nl + output_dir = "./movies/" ! destination dir for netcdf file + output_timeunits = 0, ! 1=days, 2=hours, 0=timesteps + output_frequency = 300, + output_varnames1 ='T','Th','ps','geo','Q','u','v' !,'Q2','Q3' + interp_type = 0 ! 0=native grid, 1=bilinear + output_type ='netcdf' ! netcdf or pnetcdf + num_io_procs = 1 + interp_nlat = 128 + interp_nlon = 256 + interp_gridtype = 2 ! gauss grid +/ +&prof_inparm + profile_outpe_num = 100 + profile_single_file = .true. +/ diff --git a/components/homme/test/reg_test/namelists/thetanh-moist-bubble-sl.nl b/components/homme/test/reg_test/namelists/thetanh-moist-bubble-sl.nl index ce65723e382a..7dd72926c5f8 100644 --- a/components/homme/test/reg_test/namelists/thetanh-moist-bubble-sl.nl +++ b/components/homme/test/reg_test/namelists/thetanh-moist-bubble-sl.nl @@ -27,7 +27,6 @@ hypervis_order = 2 ! 2 = hyperviscosity hypervis_subcycle = 1 ! 1 = no hypervis subcycling hypervis_subcycle_tom = 1 - se_ftype=2 limiter_option = 9 planar_slice=.true. lx = 20000.0 diff --git a/components/homme/test/reg_test/run_tests/test-list.cmake b/components/homme/test/reg_test/run_tests/test-list.cmake index 175b40678c60..dbb0fc4f98e9 100644 --- a/components/homme/test/reg_test/run_tests/test-list.cmake +++ b/components/homme/test/reg_test/run_tests/test-list.cmake @@ -106,7 +106,9 @@ IF (BUILD_HOMME_THETA_KOKKOS) thetah-nhgw-slice-kokkos.cmake thetanh-nhgw-slice-kokkos.cmake thetanh-moist-bubble-sl.cmake - thetanh-moist-bubble-sl-kokkos.cmake) + thetanh-moist-bubble-sl-kokkos.cmake + thetanh-moist-bubble-sl-pg2.cmake + thetanh-moist-bubble-sl-pg2-kokkos.cmake) IF (HOMMEXX_BFB_TESTING) LIST(APPEND HOMME_ONEOFF_CVF_TESTS thetanh-moist-bubble @@ -115,7 +117,8 @@ IF (BUILD_HOMME_THETA_KOKKOS) thetanh-nhgw thetah-nhgw-slice thetanh-nhgw-slice - thetanh-moist-bubble-sl) + thetanh-moist-bubble-sl + thetanh-moist-bubble-sl-pg2) ENDIF() #cmake/namelist will be built with create-... script diff --git a/components/homme/test/reg_test/run_tests/thetanh-moist-bubble-sl-pg2-kokkos.cmake b/components/homme/test/reg_test/run_tests/thetanh-moist-bubble-sl-pg2-kokkos.cmake new file mode 100644 index 000000000000..7df1c1bfe912 --- /dev/null +++ b/components/homme/test/reg_test/run_tests/thetanh-moist-bubble-sl-pg2-kokkos.cmake @@ -0,0 +1,6 @@ +# Kokkos version of pre-existing F90 test. +SET(TEST_NAME thetanh-moist-bubble-sl-pg2-kokkos) +SET(EXEC_NAME theta-l-nlev20-native-kokkos) +SET(NUM_CPUS 16) +SET(NAMELIST_FILES ${HOMME_ROOT}/test/reg_test/namelists/thetanh-moist-bubble-sl-pg2.nl) +SET(NC_OUTPUT_FILES planar_rising_bubble_pg21.nc) diff --git a/components/homme/test/reg_test/run_tests/thetanh-moist-bubble-sl-pg2.cmake b/components/homme/test/reg_test/run_tests/thetanh-moist-bubble-sl-pg2.cmake new file mode 100644 index 000000000000..e8d8463f97b2 --- /dev/null +++ b/components/homme/test/reg_test/run_tests/thetanh-moist-bubble-sl-pg2.cmake @@ -0,0 +1,6 @@ +# theta-l NH moist bubble test using SL transport. +SET(TEST_NAME thetanh-moist-bubble-sl-pg2) +SET(EXEC_NAME theta-l-nlev20-native) +SET(NUM_CPUS 16) +SET(NAMELIST_FILES ${HOMME_ROOT}/test/reg_test/namelists/thetanh-moist-bubble-sl-pg2.nl) +SET(NC_OUTPUT_FILES planar_rising_bubble_pg21.nc) From e1a1356544581bc38410be468c876a939ea907bd Mon Sep 17 00:00:00 2001 From: "Andrew M. Bradley" Date: Wed, 1 Nov 2023 21:55:00 -0500 Subject: [PATCH 4/6] Homme: Add gllfvremap_planar_ut unit test. --- .../homme/src/share/cxx/GllFvRemapImpl.cpp | 2 + .../homme/src/share/cxx/mpi/Connectivity.cpp | 2 +- .../homme/src/share/cxx/mpi/Connectivity.hpp | 2 +- .../homme/src/share/gllfvremap_util_mod.F90 | 117 +++++++++++++----- .../homme/test/unit_tests/CMakeLists.txt | 26 ++-- .../test_execs/share_kokkos_ut/CMakeLists.txt | 1 + .../share_kokkos_ut/geometry_interface.F90 | 22 +++- .../thetal_kokkos_ut/CMakeLists.txt | 1 + .../thetal_kokkos_ut/gllfvremap_interface.F90 | 10 +- .../thetal_kokkos_ut/gllfvremap_ut.cpp | 7 +- .../thetal_test_interface.F90 | 91 ++++++++++++-- 11 files changed, 219 insertions(+), 62 deletions(-) diff --git a/components/homme/src/share/cxx/GllFvRemapImpl.cpp b/components/homme/src/share/cxx/GllFvRemapImpl.cpp index 0736bea763ae..6148f69cfa9c 100644 --- a/components/homme/src/share/cxx/GllFvRemapImpl.cpp +++ b/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 diff --git a/components/homme/src/share/cxx/mpi/Connectivity.cpp b/components/homme/src/share/cxx/mpi/Connectivity.cpp index e11b2f65819e..2d0841c942e2 100644 --- a/components/homme/src/share/cxx/mpi/Connectivity.cpp +++ b/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, diff --git a/components/homme/src/share/cxx/mpi/Connectivity.hpp b/components/homme/src/share/cxx/mpi/Connectivity.hpp index 3cd350f40ac4..898729f86beb 100644 --- a/components/homme/src/share/cxx/mpi/Connectivity.hpp +++ b/components/homme/src/share/cxx/mpi/Connectivity.hpp @@ -156,7 +156,7 @@ class Connectivity ExecViewManaged::HostMirror h_ucon_ptr; ExecViewManaged d_ucon_dir_ptr; ExecViewManaged::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; diff --git a/components/homme/src/share/gllfvremap_util_mod.F90 b/components/homme/src/share/gllfvremap_util_mod.F90 index 4e56f1a514a4..e3f80ff95598 100644 --- a/components/homme/src/share/gllfvremap_util_mod.F90 +++ b/components/homme/src/share/gllfvremap_util_mod.F90 @@ -47,6 +47,8 @@ module gllfvremap_util_mod type (PhysgridData_t), private :: pg_data + logical :: is_sphere + public :: & ! Test gllfvremap's main API. gfr_check_api, & @@ -59,6 +61,22 @@ module gllfvremap_util_mod contains + function sphere2cart(sphere) result(cart) + use coordinate_systems_mod, only: spherical_polar_t, cartesian3D_t, change_coordinates + + type (spherical_polar_t), intent(in) :: sphere + type (cartesian3D_t) :: cart + + if (is_sphere) then + cart = change_coordinates(sphere) + else + ! See conventions established in planar_mod::coordinates_atomic. + cart%x = sphere%lon + cart%y = sphere%lat + cart%z = 0 + end if + end function sphere2cart + subroutine init(nphys) ! Init pg_data. @@ -102,8 +120,8 @@ subroutine set_gll_state(hvcoord, elem, nt1, nt2) ! function on the sphere. use dimensions_mod, only: nlev, qsize - use physical_constants, only: g, dd_pi - use coordinate_systems_mod, only: cartesian3D_t, change_coordinates + use physical_constants, only: g, dd_pi, Lx, Ly + use coordinate_systems_mod, only: cartesian3D_t use hybvcoord_mod, only: hvcoord_t use element_ops, only: get_field @@ -113,35 +131,57 @@ subroutine set_gll_state(hvcoord, elem, nt1, nt2) type (State_t) :: s1 type (cartesian3D_t) :: p - real(kind=real_kind) :: wr(np,np,nlev,2) + real(kind=real_kind) :: wr(np,np,nlev,2), fx, fy integer :: i, j, k, q, tl + if (.not. is_sphere) then + fx = 2*dd_pi/Lx + fy = 2*dd_pi/Ly + end if + elem%state%Q(:,:,:,1) = zero ! moisture tracer is 0 do j = 1,np do i = 1,np - p = change_coordinates(elem%spherep(i,j)) - do k = 1,nlev - do q = 2,qsize - elem%state%Q(i,j,k,q) = one + & - half*sin((half + modulo(q,2))*p%x)* & - sin((half + modulo(q,3))*1.5d0*p%y)* & - sin((-2.3d0 + modulo(q,5))*p%z) + p = sphere2cart(elem%spherep(i,j)) + if (is_sphere) then + do k = 1,nlev + do q = 2,qsize + elem%state%Q(i,j,k,q) = one + & + half*sin((half + modulo(q,2))*p%x)* & + sin((half + modulo(q,3))*1.5d0*p%y)* & + sin((-2.3d0 + modulo(q,5))*(0.7d0 + p%z)) + end do end do - end do - s1%ps(i,j) = 1.0d3*(one + 0.05d0*sin(two*p%x+half)*sin(p%y+1.5d0)*sin(3*p%z+2.5d0)) - s1%phis(i,j) = one + half*sin(p%x-half)*sin(half*p%y+2.5d0)*sin(2*p%z-2.5d0) - do k = 1,nlev - ! u, v have to be set carefully because they are - ! converted to contravariant velocity. Thus, at the - ! poles, we need u, v to make sense to measure OOA - ! correctly. - wr(i,j,k,1) = sin(p%x)*sin(1.5*p%y)*cos(1.7*p%z) - wr(i,j,k,2) = sin(half*p%x)*sin(1.5*p%y)*cos(half*dd_pi*p%z) - elem%derived%omega_p(i,j,k) = wr(i,j,k,1) - end do - do k = 1,nlev - s1%T(i,j,k) = one + half*sin(p%x+1.5d0)*sin(1.5d0*p%y+half)*sin(two*p%z-half) - end do + s1%ps(i,j) = 1.0d3*(one + 0.05d0*sin(two*p%x+half)*sin(p%y+1.5d0)*sin(3*p%z+2.5d0)) + s1%phis(i,j) = one + half*sin(p%x-half)*sin(half*p%y+2.5d0)*sin(2*p%z-2.5d0) + do k = 1,nlev + ! u, v have to be set carefully because they are converted to + ! contravariant velocity. Thus, at the poles, we need u, v to make + ! sense to measure OOA correctly. + wr(i,j,k,1) = sin(p%x)*sin(1.5*p%y)*cos(1.7*p%z) + wr(i,j,k,2) = sin(half*p%x)*sin(1.5*p%y)*cos(half*dd_pi*p%z) + elem%derived%omega_p(i,j,k) = wr(i,j,k,1) + end do + do k = 1,nlev + s1%T(i,j,k) = one + half*sin(p%x+1.5d0)*sin(1.5d0*p%y+half)*sin(two*p%z-half) + end do + else + do k = 1,nlev + do q = 2,qsize + elem%state%Q(i,j,k,q) = one + half*sin(fx*p%x)*sin(fy*p%y) + end do + end do + s1%ps(i,j) = 1.0d3*(one + 0.05d0*sin(two*fx*p%x+half)*sin(fy*p%y+1.5d0)) + s1%phis(i,j) = one + half*sin(fx*p%x-half)*sin(fy*p%y+2.5d0) + do k = 1,nlev + wr(i,j,k,1) = sin(fx*p%x)*sin(2*fy*fy*p%y) + wr(i,j,k,2) = sin(fx*p%x)*sin(2*fy*p%y) + elem%derived%omega_p(i,j,k) = wr(i,j,k,1) + end do + do k = 1,nlev + s1%T(i,j,k) = one + half*sin(fx*p%x+1.5d0)*sin(fy*p%y+half) + end do + end if end do end do s1%u = wr(:,:,:,1) @@ -168,6 +208,21 @@ subroutine set_gll_state(hvcoord, elem, nt1, nt2) end do end subroutine set_gll_state + function make_tendency(p) result(t) + use coordinate_systems_mod, only: cartesian3D_t + use physical_constants, only: dd_pi, Lx, Ly + + type (cartesian3D_t), intent(in) :: p + real(kind=real_kind) :: t, r + + if (is_sphere) then + r = sqrt(p%x*p%x + p%y*p%y + p%z*p%z) + t = 0.25_real_kind*sin(3.2*p%x)*sin(4.2*p%y)*sin(0.7 + 2.3*p%z) + else + t = 0.25_real_kind*sin(2*dd_pi*p%x/Lx)*sin(2*dd_pi*p%y/Ly) + end if + end function make_tendency + function run(hybrid, hvcoord, elem, nets, nete, nphys, tendency) result(nerr) ! Run 3 convergence and property-preservation whole-mesh ! tests. See below for descriptions of the three tests. @@ -175,7 +230,7 @@ function run(hybrid, hvcoord, elem, nets, nete, nphys, tendency) result(nerr) use kinds, only: iulog use hybvcoord_mod, only: hvcoord_t use dimensions_mod, only: nlev, qsize - use coordinate_systems_mod, only: cartesian3D_t, change_coordinates + use coordinate_systems_mod, only: cartesian3D_t use element_ops, only: get_temperature, get_field use prim_driver_base, only: applyCAMforcing_tracers use prim_advance_mod, only: applyCAMforcing_dynamics @@ -238,7 +293,7 @@ function run(hybrid, hvcoord, elem, nets, nete, nphys, tendency) result(nerr) do i = 1,nf col = nf*(j-1) + i call gfr_f_get_cartesian3d(ie, i, j, p) - f = 0.25_real_kind*sin(3.2*p%x)*sin(4.2*p%y)*sin(2.3*p%z) + f = make_tendency(p) pg_data%uv(col,:,:,ie) = f/dt pg_data%T(col,:,ie) = f/dt ! no moisture adjustment => no dp3d adjustment @@ -285,8 +340,8 @@ function run(hybrid, hvcoord, elem, nets, nete, nphys, tendency) result(nerr) if (tendency .and. q > 1) then do j = 1,np do i = 1,np - p = change_coordinates(elem(ie)%spherep(i,j)) - tend(i,j,:) = 0.25_real_kind*sin(3.2*p%x)*sin(4.2*p%y)*sin(2.3*p%z) + p = sphere2cart(elem(ie)%spherep(i,j)) + tend(i,j,:) = make_tendency(p) end do end do end if @@ -525,7 +580,7 @@ function gfr_check_api(hybrid, nets, nete, hvcoord, elem) result(nerr) use kinds, only: iulog use hybvcoord_mod, only: hvcoord_t - use control_mod, only: ftype + use control_mod, only: ftype, geometry use gllfvremap_mod type (hybrid_t), intent(in) :: hybrid @@ -536,6 +591,8 @@ function gfr_check_api(hybrid, nets, nete, hvcoord, elem) result(nerr) integer :: nphys, ftype_in, ftype_idx, boost_idx, nerr logical :: boost_pg1 + is_sphere = trim(geometry) /= 'plane' + ftype_in = ftype nerr = 0 diff --git a/components/homme/test/unit_tests/CMakeLists.txt b/components/homme/test/unit_tests/CMakeLists.txt index 4dd08a288cec..23c6b1a3169d 100644 --- a/components/homme/test/unit_tests/CMakeLists.txt +++ b/components/homme/test/unit_tests/CMakeLists.txt @@ -1,23 +1,25 @@ SET(UNITTESTER_DIR ${CMAKE_CURRENT_SOURCE_DIR} PARENT_SCOPE) +function(cxx_unit_test_add_test test_name target_name NUM_CPUS) + SET(EXTRA_ARGS ${ARGN}) + set(TMP ${USE_MPI_OPTIONS}) + separate_arguments(TMP) + if (USE_MPI_RUN_SCRIPT) + ADD_TEST(${test_name} ${USE_MPI_RUN_SCRIPT} ${NUM_CPUS} ./${target_name} ${EXTRA_ARGS}) + else() + ADD_TEST(${test_name} ${USE_MPIEXEC} -n ${NUM_CPUS} ${TMP} ./${target_name} ${EXTRA_ARGS}) + endif() + SET_TESTS_PROPERTIES(${test_name} PROPERTIES LABELS "unit") +endfunction() + macro(cxx_unit_test target_name target_f90_srcs target_cxx_srcs include_dirs config_defines NUM_CPUS) ADD_EXECUTABLE(${target_name} ${UNITTESTER_DIR}/tester.cpp ${target_f90_srcs} ${target_cxx_srcs}) #add exec to test_execs ,baseline, and check targets in makefile ADD_DEPENDENCIES(test-execs ${target_name}) ADD_DEPENDENCIES(baseline ${target_name}) ADD_DEPENDENCIES(check ${target_name}) - #IF (${NUM_CPUS} EQUAL 1) - # ADD_TEST(${target_name}_test ${target_name}) - #ELSE() - set(TMP ${USE_MPI_OPTIONS}) - separate_arguments(TMP) - if (USE_MPI_RUN_SCRIPT) - ADD_TEST(${target_name}_test ${USE_MPI_RUN_SCRIPT} ${NUM_CPUS} ./${target_name}) - else() - ADD_TEST(${target_name}_test ${USE_MPIEXEC} -n ${NUM_CPUS} ${TMP} ./${target_name}) - endif() - #ENDIF() - SET_TESTS_PROPERTIES(${target_name}_test PROPERTIES LABELS "unit") + + cxx_unit_test_add_test(${target_name}_test ${target_name} ${NUM_CPUS}) TARGET_LINK_LIBRARIES(${target_name} timing ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES}) target_link_libraries(${target_name} kokkos) diff --git a/components/homme/test_execs/share_kokkos_ut/CMakeLists.txt b/components/homme/test_execs/share_kokkos_ut/CMakeLists.txt index 2b84af38e68d..c0a3bb3c9d68 100644 --- a/components/homme/test_execs/share_kokkos_ut/CMakeLists.txt +++ b/components/homme/test_execs/share_kokkos_ut/CMakeLists.txt @@ -54,6 +54,7 @@ SET (BOUNDARY_EXCHANGE_UT_F90_SRCS ${SRC_SHARE_DIR}/control_mod.F90 ${SRC_SHARE_DIR}/coordinate_systems_mod.F90 ${SRC_SHARE_DIR}/cube_mod.F90 + ${SRC_SHARE_DIR}/planar_mod.F90 ${SRC_SHARE_DIR}/derivative_mod_base.F90 ${SRC_SHARE_DIR}/dimensions_mod.F90 ${SRC_SHARE_DIR}/edge_mod_base.F90 diff --git a/components/homme/test_execs/share_kokkos_ut/geometry_interface.F90 b/components/homme/test_execs/share_kokkos_ut/geometry_interface.F90 index 78acd80d7bbf..cb693e285649 100644 --- a/components/homme/test_execs/share_kokkos_ut/geometry_interface.F90 +++ b/components/homme/test_execs/share_kokkos_ut/geometry_interface.F90 @@ -42,8 +42,10 @@ end subroutine initmp_f90 subroutine init_cube_geometry_f90 (ne_in) bind(c) - use iso_c_binding, only : c_int + use iso_c_binding, only: c_int + use control_mod, only: topology use cube_mod, only: CubeTopology, CubeElemCount, CubeEdgeCount + use planar_mod, only: PlaneTopology, PlaneElemCount, PlaneEdgeCount use dimensions_mod, only: nelem, nelemd, npart, np use dimensions_mod, only: ne use gridgraph_mod, only: allocate_gridvertex_nbrs @@ -56,10 +58,18 @@ subroutine init_cube_geometry_f90 (ne_in) bind(c) ! Locals ! integer :: ie, num_edges + logical :: is_cube + + is_cube = trim(topology) == 'cube' ne = ne_in - nelem = CubeElemCount() - num_edges = CubeEdgeCount() + if (is_cube) then + nelem = CubeElemCount() + num_edges = CubeEdgeCount() + else + nelem = PlaneElemCount() + num_edges = PlaneEdgeCount() + end if if (nelem < par%nprocs) then call abortmp('Error: too many MPI tasks. Run the test with less MPI ranks, or increase ne.') @@ -74,7 +84,11 @@ subroutine init_cube_geometry_f90 (ne_in) bind(c) enddo ! Generate mesh connectivity - call CubeTopology(GridEdge, GridVertex) + if (is_cube) then + call CubeTopology(GridEdge, GridVertex) + else + call PlaneTopology(GridEdge,GridVertex) + end if ! Set number of partitions before partitioning mesh npart = par%nprocs diff --git a/components/homme/test_execs/thetal_kokkos_ut/CMakeLists.txt b/components/homme/test_execs/thetal_kokkos_ut/CMakeLists.txt index d798c13960c0..6e3b5080aecd 100644 --- a/components/homme/test_execs/thetal_kokkos_ut/CMakeLists.txt +++ b/components/homme/test_execs/thetal_kokkos_ut/CMakeLists.txt @@ -347,3 +347,4 @@ ELSE() ENDIF() cxx_unit_test (gllfvremap_ut "${GLLFVREMAP_UT_F90_SRCS}" "${GLLFVREMAP_UT_CXX_SRCS}" "${GLLFVREMAP_UT_INCLUDE_DIRS}" "${CONFIG_DEFINES}" ${NUM_CPUS}) TARGET_LINK_LIBRARIES(gllfvremap_ut thetal_kokkos_ut_lib) +cxx_unit_test_add_test(gllfvremap_planar_ut gllfvremap_ut ${NUM_CPUS} "hommexx -planar") diff --git a/components/homme/test_execs/thetal_kokkos_ut/gllfvremap_interface.F90 b/components/homme/test_execs/thetal_kokkos_ut/gllfvremap_interface.F90 index 8f54867cbe88..940b72fb7f4c 100644 --- a/components/homme/test_execs/thetal_kokkos_ut/gllfvremap_interface.F90 +++ b/components/homme/test_execs/thetal_kokkos_ut/gllfvremap_interface.F90 @@ -11,7 +11,7 @@ subroutine init_gllfvremap_f90(ne, hyai, hybi, hyam, hybm, ps0, dvv, mp, qsize_i is_sphere) bind(c) use iso_c_binding, only: c_double use hybvcoord_mod, only: set_layer_locations - use thetal_test_interface, only: init_f90 + use thetal_test_interface, only: init_f90, init_planar_f90 use theta_f2c_mod, only: init_elements_c use edge_mod_base, only: initEdgeBuffer, edge_g, initEdgeSBuffer use prim_advection_base, only: edgeAdvQminmax @@ -28,13 +28,15 @@ subroutine init_gllfvremap_f90(ne, hyai, hybi, hyam, hybm, ps0, dvv, mp, qsize_i integer :: edgesz - if (.not. is_sphere) print *, "NOT IMPL'ED YET" - qsize = qsize_in use_moisture = .true. theta_hydrostatic_mode = .false. - call init_f90(ne, hyai, hybi, hyam, hybm, dvv, mp, ps0) + if (is_sphere) then + call init_f90(ne, hyai, hybi, hyam, hybm, dvv, mp, ps0) + else + call init_planar_f90(ne+1, ne, hyai, hybi, hyam, hybm, dvv, mp, ps0) + end if call init_elements_c(nelemd) edgesz = max((qsize+3)*nlev+2,6*nlev+1) diff --git a/components/homme/test_execs/thetal_kokkos_ut/gllfvremap_ut.cpp b/components/homme/test_execs/thetal_kokkos_ut/gllfvremap_ut.cpp index faf1b57eca4c..0f14b0c3e55a 100644 --- a/components/homme/test_execs/thetal_kokkos_ut/gllfvremap_ut.cpp +++ b/components/homme/test_execs/thetal_kokkos_ut/gllfvremap_ut.cpp @@ -171,7 +171,6 @@ struct Session { printf("seed %u\n", seed); parse_command_line(); - assert(is_sphere); // planar isn't available in Hxx yet auto& c = Context::singleton(); @@ -250,7 +249,7 @@ struct Session { private: static std::shared_ptr s_session; - // gllfvremap_ut hommexx -ne NE -qsize QSIZE + // gllfvremap_ut hommexx -ne NE -qsize QSIZE [-planar] void parse_command_line () { const bool am_root = get_comm().root(); ne = 2; @@ -271,6 +270,7 @@ struct Session { } } ne = std::max(2, std::min(128, ne)); + if ( ! is_sphere) ne = std::max(4, ne); qsize = std::max(1, std::min(QSIZE_D, qsize)); if ( ! ok && am_root) printf("gllfvremap_ut> Failed to parse command line, starting with: %s\n", @@ -282,7 +282,8 @@ struct Session { #else 0; #endif - printf("gllfvremap_ut> bfb %d ne %d qsize %d\n", bfb, ne, qsize); + printf("gllfvremap_ut> bfb %d ne %d qsize %d sphere %d\n", + bfb, ne, qsize, int(is_sphere)); } } }; diff --git a/components/homme/test_execs/thetal_kokkos_ut/thetal_test_interface.F90 b/components/homme/test_execs/thetal_kokkos_ut/thetal_test_interface.F90 index 1065580235f9..458a75628df6 100644 --- a/components/homme/test_execs/thetal_kokkos_ut/thetal_test_interface.F90 +++ b/components/homme/test_execs/thetal_kokkos_ut/thetal_test_interface.F90 @@ -27,10 +27,6 @@ subroutine init_f90 (ne, hyai, hybi, hyam, hybm, dvv, mp, ps0) bind(c) use dimensions_mod, only: nelemd, nlev, nlevp, np use geometry_interface_mod, only: initmp_f90, init_cube_geometry_f90, init_connectivity_f90 use geometry_interface_mod, only: par, elem - use hybvcoord_mod, only: set_layer_locations - use element_state, only: allocate_element_arrays, setup_element_pointers_ie, & - nlev_tom, nu_scale_top - use mass_matrix_mod, only: mass_matrix use quadrature_mod, only: gausslobatto, quadrature_t use physical_constants, only: scale_factor, scale_factor_inv, laplacian_rigid_factor, rearth, rrearth ! @@ -43,7 +39,7 @@ subroutine init_f90 (ne, hyai, hybi, hyam, hybm, dvv, mp, ps0) bind(c) ! ! Locals ! - integer :: ie, k + integer :: ie type (quadrature_t) :: gp scale_factor = rearth @@ -64,6 +60,88 @@ subroutine init_f90 (ne, hyai, hybi, hyam, hybm, dvv, mp, ps0) bind(c) do ie=1,nelemd call cube_init_atomic(elem(ie),gp%points) enddo + + call init_common(hyai, hybi, hyam, hybm, dvv, mp, ps0) + end subroutine init_f90 + + subroutine init_planar_f90 (ne_x_in, ne_y_in, hyai, hybi, hyam, hybm, dvv, mp, ps0) bind(c) + ! Alternative to init_f90 to create a planar model. + use control_mod, only: cubed_sphere_map, geometry, topology + use planar_mod, only: plane_init_atomic, plane_set_corner_coordinates + use derivative_mod, only: derivinit + use dimensions_mod, only: nelemd, nlev, nlevp, np, ne_x, ne_y + use geometry_interface_mod, only: initmp_f90, init_cube_geometry_f90, init_connectivity_f90, & + par, elem + use quadrature_mod, only: gausslobatto, quadrature_t + use physical_constants, only: scale_factor, scale_factor_inv, laplacian_rigid_factor, & + Sx, Sy, Lx, Ly, dx, dy, dx_ref, dy_ref, domain_size + ! + ! Inputs + ! + integer (kind=c_int), intent(in) :: ne_x_in, ne_y_in + real (kind=real_kind), intent(in) :: hyai(nlevp), hybi(nlevp), hyam(nlev), hybm(nlev) + real (kind=real_kind), intent(in) :: ps0 + real (kind=real_kind), intent(out) :: dvv(np,np), mp(np,np) + ! + ! Locals + ! + integer :: ie + type (quadrature_t) :: gp + + ne_x = ne_x_in + ne_y = ne_y_in + + scale_factor = 1 + scale_factor_inv = 1 + laplacian_rigid_factor = 0 + topology = 'plane' + geometry = 'plane' + Lx = 10000 + Ly = 5000 + Sx = -5000 + Sy = 2500 + + domain_size = Lx * Ly + dx = Lx/ne_x + dy = Ly/ne_y + dx_ref = 1.0D0/ne_x + dy_ref = 1.0D0/ne_y + + call derivinit(deriv) + + call initmp_f90() + call init_cube_geometry_f90(ne_x) ! ne_x is unused + call init_connectivity_f90() + + cubed_sphere_map = 2 + gp=gausslobatto(np) ! GLL points + do ie=1,nelemd + call plane_set_corner_coordinates(elem(ie)) + end do + do ie=1,nelemd + call plane_init_atomic(elem(ie),gp%points) + enddo + + call init_common(hyai, hybi, hyam, hybm, dvv, mp, ps0) + end subroutine init_planar_f90 + + subroutine init_common(hyai, hybi, hyam, hybm, dvv, mp, ps0) + use element_state, only: allocate_element_arrays, setup_element_pointers_ie + use mass_matrix_mod, only: mass_matrix + use hybvcoord_mod, only: set_layer_locations + use dimensions_mod, only: nelemd, nlev, nlevp, np + use geometry_interface_mod, only: par, elem + ! + ! Inputs + ! + real (kind=real_kind), intent(in) :: hyai(nlevp), hybi(nlevp), hyam(nlev), hybm(nlev) + real (kind=real_kind), intent(in) :: ps0 + real (kind=real_kind), intent(out) :: dvv(np,np), mp(np,np) + ! + ! Locals + ! + integer :: ie, k + call allocate_element_arrays(nelemd) call mass_matrix(par,elem) @@ -89,8 +167,7 @@ subroutine init_f90 (ne, hyai, hybi, hyam, hybm, dvv, mp, ps0) bind(c) call set_layer_locations (hvcoord,.false.,.false.) deriv%dvv = dvv - - end subroutine init_f90 + end subroutine init_common subroutine init_geo_views_f90 (d_ptr, dinv_ptr, & phis_ptr, gradphis_ptr, fcor_ptr, & From 50d641ddd145c1188d857114a22f3a6bd5c74e2a Mon Sep 17 00:00:00 2001 From: "Andrew M. Bradley" Date: Mon, 6 Nov 2023 12:49:15 -0600 Subject: [PATCH 5/6] Homme: Fix latitude in dcmip2016_test1_forcing. This is an old bug in the original forcing function. The physgrid variant, dcmip2016_test1_forcing, is fine. --- components/homme/src/test_src/dcmip16_wrapper.F90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/components/homme/src/test_src/dcmip16_wrapper.F90 b/components/homme/src/test_src/dcmip16_wrapper.F90 index a71a71b2090a..de26960d5b69 100644 --- a/components/homme/src/test_src/dcmip16_wrapper.F90 +++ b/components/homme/src/test_src/dcmip16_wrapper.F90 @@ -528,7 +528,11 @@ subroutine dcmip2016_test1_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl) th_c = theta_kess(i,j,nlev:1:-1) ! get forced versions of u,v,p,qv,qc,qr. rho is constant - lat=0.0 + if (case_planar_bubble) then + lat = 0 + else + lat = elem(ie)%spherep(i,j)%lat + end if call DCMIP2016_PHYSICS(test, u_c, v_c, p_c, th_c, qv_c, qc_c, qr_c, rho_c, dt, z_c, zi_c, lat, nlev, & precl(i,j,ie), pbl_type, prec_type) @@ -951,7 +955,7 @@ subroutine dcmip2016_test2_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl, t zi_c = zi (i,j,nlevp:1:-1) th_c = theta_kess(i,j,nlev:1:-1) - lat=0.0 + lat=0.0 ! unused in test 2 ! get forced versions of u,v,p,qv,qc,qr. rho is constant call DCMIP2016_PHYSICS(test, u_c, v_c, p_c, th_c, qv_c, qc_c, qr_c, rho_c, dt, z_c, zi_c, lat, nlev, & precl(i,j,ie), pbl_type, prec_type) From 8324c66b8ab506cd9e537309e476518381661ab1 Mon Sep 17 00:00:00 2001 From: "Andrew M. Bradley" Date: Mon, 13 Nov 2023 16:09:45 -0600 Subject: [PATCH 6/6] Homme/SL: Address a warning. --- components/homme/src/share/compose/compose_cedr.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/homme/src/share/compose/compose_cedr.cpp b/components/homme/src/share/compose/compose_cedr.cpp index 2b497db6647b..0aeaebf9487c 100644 --- a/components/homme/src/share/compose/compose_cedr.cpp +++ b/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) {