Several fixes for RK3, including tiling when using AGRIF
Context
- Branches impacted: main
- Reference configuration/test case: AGRIF_DEMO, VORTEX
Analysis
AGRIF issues
AGRIF_DEMO and VORTEX are not reproducible when the tiling is used, as noted in !387 (merged).
There are three issues:
-
There is overlap in the calculation of
zFu
/zFv
between tiles#151 (closed) provides a general description of this fundamental problem with the current tiling implementation.
In this specific case, the overlap occurs as follows:
-
zFu
/zFv
are defined in stprk3_stg.F90 within the tiling loop:IF( ln_tile ) CALL dom_tile_start ! [tiling] DYN tiling loop (2) DO jtile = 1, nijtile IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) ... DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) ! advective transport zFu(ji,jj,jk) = e2u(ji,jj)*e3u(ji,jj,jk,Kmm) * ( uu(ji,jj,jk,Kmm) + zub(ji,jj)*umask(ji,jj,jk) ) zFv(ji,jj,jk) = e1v(ji,jj)*e3v(ji,jj,jk,Kmm) * ( vv(ji,jj,jk,Kmm) + zvb(ji,jj)*vmask(ji,jj,jk) ) END_3D
-
In the AGRIF_DEMO case where
ln_dynadv_vec = .true.
andln_shuman = .false.
,zFu
/zFv
are next accessed by tra_adv_trp where they are updated, within the same tiling loop instp_rk3_stg
:... IF( .NOT.ln_shuman ) THEN ! No Shuman averaging- tra_adv_trp can be in tiling loop ! !* Update or Compute (VIF) advective transport CALL tra_adv_trp( kstp, kstg, nit000, Kbb, Kmm, Kaa, Krhs, zFu, zFv, zFw ) ! ! BBL coefficients required for both passive- and active-tracer transport within ! ! the BBL (stage 3 only, requires uu, vv, gdept at Kmm) IF( kstg == 3 .AND. ln_trabbl ) CALL bbl( kstp, nit000, Kbb, Kmm ) ! ENDIF ! END DO ! IF( ln_tile ) CALL dom_tile_stop
-
These operations on
zFu
/zFv
can be expressed as:IF( ln_tile ) CALL dom_tile_start DO jtile = 1, nijtile IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) ! stp_rk3_stg DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) zFu(ji,jj,jk) = ... zFv(ji,jj,jk) = ... END_3D ! tra_adv_trp DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) zFu(ji,jj,jk) = zFu(ji,jj,jk) + ... zFv(ji,jj,jk) = zFv(ji,jj,jk) + ... END_3D END DO IF( ln_tile ) CALL dom_tile_stop
The DO loops include halo points, meaning that the calculations will overlap between tiles. Because
zFu
/zFv
are updated in-place, this would normally mean that some points are updated more than once per timestep, giving the incorrect answer (see the first slide in #151 (closed)). The results do not change however, becausezFu
/zFv
are reset at the start of every tiling loop iteration by the calculation of their initial value. -
When
key_agrif
is activated, this tiling loop must be split around the call toAgrif_Sponge_dyn
since this subroutine will not work with the tiling:# if defined key_agrif END DO IF( ln_tile ) CALL dom_tile_stop IF(.NOT. Agrif_Root() ) THEN ! AGRIF: sponge ==> momentum and tracer RHS IF( kstg == 3 ) CALL Agrif_Sponge_dyn ENDIF IF( ln_tile ) CALL dom_tile_start ! [tiling] DYN tiling loop (2, continued) DO jtile = 1, nijtile IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) # endif
-
The two operations on
zFu
/zFv
are now in separate tiling loops and can be expressed as:IF( ln_tile ) CALL dom_tile_start DO jtile = 1, nijtile IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) ! stp_rk3_stg DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) zFu(ji,jj,jk) = ... zFv(ji,jj,jk) = ... END_3D END DO IF( ln_tile ) CALL dom_tile_stop IF( ln_tile ) CALL dom_tile_start DO jtile = 1, nijtile IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) ! tra_adv_trp DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) zFu(ji,jj,jk) = zFu(ji,jj,jk) + ... zFv(ji,jj,jk) = zFv(ji,jj,jk) + ... END_3D END DO IF( ln_tile ) CALL dom_tile_stop
zFu
/zFv
are now updated in-place within the 2nd tiling loop, while the calculation of their initial value is performed in the 1st tiling loop. Tiling now changes the result, becausezFu
/zFv
are no longer reset to their initial values at the start of each tiling loop iteration. This particular issue is demonstrated by the first slide in #151 (closed).
-
-
An incorrect workaround for tiling was used for the call to
Agrif_avm
inzdf_phy
In zdf_phy we have:
#if defined key_agrif ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2) IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile IF( l_zdfsh2 ) CALL Agrif_avm ENDIF #endif ! ! !* start from turbulent closure values DO_3D( 0, 0, 0, 0, 2, jpkm1 ) avt(ji,jj,jk) = avt_k(ji,jj,jk) avm(ji,jj,jk) = avm_k(ji,jj,jk) END_3D
avm_k
is not updated on the child zoom until the last tile, butavm
is updated fromavm_k
for every tile. This means that values ofavm
on the child zoom are mostly those from the previous timestep. -
An incorrect workaround for tiling was used for the zeroing of
ww
on AGRIF child zooms inwzv
In wzv (both RK3 and MLF versions) we have:
#if defined key_agrif IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile IF( .NOT. AGRIF_Root() ) THEN ! ! Mask vertical velocity at first/last columns/row ! inside computational domain (cosmetic) DO jk = 1, jpkm1 IF( lk_west ) THEN ! --- West --- ! DO ji = mi0(2+nn_hls,nn_hls), mi1(2+nn_hls,nn_hls) DO jj = 1, jpj pww(ji,jj,jk) = 0._wp END DO END DO ENDIF ...
While described as "cosmetic", this zeroing does in fact affect the results. The present solution for the tiling therefore doesn't work- the zeroing must be done for each tile.
Other issues
A number of other bugs were identified during the investigation of this issue:
-
The
lbc_lnk
associated withll_zAimp2
in tra_adv_fct may result in an MPI error due to not being called on all MPI processes.IF( ln_zad_Aimp ) THEN IF( MAXVAL( ABS( wi(T2D(1),:) ) ) > 0._wp ) THEN IF( kn_fct_imp == 1 ) ll_zAimp1 = .TRUE. IF( kn_fct_imp == 2 ) ll_zAimp2 = .TRUE. ENDIF END IF
The condition
MAXVAL( ABS( wi(T2D(1),:) ) ) > 0._wp
may not be met on all processes, in which casell_zAimp2
will not be true on all MPI processes and they will not all call the associatedlbc_lnk
. In this case, the processes that do not call it will move on to the nextlbc_lnk
and data will be exchanged from completely different arrays. As a result, an MPI error will be raised due to inconsistent message sizes. -
A similar MPI communications issue can occur with the tiling when certain tile sizes are used.
The tile size is limited by the size of the MPI domain, such that the former must be smaller than the latter when accounting for haloes. The maximum tile size is thus:
imax = Ni_0 - (2 * nn_hls) - 1 jmax = Nj_0 - (2 * nn_hls) - 1
In the event that
Ni_0
andNj_0
are too small to permit tiling (imax < 1
/jmax < 1
), then one tile is used (nn_ltile_i = Ni_0
/nn_ltile_j = Nj_0
). Otherwise, the tile size will be limited to the maximum tile size.Ni_0
andNj_0
will not be the same for all MPI processes if the global domain does not divide equally between them. This can result in situations where some MPI processes have domains that are too small for tiling, and others do not. In this case, some processes will only have one tile (nijtile = 1
) and others will have several (nijtile > 1
).There are several parts of the code that run only for the final tile, e.g. from
zdf_phy
:IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN CALL lbc_lnk( 'zdfphy', avm, 'W', 1.0_wp ) ENDIF
If some MPI processes have
nijtile = 1
and othersnijtile > 1
, then on the first tile (ntile = 1
) the former processes will call thislbc_lnk
but the latter will move onto the nextlbc_lnk
call. This means data from two differentlbc_lnk
calls will be exchanged, resulting in the same MPI error as the previous issue.This issue was occuring for the first child zoom of AGRIF_DEMO when run with an 8 x 4 processor decomposition- the southern processes had
jpj = 10
(jmax = 1
) and the northern processes hadjpj = 9
(jmax = 0
). -
There are a couple of uses of
ln_ldfeiv_dia
in ldf_eiv_trp_RK3 that should bel_ldfeiv_dia
.While this did not affect the diagnostic results, it meant some unnecessary allocation and computation was being performed on RK3 steps 1 and 2.
Fix
AGRIF issues
-
dom_tile_copy
calls have been added forzFu
/zFv
andwi
The
dom_tile_copy
functionality was introduced by #151 (closed) to address tile overlap issues without requiring extensive refactoring. It is used to address the issue reported here (for arrayszFu
/zFv
) and a similar issue in the calculation ofwi
when using the adaptive-implicit vertical advection (ln_zad_Aimp = .true.
).Note that
zFu
/zFv
have been redeclared asALLOCATABLE
, as this attribute is required for arrays passed in todom_tile_copy
. -
The tiling loop has been split around the call to
Agrif_avm
, as for other calls to AGRIF routines:#if defined key_agrif END DO IF( ln_tile ) CALL dom_tile_stop ! ! !== ocean Kz ==! (avt, avs, avm) ! ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2) IF( l_zdfsh2 .AND. .NOT. Agrif_Root() ) CALL Agrif_avm IF( ln_tile ) CALL dom_tile_start ! [tiling] ZDF tiling loop (continued) DO jtile = 1, nijtile IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) #else ! !== ocean Kz ==! (avt, avs, avm) #endif
This means that the tiling loop around the call to
zdf_phy
in the main timestepping subroutines has had to be moved intozdf_phy
itself. -
Calls to
in_hdom
are used to ensure that zeroing is only performed for points within the current tile- IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile IF( .NOT. AGRIF_Root() ) THEN ! ! Mask vertical velocity at first/last columns/row ! inside computational domain (cosmetic) DO jk = 1, jpkm1 IF( lk_west ) THEN ! --- West --- ! - DO ji = mi0(2+nn_hls,nn_hls), mi1(2+nn_hls,nn_hls) + DO ji = mi0(1+nn_hls,nn_hls), mi1(1+nn_hls,nn_hls) - DO jj = 1, jpj + DO jj = 2, jpj-1 - pww(ji,jj,jk) = 0._wp + IF( in_hdom(ji, jj, khls=1) ) pww(ji,jj,jk) = 0._wp END DO END DO ENDIF ...
The DO loop bounds have also been modified to reflect that
ww
is now only updated on the first halo points.The
in_hdom
function was introduced by #367 (closed) and is used for a similar purpose in two BDY subroutines.
Other issues
-
The conditional statement for the
lbc_lnk
intra_adv_fct
has been modified to be more consistent across MPI processes:- IF( ll_zAimp2 ) THEN + IF( ln_zad_Aimp .AND. kn_fct_imp == 2 ) THEN CALL lbc_lnk( 'traadv_fct', ztFu, 'U', -1.0_wp, ztFv, 'V', -1.0_wp, ztFw, 'T', 1.0_wp, zta_up1, 'T', 1.0_wp ) ELSE CALL lbc_lnk( 'traadv_fct', ztFu, 'U', -1.0_wp, ztFv, 'V', -1.0_wp, zta_up1, 'T', 1.0_wp ) ENDIF
-
Tiling is now disabled if any MPI process has only one tile (
nijtile = 1
)This has been implemented in
dom_tile_init
indomtile.F90
. -
ln_ldfeiv_dia
has been replaced byl_ldfeiv_dia
inldf_eiv_trp_RK3
A note on the AGRIF issues and tile overlap
It is important to note that diagnosing problems with the tiling can be very time consuming. Tile overlap is by far the most significant issue with the current tiling implementation. There is no easy way of identifying which arrays are affected by tile overlap or which code changes may have caused it, so the developer must have a decent understanding of how the tiling may be changing the path through the code.
This is not something we should require of developers- the goal has always been for the tiling to be as transparent and non-intrusive as possible. As things stand, this is not the case. While dom_tile_copy
provides a reasonably user-friendly solution to tile overlap, diagnosing it remains very difficult. This raises the question of whether the tiling is fit for purpose, and whether it should be turned on in SETTE if it still has the potential to cause developer headaches.