Opened 6 months ago

Closed 5 months ago

#2696 closed Defect (fixed)

Investigate halo1 - halo2 differences with ICB

Reported by: acc Owned by: acc
Priority: low Milestone:
Component: ICB Version: trunk
Severity: minor Keywords:
Cc: Branch review:
MP ready?: Task progress:

Description (last modified by nemo)

Context

SETTE tests with ORCA2_ICE_PISCES give different results with nn_hls=1 and nn_hls=2 if icebergs are active. This ticket is to record investigations and potential solutions.

Analysis

This is a tricky one to pin down because differences in run.stat outputs don't show for many timesteps (84 on my test system), the initial differences are very small and much of the calculation in icb routines is done through function calls that are difficult to instrument. First observations are:

  • Running with icebergs on but ln_passive_mode =.true. will give identical results. This isn't surprising since the icebergs do not affect the ocean in this mode. It does, however, provide an option for running with icebergs on for the purposes of assessing computational performance with nn_hls=2 vs nn_hls=1.
  • With ddefault settings, the earliest detectable diffferences can be found after 6 or 7 timesteps in the values returned from icb_accel (icbdyn.F90). These are caused by 14th or 15th decimal place differences in the interpolation weights returned by icb_utl_pos which is called from icb_utl_interp which is called by icb_accel. These differences occur despite the iceberg positions appearing to be 'identical' in the two cases. iceberg positions are given as a continuous local gridcell value so they are offset by one between the two cases but appear otherwise identical to 15 decimal places. Nevertheless the interpolation weights in the two cases have small differences which eventually lead to significant differences. The onset of differences can be delayed by ensuring the calculation of weights is identical in both cases by removing the additional offset in the nn_hls=2 case. For example:
  • icbutl.F90

    svn diff icbutl.F90
     
    279365         ! note that here there is no +0.5 added 
    280366         ! since we're looking for four T points containing quadrant we're in of 
    281367         ! current T cell 
    282          kii = MAX(0, INT( pi        )) 
    283          kij = MAX(0, INT( pj        ))    ! T-point 
    284          zwi = pi - REAL(kii,wp) 
    285          zwj = pj - REAL(kij,wp) 
     368         kii = MAX(0, INT( pi-(nn_hls-1)        )) 
     369         kij = MAX(0, INT( pj-(nn_hls-1)        ))    ! T-point 
     370         zwi = pi-(nn_hls-1) - REAL(kii,wp) 
     371         zwj = pj-(nn_hls-1) - REAL(kij,wp) 
    286372      CASE ( 'U' ) 
    287          kii = MAX(0, INT( pi-0.5_wp )) 
    288          kij = MAX(0, INT( pj        ))    ! U-point 
    289          zwi = pi - 0.5_wp - REAL(kii,wp) 
    290          zwj = pj - REAL(kij,wp) 
     373         kii = MAX(0, INT( pi-(nn_hls-1)-0.5_wp )) 
     374         kij = MAX(0, INT( pj-(nn_hls-1)        ))    ! U-point 
     375         zwi = pi-(nn_hls-1) - 0.5_wp - REAL(kii,wp) 
     376         zwj = pj-(nn_hls-1) - REAL(kij,wp) 
    291377      CASE ( 'V' ) 
    292          kii = MAX(0, INT( pi        )) 
    293          kij = MAX(0, INT( pj-0.5_wp ))    ! V-point 
    294          zwi = pi - REAL(kii,wp) 
    295          zwj = pj - 0.5_wp - REAL(kij,wp) 
     378         kii = MAX(0, INT( pi-(nn_hls-1)        )) 
     379         kij = MAX(0, INT( pj-(nn_hls-1)-0.5_wp ))    ! V-point 
     380         zwi = pi-(nn_hls-1) - REAL(kii,wp) 
     381         zwj = pj-(nn_hls-1) - 0.5_wp - REAL(kij,wp) 
    296382      CASE ( 'F' ) 
    297          kii = MAX(0, INT( pi-0.5_wp )) 
    298          kij = MAX(0, INT( pj-0.5_wp ))    ! F-point 
    299          zwi = pi - 0.5_wp - REAL(kii,wp) 
    300          zwj = pj - 0.5_wp - REAL(kij,wp) 
     383         kii = MAX(0, INT( pi-(nn_hls-1)-0.5_wp )) 
     384         kij = MAX(0, INT( pj-(nn_hls-1)-0.5_wp ))    ! F-point 
     385         zwi = pi-(nn_hls-1) - 0.5_wp - REAL(kii,wp) 
     386         zwj = pj-(nn_hls-1) - 0.5_wp - REAL(kij,wp) 
    301387      END SELECT 
     388      kii = kii + (nn_hls-1) 
     389      kij = kij + (nn_hls-1) 
    302390      ! 
    303391      ! compute weight 
    304392      pw(1) = (1._wp-zwi) * (1._wp-zwj) 

