VP solver robustness issues ("bad departure points") (was: Verify b4b-ness of different MPI decompositions for the VP solver / performance evaluation of `repro-vp` branch)

[phil-blain]

Tried a new suite with dynpicard, no OpenMP

# Test         Grid    PEs        Sets    BFB-compare
smoke          gx3     1x1        diag1,dynpicard
sleep 180
smoke          gx3     8x1        diag1,dynpicard   smoke_gx3_1x1_diag1_dynpicard
smoke          gx3     2x1        diag1,dynpicard   smoke_gx3_1x1_diag1_dynpicard
smoke          gx3     4x1        diag1,dynpicard   smoke_gx3_1x1_diag1_dynpicard

none of the three MPI cases are bfb with the 1x1 case.

[phil-blain]

Tried setting b4bflag=reprosum in the 2x1 case. Still not b4b with the 1x1 case.

[phil-blain]

Tried with maxits_nonlin=500 Still not b4b with the 1x1 case.

[phil-blain]

OK, looking at the differences between 1x1 and {2x1, 4x1, 8x1} there are some differences, bigger in 4x1 and 8x1.

• try with more nonlinear iterations
• try disabling advection and thermo

[phil-blain]

With 5000 nonlinear iterations, differences are in the same range for 2, 4, 8 procs (vs 1). i.e. [-1E-6, 1E6].

[phil-blain]

note: this tests runs 1 day, and over these 24 time steps, only 4 need over 500 iterations to reach reltol_nonlin=1E-8 (default value in the namelist). With 5000 as maxits_nonlin, we still do not get more than max 1398 iterations.

[phil-blain]

With reltol_nonlin=1E-12 and maxits_nonlin=5000, the difference vs 1 proc is more or less the same for 2, 4, 8 procs, and looks something like this:

interestingly :

