Precision truncation to facilitate solution reproducibility with implicit time stepping on unstructured mesh

model/inp/ww3_grid.inp

@@ -358,6 +358,8 @@ $                           JGS_LIMITER  : TRUE: Use total (quasi-steady: limits
 $                                          FALSE: default 
 $                           JGS_LIMITER_FUNC  : 1 - old limiter (default)
 $                                               2 - alternatnive limiter
+$                           JGS_TRUNK_DIGITS  : [Only with switch TRNK] Number of digits of precision to truncate solution to at the end of iterative solver. 
+$                                               Can be used to force bit-for-bit reproducibility.
 $                           SETUP_APPLY_WLV  : Compute wave setup (TRUE/FALSE, default TRUE) 
 $                           SOLVERTHR_SETUP  : Solver threshold for setup computations (default 1E-6)
 $                           CRIT_DEP_SETUP  : Critical depth for setup computations (default 0.1) 

model/src/w3gdatmd.F90

@@ -1055,6 +1055,7 @@ MODULE W3GDATMD
     LOGICAL :: B_JGS_LIMITER
     LOGICAL :: B_JGS_USE_JACOBI
     LOGICAL :: B_JGS_BLOCK_GAUSS_SEIDEL
+    INTEGER :: B_JGS_TRUNK_DIGITS
     INTEGER :: B_JGS_MAXITER
     INTEGER :: B_JGS_LIMITER_FUNC
     REAL*8  :: B_JGS_PMIN
@@ -1418,6 +1419,7 @@ MODULE W3GDATMD
   LOGICAL, POINTER :: B_JGS_BLOCK_GAUSS_SEIDEL
   INTEGER, POINTER :: B_JGS_MAXITER
   INTEGER, POINTER :: B_JGS_LIMITER_FUNC
+  INTEGER, POINTER :: B_JGS_TRUNK_DIGITS
   REAL(8), POINTER :: B_JGS_PMIN
   REAL(8), POINTER :: B_JGS_DIFF_THR
   REAL(8), POINTER :: B_JGS_NORM_THR
@@ -2864,6 +2866,7 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST )
     B_JGS_NORM_THR => MPARS(IMOD)%SCHMS%B_JGS_NORM_THR
     B_JGS_NLEVEL => MPARS(IMOD)%SCHMS%B_JGS_NLEVEL
     B_JGS_SOURCE_NONLINEAR => MPARS(IMOD)%SCHMS%B_JGS_SOURCE_NONLINEAR
+    B_JGS_TRUNK_DIGITS => MPARS(IMOD)%SCHMS%B_JGS_TRUNK_DIGITS
     RETURN
     !
     ! Formats

model/src/w3gridmd.F90

@@ -942,6 +942,7 @@ MODULE W3GRIDMD
   REAL*8  :: JGS_PMIN
   REAL*8  :: JGS_DIFF_THR
   REAL*8  :: JGS_NORM_THR
+  INTEGER :: JGS_TRUNK_DIGITS
   REAL*8  :: SOLVERTHR_SETUP
   REAL*8  :: CRIT_DEP_SETUP
   !
@@ -1109,6 +1110,7 @@ MODULE W3GRIDMD
        JGS_NORM_THR,                              &
        JGS_NLEVEL,                                &
        JGS_SOURCE_NONLINEAR,                      &
+       JGS_TRUNK_DIGITS,                          &
        SETUP_APPLY_WLV, SOLVERTHR_SETUP,          &
        CRIT_DEP_SETUP
   NAMELIST /MISC/ CICE0, CICEN, LICE, XSEED, FLAGTR, XP, XR, &
@@ -2480,6 +2482,7 @@ SUBROUTINE W3GRID()
     JGS_NORM_THR = 1.E-20
     JGS_NLEVEL = 0
     JGS_SOURCE_NONLINEAR = .FALSE.
+    JGS_TRUNK_DIGITS = 5
     ! read data from the unstructured devoted namelist
     CALL READNL ( NDSS, 'UNST', STATUS )
 
@@ -2496,6 +2499,7 @@ SUBROUTINE W3GRID()
     B_JGS_NORM_THR = JGS_NORM_THR
     B_JGS_NLEVEL = JGS_NLEVEL
     B_JGS_SOURCE_NONLINEAR = JGS_SOURCE_NONLINEAR
+    B_JGS_TRUNK_DIGITS = JGS_TRUNK_DIGITS
 
     nbSel=0
 
@@ -3363,6 +3367,9 @@ SUBROUTINE W3GRID()
            JGS_DIFF_THR,                               &
            JGS_NORM_THR,                               &
            JGS_NLEVEL,                                 &