There is further complication in that icb_utl_interp also estimates the local ssh gradient using:

         CALL icb_utl_pos( pi+0.1_wp, pj, 'T', iiTp, ijTp, zwTp, zmskTp )
         CALL icb_utl_pos( pi-0.1_wp, pj, 'T', iiTm, ijTm, zwTm, zmskTm )
         !
         IF ( .NOT. PRESENT(pe1) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )
         pssh_i = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) -   &
            &       icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. )  ) / ( 0.2_wp * pe1 )

To make these calculations identical, the (nn_hls-1) factor has to be subtracted before the offset is applied. An alternative version of icb_utl_pos can be created for this purpose.

These changes further delay the onset of differences but do not eliminate them.

  • A test has been carried out where the least signifcant decimal places of the iceberg positions and speeds are set to zero. I.e.:
  • icbdyn.F90

    svn diff icbdyn.F90
     
    9696         ! position using di/dt & djdt   !   V2  in m/s 
    9797         zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zuvel1 + zdt_2 * zax1 
    9898         zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1 
     99         zuvel2 = int(zuvel2 * 1.e10_wp,KIND=8)*1.e-10_wp 
     100         zvvel2 = int(zvvel2 * 1.e10_wp,KIND=8)*1.e-10_wp 
     101         zxi2 = int(zxi2 * 1.e10_wp,KIND=8)*1.e-10_wp 
     102         zyj2 = int(zyj2 * 1.e10_wp,KIND=8)*1.e-10_wp 
    99103         ! 
    100104         CALL icb_ground( berg, zxi2, zxi1, zu1,   & 
    101105            &                   zyj2, zyj1, zv1, ll_bounced ) 
     
    112116         !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3) 
    113117         zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2 
    114118         zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2 
     119         zuvel3 = int(zuvel3 * 1.e10_wp,KIND=8)*1.e-10_wp 
     120         zvvel3 = int(zvvel3 * 1.e10_wp,KIND=8)*1.e-10_wp 
     121         zxi3 = int(zxi3 * 1.e10_wp,KIND=8)*1.e-10_wp 
     122         zyj3 = int(zyj3 * 1.e10_wp,KIND=8)*1.e-10_wp 
    115123         ! 
    116124         CALL icb_ground( berg, zxi3, zxi1, zu2,   & 
    117125            &                   zyj3, zyj1, zv2, ll_bounced ) 
     
    128136         !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3 
    129137         zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3 
    130138         zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3 
     139         zuvel4 = int(zuvel4 * 1.e10_wp,KIND=8)*1.e-10_wp 
     140         zvvel4 = int(zvvel4 * 1.e10_wp,KIND=8)*1.e-10_wp 
     141         zxi4 = int(zxi4 * 1.e10_wp,KIND=8)*1.e-10_wp 
     142         zyj4 = int(zyj4 * 1.e10_wp,KIND=8)*1.e-10_wp 
    131143 
    132144         CALL icb_ground( berg, zxi4, zxi1, zu3,   & 
    133145            &                   zyj4, zyj1, zv3, ll_bounced ) 
     
    151163         CALL icb_ground( berg, zxi_n, zxi1, zuvel_n,   & 
    152164            &                   zyj_n, zyj1, zvvel_n, ll_bounced ) 
    153165 
    154          pt%uvel = zuvel_n                        !** save in berg structure 
    155          pt%vvel = zvvel_n 
    156          pt%xi   = zxi_n 
    157          pt%yj   = zyj_n 
     166         !pt%uvel = zuvel_n                        !** save in berg structure 
     167         !pt%vvel = zvvel_n 
     168         pt%uvel = int(zuvel_n * 1.e10_wp,KIND=8)*1.e-10_wp                        !** save in berg structure 
     169         pt%vvel = int(zvvel_n * 1.e10_wp,KIND=8)*1.e-10_wp 
     170         pt%xi   = int( zxi_n * 1.e10_wp,KIND=8)*1.e-10_wp 
     171         pt%yj   = int( zyj_n * 1.e10_wp,KIND=8)*1.e-10_wp 