• there are more differences in the South hemisphere (likely due to the fact that the concentration is very close to 1 in all the North hemisphere
• the differences are not only near the coast, and do not appear to be near only the ice edge either

Note that we do not nearly reach 5000 iterations with this value of reltol_nonlin.

[phil-blain]

To get the number of iterations to reach the required tolerance:

\grep Arctic -B1 path_to_case_directory/logs/cice.runlog.220606-19* |\grep monitor

this is with diagfreq = 1 and monitor_nonlin=.true.

[phil-blain]

• If we disable thermo, ridging and transport, then with the same settings as above (reltol_nonlin=1E-12), we get differences in the range [-1E-16,1E-17] for uvel,vvel which is much more in line with what is expected. aice and hi are b4b.
• If we only disable thermo (and same settings again), then we get slightly higher differences in the veloc components:

    -1 :       Date     Time   Level Gridsize    Miss :     Minimum        Mean     Maximum : Parameter name
     9 : 2005-01-02 00:00:00       0    11600    3594 : -3.7470e-16  8.6176e-20  3.4478e-16 : aice     
    11 : 2005-01-02 00:00:00       0    11600    3594 : -4.4409e-16  7.4066e-20  4.4409e-16 : hi   
    13 : 2005-01-02 00:00:00       0    11600    3594 : -5.2082e-17  8.6478e-17  6.9237e-13 : uvel          
    14 : 2005-01-02 00:00:00       0    11600    3594 : -5.5376e-13 -6.9204e-17  1.0734e-17 : vvel    

and now aice and hi start to not be b4b, but are still small

[phil-blain]

• try the diag (and maybe ident) precond
• do a long (5 years) gx1 run with different MPI decomp and compare the thickness field at the end (mid-Jan).

[phil-blain]

With precond=diag, no thermo, no transport:

    -1 :       Date     Time   Level Gridsize    Miss :     Minimum        Mean     Maximum : Parameter name    
    13 : 2005-01-02 00:00:00       0    11600    3594 : -5.6379e-18  6.6466e-20  3.9514e-16 : uvel          
    14 : 2005-01-02 00:00:00       0    11600    3594 : -2.1605e-16 -5.4175e-21  2.0172e-16 : vvel 

[phil-blain]

With `precond='diag', with thermo, with transport, we get a model abort:

istep1:         2    idate:  20050101    sec:      7200
 (JRA55_data) reading forcing file 1st ts = /home/ords/cmdd/cmde/sice500//CICE_data/forcing/gx3/JRA55/8XDAILY/JRA55_gx3_03hr_forcing_2005.nc
Rank 2 [Mon Jun 13 21:08:20 2022] [c0-0c0s9n1] application called MPI_Abort(MPI_COMM_WORLD, 128) - process 2
    (icepack_warnings_setabort) T :file icepack_itd.F90 :line          900
 (cleanup_itd) aggregate ice area out of bounds
  (cleanup_itd)aice:   1.00245455360081
  (cleanup_itd)n, aicen:           1  0.676531837003209
  (cleanup_itd)n, aicen:           2  0.224493247425031
  (cleanup_itd)n, aicen:           3  4.769818624818129E-002
  (cleanup_itd)n, aicen:           4  3.766552529401467E-002
  (cleanup_itd)n, aicen:           5  1.606575763037559E-002
 (icepack_warnings_aborted) ... (icepack_step_therm2)

weird as I did his test two years ago (#33 (comment)), although with reltol_nonlin=1E-8...)

@JFLemieux73 je garde une trace de mes expériences MPI dans cette issue si tu veux rester au courant.

[phil-blain]

OK, it's because I forgot to also re-enable ridging, the model did not like that (is that expected?...)

EDIT after discussing with JF, yes it is expected, convergence can cause that.

[phil-blain]

OK, with ridging, advection, and transport, reltol=1E-12, precond diag;

    -1 :       Date     Time   Level Gridsize    Miss :     Minimum        Mean     Maximum : Parameter name
     9 : 2005-01-02 00:00:00       0    11600    3594 : -3.9378e-06 -4.8300e-10  4.4782e-06 : aice          
    11 : 2005-01-02 00:00:00       0    11600    3594 : -1.9175e-06 -7.3383e-11  4.1247e-06 : hi            
    13 : 2005-01-02 00:00:00       0    11600    3594 : -2.8908e-07  7.0716e-09  3.3297e-05 : uvel          
    14 : 2005-01-02 00:00:00       0    11600    3594 : -5.4352e-07  6.0686e-09  2.5828e-05 : vvel  

[phil-blain]

idem with pgmres:

phb001@xc4elogin1(daley): [17:06:02] $ cdo infov diff.nc 2>/dev/null | \grep -E 'name|aice|hi|vel'
    -1 :       Date     Time   Level Gridsize    Miss :     Minimum        Mean     Maximum : Parameter name
     9 : 2005-01-02 00:00:00       0    11600    3594 : -3.8860e-06  3.3100e-09  2.7253e-05 : aice          
    11 : 2005-01-02 00:00:00       0    11600    3594 : -9.5537e-06 -1.0926e-09  2.0109e-06 : hi            
    13 : 2005-01-02 00:00:00       0    11600    3594 : -4.9737e-07  5.9547e-10  2.8301e-06 : uvel          
    14 : 2005-01-02 00:00:00       0    11600    3594 : -2.5361e-06 -5.8569e-10  1.2014e-06 : vvel  

I don't think it's only the preconditoner since these results are similar as the 'diag' precond.

[phil-blain]

So I did side by side, step by step debugging of 1x1 vs 2x1. The values are the same on both side until the first normalization in the FGMRES algorithm. Since we do a global sum of different numbers, in a different order (that mathematically sum to the same result on all decompositions), then because of floating point arithmetic we get a different norm (of the residual), and then we propagate that to the whole of the vectors by normalizing.

So in the end it is not surprising that we get different results. We will run a QC test of different decompositions against each other to ensure we get the same climate.