+#ifdef W3_TRNK
+           JGS_TRUNK_DIGITS,                           &
+#endif
            JGS_SOURCE_NONLINEAR
       !
       WRITE (NDSO,2976)    P2SF, I1P2SF, I2P2SF,                    &
@@ -6342,7 +6349,7 @@ SUBROUTINE W3GRID()
 2922 FORMAT ( '  &SNL1 LAMBDA =',F7.3,', NLPROP =',E10.3,       &
          ', KDCONV =',F7.3,', KDMIN =',F7.3,','/           &
          '        SNLCS1 =',F7.3,', SNLCS2 =',F7.3,        &
-         ', SNLCS3 = ',F7.3,','/                            &
+         ', SNLCS3 = ',F7.3,','/                           &
          '        IQTYPE =',I2,', TAILNL =',F5.1,','/      &
          '        GQMNF1 =',I2,', GQMNT1 =',I2,',',        &
          ' GQMNQ_OM2 =',I2,', GQMTHRSAT =',E11.4,', GQMTHRCOU =',F4.3,','/ &
@@ -6714,7 +6721,12 @@ SUBROUTINE W3GRID()
          ',  JGS_DIFF_THR=', F8.3,                              &
          ',  JGS_NORM_THR=', F8.3,                              &
          ',  JGS_NLEVEL=', I3,                                  &
+#ifdef W3_TRNK
+         ',  JGS_TRUNK_DIGITS=', I3,                            &
+#endif
          ',  JGS_SOURCE_NONLINEAR=', L3 / )
+
+
     !
 960 FORMAT (/'  Miscellaneous ',A/                                   &
          ' --------------------------------------------------')

model/src/w3iogrmd.F90

@@ -805,6 +805,9 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT &
              B_JGS_DIFF_THR,                                  &
              B_JGS_NORM_THR,                                  &
              B_JGS_NLEVEL,                                    &
+#ifdef W3_TRNK
+             B_JGS_TRUNK_DIGITS,                              & 
+#endif
              B_JGS_SOURCE_NONLINEAR
 #ifdef W3_ASCII
         WRITE (NDSA,*)                                        &
@@ -839,6 +842,9 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT &
              B_JGS_DIFF_THR,                                  &
              B_JGS_NORM_THR,                                  &
              B_JGS_NLEVEL,                                    &
+#ifdef W3_TRNK
+             B_JGS_TRUNK_DIGITS,                              & 
+#endif
              B_JGS_SOURCE_NONLINEAR
 #endif
         !Init COUNTCON and IOBDP to zero, it needs to be set somewhere or
@@ -1004,6 +1010,9 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT &
              B_JGS_DIFF_THR,                                  &
              B_JGS_NORM_THR,                                  &
              B_JGS_NLEVEL,                                    &