With these changes differences in run.stat values do not occur until timestep 255.

Recommendation

It seems the differences are due to numerical differences arising from the offset in local gridcell number. A solution is to remove the factor of (nn_hls-1) from the initial iceberg positions and to add this factor back on whenever continuous positions are converted back to gridcell indices. Doing this ensures that all positional calculations are done with exactly the same numbers and results in identical solutions with either nn_hls=1 or nn_hls=2. Changes have been made and tested on a ticket branch pending review and independent testing.

Commit History (3)

ChangesetAuthorTimeChangeLog
15088acc2021-07-06T15:03:34+02:00

Check in ICB changes from 2021/ticket2696_icb_halo1_halo2_compatibility (#2696) to restore halo1 - halo2 equivalence with ORCA2_ICE_PISCES and icebergs. The benefit won't be apparent until the reproducibility issue with icedyn_rhg_evp.F90 is sorted. This fix closes #2696

15043acc2021-06-22T13:26:16+02:00

#2696 .A possible fix for obtaining identical results with icebergs active and with either nn_hls=1 or nn_hls=2. This solution works by ensuring the real-valued iceberg position calculations are done with identical numbers in both cases. Iceberg positions are stored as a continuous local grid-cell value and a factor of (nn_hls-1) is subtracted from initial positions to ensure equivalence in both cases. This factor needs to be added back whenever iceberg positions are converted to local indices. The changes in this commit are enough to ensure a full ORCA2_ICE_PISCES SETTE test is passed in both case (992 timesteps) and provides identical results. Suitability of this solution for inclusion in the trunk needs to be discussed and verified on other platforms and compilers.

15040acc2021-06-22T12:32:58+02:00

Development branch to test fixes for halo1 - halo2 compatibility with active icebergs. #2696

Change History (6)

comment:1 Changed 6 months ago by acc

In 15040:

Development branch to test fixes for halo1 - halo2 compatibility with active icebergs. #2696

comment:2 Changed 6 months ago by acc

In 15043:

#2696 .A possible fix for obtaining identical results with icebergs active and with either
nn_hls=1 or nn_hls=2. This solution works by ensuring the real-valued iceberg position
calculations are done with identical numbers in both cases. Iceberg positions are stored
as a continuous local grid-cell value and a factor of (nn_hls-1) is subtracted from
initial positions to ensure equivalence in both cases. This factor needs to be added
back whenever iceberg positions are converted to local indices. The changes in this
commit are enough to ensure a full ORCA2_ICE_PISCES SETTE test is passed in both case
(992 timesteps) and provides identical results. Suitability of this solution for inclusion
in the trunk needs to be discussed and verified on other platforms and compilers.

comment:3 Changed 6 months ago by acc

  • Description modified (diff)

comment:4 Changed 6 months ago by nemo

  • Description modified (diff)

comment:5 Changed 6 months ago by nemo

  • Description modified (diff)

comment:6 Changed 5 months ago by acc

  • Owner set to acc
  • Resolution set to fixed
  • Status changed from new to closed

In 15088:

Check in ICB changes from 2021/ticket2696_icb_halo1_halo2_compatibility (#2696) to restore halo1 - halo2 equivalence with ORCA2_ICE_PISCES and icebergs. The benefit won't be apparent until the reproducibility issue with icedyn_rhg_evp.F90 is sorted. This fix closes #2696

Note: See TracTickets for help on using tickets.