+#ifdef W3_TRNK
+             B_JGS_TRUNK_DIGITS,                              & 
+#endif
              B_JGS_SOURCE_NONLINEAR
         IF (IERR.NE.0) CALL EXTIOF(NDSE,IERR,'W3IOGR','mod_def.'//FILEXT(:IEXT),51)
         IF (.NOT. GUGINIT) THEN

model/src/w3profsmd_pdlib.F90

@@ -5235,6 +5235,7 @@ SUBROUTINE APPLY_BOUNDARY_CONDITION(IMOD)
     REAL    :: RD1, RD2, RD10, RD20
     INTEGER :: IK, ITH, ISEA
     INTEGER :: IBI, IP_glob, ISP, JX
+
 #ifdef W3_S
     CALL STRACE (IENT, 'APPLY_BOUNDARY_CONDITION')
 #endif
@@ -5538,6 +5539,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
     USE W3PARALL, only : ListISPprevDir, ListISPnextDir
     USE W3PARALL, only : JX_TO_JSEA
     USE W3GDATMD, only: B_JGS_NLEVEL, B_JGS_SOURCE_NONLINEAR
+
     USE yowfunction, only : pdlib_abort
     USE yowNodepool, only: np_global
     USE W3DISPMD, only : WAVNU_LOCAL
@@ -5549,6 +5551,10 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
     USE W3GDATMD, only: REFPARS
 #endif
 
+#ifdef W3_TRNK
+    USE W3GDATMD, only: B_JGS_TRUNK_DIGITS
+#endif
+
     implicit none
     LOGICAL, INTENT(IN) :: LCALC
     INTEGER, INTENT(IN) :: IMOD
@@ -5603,6 +5609,11 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
     INTEGER JP_glob
     INTEGER is_converged, itmp
 
+#ifdef W3_TRNK
+    integer :: expVA
+    real    :: trVA
+#endif
+
     INTEGER :: TESTNODE = 923
 
     LOGICAL :: LSIG = .FALSE.
@@ -6327,6 +6338,24 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
       endif
     endif
 #endif
+
+#ifdef W3_TRNK
+    !
+    ! Truncate precision to B_JGS_TRUNK_DIGITS to enforce bit for bit reproducability
+    ! across repeated wavewatch runs.
+    !
+    DO IP = 1, npa
+      DO ISP=1,NSPEC
+        if (VA(ISP,IP) .gt. tiny(1.0) )then
+          expVA=nint(log10( VA(ISP,IP) ) )
+          trVA = 10.**( expVA - B_JGS_TRUNK_DIGITS )
+          if (trVA .gt. tiny(1.0)) VA(ISP,IP) = ANINT(VA(ISP,IP) / trVA) * trVA
+        endif
+      ENDDO
+    ENDDO
+#endif
+
+
     !
     call print_memcheck(memunit, 'memcheck_____:'//' WW3_PROP SECTION LOOP 7')
     !

regtests/bin/matrix.base

@@ -2175,6 +2175,8 @@
      echo "$rtst -s MPI -s PDLIB -i input_unstr -w work_unstr_b -g b -f -p $mpi -n $np $ww3 ww3_ufs1.1"   >> matrix.body
  # Domain Decomposition Implicit
      echo "$rtst -s MPI -s PDLIB -i input_unstr -w work_unstr_c -g c -f -p $mpi -n $np $ww3 ww3_ufs1.1"   >> matrix.body
+ # Domain Decomposition Implicit with precision truncation to force bit for bit reproducability
+     echo "$rtst -s MPI -s PDLIB_TRNK -i input_unstr -w work_unstr_d -g d -f -p $mpi -n $np $ww3 ww3_ufs1.1"  >>matrix.body  
   fi
 
  #Test of UFS applications with ww3_multi_esmf and grib2 output

regtests/ww3_ufs1.1/input_unstr/switch_PDLIB_TRNK

@@ -0,0 +1 @@
+TRNK NCO PDLIB SCOTCH NOGRB BIN2NC DIST MPI PR3 UQ FLX0 SEED ST4 STAB0 NL1 BT1 DB1 MLIM FLD2 TR0 BS0 RWND WNX1 WNT1 CRX1 CRT1 O0 O1 O2 O3 O4 O5 O6 O7 O14 O15 IC0 IS0 REF0

regtests/ww3_ufs1.1/input_unstr/ww3_grid_d.inp

@@ -0,0 +1,160 @@
+$ -------------------------------------------------------------------- $
+$ WAVEWATCH III Grid preprocessor input file                           $
+$ -------------------------------------------------------------------- $
+$
+  'Global Unstructured'
+$
+$ Frequency increment factor and first frequency (Hz) ---------------- $
+$ number of frequencies (wavenumbers) and directions, relative offset
+$ of first direction in terms of the directional increment [-0.5,0.5].
+$ In versions 1.18 and 2.22 of the model this value was by definiton 0,
+$ it is added to mitigate the GSE for a first order scheme. Note that
+$ this factor is IGNORED in the print plots in ww3_outp.
+$
+$gfs:  1.07  0.035  50  36  0.5
+$gefs: 1.1 0.035  33 36 0.5
+ 1.1 0.035  33 36 0.5
+$
+$ Set model flags ---------------------------------------------------- $
+$  - FLDRY         Dry run (input/output only, no calculation).
+$  - FLCX, FLCY    Activate X and Y component of propagation.
+$  - FLCTH, FLCK   Activate direction and wavenumber shifts.
+$  - FLSOU         Activate source terms.
+$
+   F T T T T T
+$
+$ Set time steps ----------------------------------------------------- $
+$ - Time step information (this information is always read)
+$     maximum global time step, maximum CFL time step for x-y and
+$     k-theta, minimum source term time step (all in seconds).
+$
+$ 
+  900. 900. 900. 900.
+$
+$ Start of namelist input section ------------------------------------ $
+$
+&OUTS 
+  USSP = 1, 
+  IUSSP = 3, 
+  STK_WN = 0.04, 0.110, 0.3305 /
+&SIN4 
+  BETAMAX = 1.315,
+  TAUWSHELTER = 1.0,
+  SWELLF = 0.798,
+  SWELLF2 = -0.0127,
+  SWELLF3 = 0.0151,
+  SWELLF4 = 100025.0,
+  SWELLF5 = 1.1999,
+  SWELLF7 = 235500.0 /
+&SNL1
+  NLPROP = 2.502E7 /
+&SDS4
+  FXFM3 = 2.501,
+  SDSC2 = -2.1975e-05,
+  SDSCUM = -0.4032,
+  SDSC6 = 0.2978,
+  SDSBR = 0.0009035 /
+&MISC 
+  CICE0 = 0.75,
+  CICEN = 0.75 /
+&SBT1 
+  GAMMA = -0.038 /
+&PRO3 
+  WDTHCG=1.5,
+  WDTHTH=1.5 /
+&UNST 
+  UGOBCAUTO = F,
+  UGOBCDEPTH = -10.,
+  EXPFSN = F,
+  EXPFSPSI = F,
+  EXPFSFCT = F,
+  IMPFSN = F,
+  EXPTOTAL = F,
+  IMPTOTAL = T,
+  IMPREFRACTION = T,
+  IMPFREQSHIFT = T,
+  IMPSOURCE = T,
+  SETUP_APPLY_WLV = F,
+  SOLVERTHR_SETUP=1E-14,
+  CRIT_DEP_SETUP=0.1,
+  JGS_NLEVEL = 0,
+  JGS_USE_JACOBI = T, 
+  JGS_BLOCK_GAUSS_SEIDEL = F,
+  JGS_TERMINATE_MAXITER = T,
+  JGS_MAXITER = 1000, 
+  JGS_TERMINATE_NORM = F,
+  JGS_TERMINATE_DIFFERENCE = T,
+  JGS_DIFF_THR = 1.E-6,
+  JGS_PMIN = 3.0,
+  JGS_LIMITER = F,
+  JGS_TRUNK_DIGITS = 5,
+  JGS_NORM_THR = 1.E-6 /
+$
+$ Mandatory string to identify end of namelist input section.
+$
+END OF NAMELISTS
+$
+$ FLAG for grid features
+$ 1 Type of grid 'UNST' 'RECT' 'CURV'
+$ 2 Flag for geographical coordinates (LLG)
+$ 3 Flag for periodic grid
+$
+$ Define grid -------------------------------------------------------- $
+$ Four records containing :
+$  1 NX, NY. As the outer grid lines are always defined as land
+$    points, the minimum size is 3x3.
+$  2 Grid increments SX, SY (degr.or m) and scaling (division) factor.
+$    If NX*SX is 360., latitudinal closure is applied.
+$  3 Coordinates of (1,1) (degr.) and scaling (division) factor.
+$  4 Limiting bottom depth (m) to discriminate between land and sea
+$    points, minimum water depth (m) as allowed in model, unit number
+$    of file with bottom depths, scale factor for bottom depths (mult.),
+$    IDLA, IDFM, format for formatted read, FROM and filename.
+$      IDLA : Layout indicator :
+$                  1   : Read line-by-line bottom to top.
+$                  2   : Like 1, single read statement.
+$                  3   : Read line-by-line top to bottom.
+$                  4   : Like 3, single read statement.
+$      IDFM : format indicator :
+$                  1   : Free format.
+$                  2   : Fixed format with above format descriptor.
+$                  3   : Unformatted.
+$      FROM : file type parameter
+$               'UNIT' : open file by unit number only.
+$               'NAME' : open file by name and assign to unit.
+$
+$ Example for longitude-latitude grid (switch !/LLG), for Cartesian
+$ grid the unit is meters (NOT km).
+$
+$
+ 'UNST' T T
+$
+   4.0 0.30  20  -1. 4 1 '(20f10.2)'  'NAME' '../input_unstr/global_1deg_unstr.msh'
+$
+$
+   10 3 1 '(....)' 'PART' 'mapsta.inp'
+$ Input boundary points ---------------------------------------------- $
+$  Close list by defining point (0,0) (mandatory)
+$
+      0   0   F
+$
+$ Excluded grid points from segment data ( FROM != PART )
+$   First defined as lines, identical to the definition of the input
+$   boundary points, and closed the same way.
+$
+      0   0   F
+$
+$   Second, define a point in a closed body of sea points to remove
+$   the entire body os sea points. Also close by point (0,0)
+$
+      0   0
+$
+$ Output boundary points --------------------------------------------- $
+$
+$  Close list by defining line with 0 points (mandatory)
+$
+       0.    0.    0.    0.       0
+$
+$ -------------------------------------------------------------------- $
+$ End of input file                                                    $
+$ -------------------------------------------------------------------- $