Electromechanics results diverge in parallel when using GMRES

[javijv4]

Description

After several time steps, the electromechanical tests (implemented in PR #596) using GMRES for the linear solver show noticeable differences between 3 and 4 cores. The differences are significantly smaller when using CG instead.

4-proc (GMRES)

3-proc (GMRES)

=== GMRES: 3 vs 4 ranks at t=50 (max abs diff) ===

model	dVm	dTact	dCa	dDisp
Regazzoni	5.0823e+00	3.1990e-04	1.0133e-04	8.5070e-04
NashPanfilov	5.0823e+00	1.0954e-02	1.0133e-04	3.4885e-02

=== CG: 3 vs 4 ranks at t=50 (max abs diff) ===

model	dVm	dTact	dCa	dDisp
Regazzoni	2.4833e-06	1.5796e-11	8.5102e-11	1.8047e-11
NashPanfilov	2.4833e-06	1.5706e-09	8.5102e-11	1.9004e-09

Reproduction

Run tests/cases/electromechanics/slab/solver_Regazzoni.xml or tests/cases/electromechanics/slab/solver_NashPanfilov.xml (in branch https://github.com/samibismar/svMultiPhysics/tree/feature-rdq20-mf, about to be merged) for several timesteps (above results are in timestep 50) using 3 or 4 cores.

Expected behavior

Serial and parallel runs with different cores should match.

Additional context

No response

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[michelebucelli]

Thanks for these tests @javijv4! One possible interpretation of this is that there is something that is particularly sensitive to changes in the order of operations, and this something is related to the GMRES implementation.

I have found some other clues that this is happening in previous issues/PRs, and I'm reporting them below.

1. Changing the order of some operations breaks some tests (Object oriented refactoring of fcType #577).

   Details

   I tried changing an expression from x = x + a + b to x += a + b. This changed the order of operations (in the new code expression a + b would be evaluated first, while in the old code x + a would be evaluated first). This broke the fluid/carreau_yasuda test, only with 4 processors, due to mismatch with the reference solution.

   Because of this, I had reverted that change.

   Original thread

2. Fast math optimization breaks all tests (Fixes for Intel OneAPI compilers and Intel MPI #584).

   Details

   The Intel OneAPI compilers turn on fast math optimization by default. This is a flag that allows the compiler to change the order of operations (among other things) in a way that is not IEEE-standard-compliant (that is, it can change the order of operations even if this changes the finite arithmetic result). This broke all tests in serial, meaning that all tests that compare to a reference solution would return a result different from that reference solution.

   I addressed this by disabling fast math optimization through the CMake configuration (this is a sensible decision regardless, in my opinion, but this behavior still points to something being very sensitive to finite arithmetic errors).

   Original thread

3. Sensitivity to Krylov space dimension.

   Details

   This I have investigated and documented less systematically, but: in an idealized ventricular electromechanics case that I've been working on as part of Coupling to 0D solver only supported for first equation #592 (solver.xml), the convergence behavior of the nonlinear solver for the structural equation is quite strongly dependent on the dimension of the Krylov subspace used in the linear solver.

   Examples of this, for the first time step:

   • Krylov dimension 50, 1 MPI process:

     ST 1-1  3.505e-01  [0 1.000e+00 1.000e+00 8.623e-14]  [53 -137 8]
     ST 1-2  5.966e-01  [-147 4.439e-08 4.439e-08 2.431e-05]  [17 -106 2]
     

   • Krlyov dimension 100, 1 MPI process:

     ST 1-1  3.958e-01  [0 1.000e+00 1.000e+00 7.746e-10]  [96 -210 23]
     ST 1-2  7.897e-01  [-28 3.825e-02 3.825e-02 7.178e-10]  [199 -19 38]
     ST 1-3  1.044e+00  [-78 1.173e-04 1.173e-04 1.072e-08]  [32 -184 4]
     

   • Krylov dimension 200, 1 MPI process:

     ST 1-1  4.008e-01  [0 1.000e+00 1.000e+00 7.746e-10]  [96 -210 22]
     ST 1-2  7.279e-01  [-28 3.825e-02 3.825e-02 8.584e-10]  [106 -209 25]
     ST 1-3  1.335e+00  [-43 6.437e-03 6.437e-03 9.929e-10]  [303 -48 60]
     ST 1-4  1.658e+00  [-47 4.114e-03 4.114e-03 8.932e-10]  [103 -208 25]
     ST 1-5  2.184e+00  [-70 3.138e-04 3.138e-04 1.931e-15]  [203 -167 54]
     

   • Krylov dimension 400, 1 MPI process:

     ST 1-1  3.991e-01  [0 1.000e+00 1.000e+00 7.746e-10]  [96 -210 23]
     ST 1-2  7.278e-01  [-28 3.825e-02 3.825e-02 8.584e-10]  [106 -209 26]
     ST 1-3  1.299e+00  [-43 6.437e-03 6.437e-03 8.793e-10]  [217 -209 57]
     ST 1-4  1.608e+00  [-11 2.532e-01 2.532e-01 8.019e-10]  [90 -209 20]
     ST 1-5  2.009e+00  [-41 8.679e-03 8.679e-03 8.911e-10]  [148 -208 39]
     ST 1-6  2.269e+00  [2 1.311e+00 1.311e+00 8.665e-10]  [37 -209 5]
     ST 1-7  2.535e+00  [-126 4.796e-07 4.796e-07 1.874e-06]  [23 -132 2]
     

   • higher Krlyov dimensions behave as for dimension 400 (at least with 1 MPI process), most likely because the max Krylov dimension is never reached from here on.
   • Krylov dimension 400, 2 MPI processes:

     ST 1-1  1.879e-01  [0 1.000e+00 1.000e+00 9.057e-10]  [36 -208 4]
     ST 1-2  3.256e-01  [-144 5.958e-08 5.958e-08 1.204e-05]  [19 -113 1]
     

   • Krylov dimension 400, 3 MPI processes:

      ST 1-1  1.408e-01  [0 1.000e+00 1.000e+00 8.237e-10]  [36 -209 5]
      ST 1-2  2.406e-01  [-144 5.661e-08 5.661e-08 2.268e-05]  [18 -107 2]
     

   • Krylov dimension 400, 4 MPI processes:

      ST 1-1  4.583e-01  [0 1.000e+00 1.000e+00 8.811e-10]  [397 -208 83]
      ST 1-2  5.529e-01  [1 1.255e+00 1.255e+00 7.253e-10]  [97 -210 23]
      ST 1-3  6.285e-01  [-19 1.054e-01 1.054e-01 9.340e-10]  [37 -208 5]
      ST 1-4  7.085e-01  [-173 2.225e-09 2.225e-09 3.320e-04]  [15 -80 1]
     

   • Kyrlov dimension 400, 6 MPI processes:

      ST 1-1  1.500e-01  [0 1.000e+00 1.000e+00 9.402e-10]  [133 -208 44]
      ST 1-2  4.411e-01  [-27 4.272e-02 4.272e-02 8.238e-10]  [291 -209 79]
      ST 1-3  5.104e-01  [1 1.163e+00 1.163e+00 8.939e-10]  [31 -208 7]
      ST 1-4  5.846e-01  [-133 2.234e-07 2.234e-07 4.219e-06]  [22 -124 4]
     

   Notice that the number of nonlinear iteration changes quite a lot (2 -> 7), and the convergence history is wildly inconsistent even for setups with the same number of nonlinear iterations. Nonlinear convergence is definitely not quadratic, as one would expect for Newton iterations.

   To my understanding of the interplay of linear and nonlinear solver this should not be happening: the Krylov space dimension should mostly affect the number of iterations (and memory occupation) the linear solver needs for convergence, but it should always converge to the same solution (up to tolerance), which is the only thing the nonlinear solver should see.

   One possible interpretation of this behavior is that the linear solver is converging to different things with different Krylov space dimensions, which might in turn suggest the accumulation of round-off errors in the linear solver, or early loss of orthogonality.

   (This behavior cannot be reproduced with CG, because CG does not use the Krylov space dimension parameter. I did however verify that the convergence history with CG is independent as the number of MPI processes varies between 1 and 8)

When encountering each of these issues, I have asked Claude to scan the codebase for areas potentially sensitive to round-off errors, and each time it has pointed me to the GMRES implementation (I think I mentioned this in one of the threads linked above), but I had never investigated this in more depth.

In light of @javijv4's findings, I think it is now appropriate to take a closer look at it. I'll post here if I find something useful.

[michelebucelli]

After investigation of the GMRES source code, I found that the current implementation is using unmodified Gram-Schmidth, which is known to be numerically unstable. This is likely leading to amplification of round-off errors, which would explain to some extent all of the issues that we are observing.

I am in the process of writing a fix that uses modified Gram-Schmidt, then I'll run the test that @javijv4 ran and see if there still are observable differences. I can already confirm that modified Gram-Schmidt addresses issues 1 and 3 from my comment above, so it's promising.

One possible reason why the classical Gram-Schmidt might be preferable to the modified Gram-Schmidt is that the former allows more compact parallel communication (one vector-valued MPI reduce for all the h(j,i) coefficients vs. one separate scalar reduce for every one of them), which might have a (probably small) impact on performance. Nonetheless, I think that the modified version produces more correct results and is the right choice in spite of this tradeoff.

A similar issue was addressed in deal.II, where they once had only modified Gram-Schmidt but then introduced the option for the user to switch between the two. I don't think it is worth the effort for svmp, but it's still worth mentioning.

[michelebucelli]

I have pushed the changes I made to this branch.

The most important difference is the use of modified Gram-Schmidt (which I refactored into its own helper function, which I think helps readability a bit and slightly reduces code duplication).

I have re-run the test that @javijv4 had run (only for Nash-Panfilov, since I'm basing this off of the main branch, with the following results:

=== 3 vs 4 ranks at t=50 (max abs diff) ===

CEP solver	struct solver	dVm	dTact	dCa	dDisp
GMRES	GMRES	1.2e-5	7.4e-9	4.0e-10	9.0e-9
GMRES	CG	1.2e-5	7.4e-9	4.0e-10	9.0e-9
CG	GMRES	3.9e-6	2.5e-9	1.3e-10	3.0e-9
CG	CG	3.9e-6	2.5e-9	1.3e-10	3.0e-9

I'm not sure why I'm getting slightly different vaules than @javijv4 for the CG/CG case, but now the discrepancy is small even with the GMRES solver (for some reason, pretty precisely three times as much as with CG).

It is also interesting (although maybe unimportant) that the mismatch for this test case appears to be dominated by what solver is used for the electrophysiology subproblem.

I am opening a PR shortly.

[michelebucelli]

I also observed that regression tests are all passing even after the change. I take it to mean one of two things:

1. most (if not all) tests only run for 1 timestep, so round-off errors do not have time to accumulate beyond the test tolerance, and for the changes to orthogonalization to become relevant;
2. the nonlinear solver stops when its residual falls below a certain tolerance, and this condition is independent of the accuracy of the linear solver; therefore, once the nonlinear solver converges, it always converges to the same solution (approximately), regardless of the linear solver behavior; at worst, the linear solver behavior should affect the convergence history.

I think both are happening at the same time: 1 is the probably why the mismatch between 3 and 4 cores was not caught by the automatic test, while I have personally observed 2.

In any case, there should be no need to regenerate the reference solutions, unless we want them to be byte-identical to the solution the code computes in its current state.

[ktbolt]

@michelebucelli Good work on this !

However, I think a PR is premature. The original was implemented using an unmodified Gram-Schmidth for a good reason so we need to better understand the implications of these changes. The changes work well for this toy CEP problem but we do need to see what affect they have on the run times of real million-element meshes for other types of problems, like FSI and ustruct.

[michelebucelli]

@ktbolt Do you have any insight about why the unmodified Gram-Schmidt might have been chosen? Or, put otherwise: do you know what that good reason was? As far as I know, its numerical instability problems are severe and well established, and there rarely is a good reason to use it, but I'd be interested in hearing otherwise.

The modified version is identical in exact arithmetic, and significantly cuts down errors. The only downside that I can see is that it requires some extra communication. This extra communication cost, by the way, might very well be offset by the improved convergence of both the linear and nonlinear solvers (which translates to less iterations). As far as I can tell, this is a case of quickly computing the wrong thing vs. computing the right thing, maybe a bit more slowly.

[ktbolt]

@michelebucelli Because much of the implementation details are obscure I do hesitate to modify core code that has been successfully used for several (15+) years.

I now recall that the original Fortran code that svFSI was based on used OpenMP for parallelism; the unmodified Gram-Schmidt inner products could be parallelized. The code was later modified to use MPI keeping the original sequence of steps I guess.

The changes you made are most likely the right thing to do but I would like to see if the extra communication cost is indeed offset by improved convergence.

[aabrown100-git]

@michelebucelli, could you try to analyze the output of the integration tests that use GMRES, before and after these changes? I guess we're interested in linear iterations and wall time?

[michelebucelli]

@ktbolt Thanks for the feedback. I understand your concern, but I feel that having visibly different results with different number of cores is a serious bug that, to me, should take priority over performance concerns (as of yet unverified).

Having said all this, as @aabrown100-git suggests, I think it is reasonable to analyze more in depth the performance (time and iterations) of the solver before and after this change for the integration tests. I will do that and report on iteration counts (linear and nonlinear) and wall time.

I will not be able to test this on a larger scale, but if someone has a larger simulation setup ready at hand and wanted to give it a try it wouldn't be a bad idea. It should be enough to switch to my branch and run for a few timesteps (enough to offset the setup time I'd say).

I now recall that the original Fortran code that svFSI was based on used OpenMP for parallelism; the unmodified Gram-Schmidt inner products could be parallelized. The code was later modified to use MPI keeping the original sequence of steps I guess.

Even in the current implementation the inner products are parallelized, in the sense that every process only computes the part of vector that it owns. The difference is that before all the process-local partial inner products would be stored in a vector, and then the parallel reduction operation would work on that vector as a single unit, whereas now each individual inner product needs to be reduced on its own, before moving on to the next inner product. The amount of data transferred between processes is the same, except that before it was transferred in a single call, while now it's transferred in multiple calls.

(I will also note, somewhat unrelatedly, that the code has leftover traces of OpenMP, e.g. in the name of the name space omp_la, although OpenMP is not used in the current implementation.)

For reference, I also looked into a few of the papers from Esmaily-Moghadam that discuss the linear solver and preconditioner, but I didn't find any mention of a deliberate choice of orthogonalization strategy. Not sure if this might be documented elsewhere.

[sujaldave]

Thanks @michelebucelli for the detailed investigation. I can help test this on a more representative cardiovascular problem.

As a starting point, I have an AAA CFD case with approximately 100k elements that I can run locally (also on HPC later for more time steps). I am planning to test it using 4, 8, and 10 MPI ranks and compare three configurations:

1. the current native GMRES,
2. GMRES with the proposed MGS change, and
3. Trilinos/Belos GMRES.

For the initial comparison, I'll try to keep the Krylov-space dimension and other solver settings the same/comparable across the three configurations and look at the convergence behavior, linear iteration counts, and run times as the number of MPI ranks changes.

As a follow-up, I can vary the Krylov-space dimension and see how sensitive each implementation is to that parameter.

I'll report the results here once I have the initial comparisons.

[alisonmarsden]

Thanks @sujaldave and @michelebucelli for digging into this! I agree that further testing is needed before we proceed. Also feel free to reach out to Mahdi for input. I will send an email to introduce him to you both.

[michelebucelli]

Thank you @alisonmarsden! And thank you @sujaldave for looking into this.

Here are some results from running the test suite on TACC Lonestar6. The following plots compare wall times, total non-linear iterations and total linear iterations with the following three orthogonalization versions:

1. classical Gram-Schmidt (blue dots, CGS), from the main branch;
2. modified Gram-Schmidt (red squares, MGS), from this branch;
3. low-synchronization MGS (green triangles, LS-MGS), from this experimental branch: this is a variant of MGS that uses a single MPI_Allreduce call, taken from this paper, at the price of a somewhat more complicated logic.

Both MGS and LS-MGS give coherent results for the electromechanics slab benchmark as the number of ranks changes, while CGS has the less stable behavior observed by @javijv4.

Plots

Key takeaways

1. MGS is sensibly slower than CGS for more than one rank; this is consistent with the extra reduction calls (one vector-valued reduction vs. many scalar-valued reductions transfer the same amount of data, but the MPI overhead is larger in the latter);
2. for the test suite, this increase is quite large (up to 76% with 4 cores);
3. there are a few test cases for which MGS reduces the total number of linear iterations, but also many test cases in which there's no changes;
4. for the test suite, LS-MGS is only slightly slower than CGS (around 5%, depending on the number of ranks), while maintaining the numerical stability of MGS (and its occasional reductions in iteration count);

Possible confounders

1. These results might be confounded by the test suite using generally small meshes; the tests run by @sujaldave will give better quantitative indication of this.
2. The tests are running with the same max Krylov space dimension (which defaults to 50, and is occasionally increased to 100). The more stable variants (MGS and LS-MGS) allow for a higher Krylov space size, and should that parameter be increased the improved convergence might partially offset the overhead. I have not yet looked into this aspect in depth, but I don't expect this will offset the 75% extra cost of MGS.

Should the additional tests confirm these results, I think that the small overhead of LS-MGS would be an acceptable tradeoff for its increased stability. At the same time, I think this is a judgement call and it might be a good idea to leave it to the user. Accordingly, I think it might be a good compromise to introduce a new parameter in <LS> to control the orthogonalization strategy, offering the user the option to switch between these three options. The default could remain CGS, for backward compatibility.

@sujaldave would it be possible for you to test the LS-MGS variant as well?

[sujaldave]

Thanks @michelebucelli . Of course, I can include the LS-MGS experimental branch in these comparisons, and I have done so here. I also included Trilinos-based runs. For the NS/bipartition solver, the Trilinos results are from my local experimental work shifting the bipartition solver to Trilinos with the resistance preconditioner.

I ran a rigid-wall CFD test on an abdominal aortic aneurysm (AAA) case with - NumberOfPoints="98523" and NumberOfCells="570927"

I ran 5 time steps for this initial comparison. Local runs were performed on a MacBook M4 Pro, and Sherlock runs were performed on Intel 5118 processors. I tested 4, 8, and 10 MPI ranks. I ran both the pure GMRES solver and the bipartition solver, labeled (NS) in the plots.

The GMRES/fsils settings were:

<LS type="GMRES" >
  <Linear_algebra type="fsils" >
    <Preconditioner> fsils </Preconditioner>
  </Linear_algebra>
  <Max_iterations> 5000 </Max_iterations>
  <Tolerance> 1e-4 </Tolerance>
  <Krylov_space_dimension> 300 </Krylov_space_dimension>
</LS>

The GMRES/Trilinos settings were:

<LS type="GMRES" >
  <Linear_algebra type="trilinos" >
    <Preconditioner> trilinos-diagonal </Preconditioner>
  </Linear_algebra>
  <Max_iterations> 1500000 </Max_iterations>
  <Tolerance> 1e-4 </Tolerance>
  <Krylov_space_dimension> 300 </Krylov_space_dimension>
</LS>

The NS/fsils settings were:

<LS type="NS" >
  <Linear_algebra type="fsils" >
    <Preconditioner> fsils </Preconditioner>
  </Linear_algebra>
  <Max_iterations> 10 </Max_iterations>
  <NS_GM_max_iterations> 200 </NS_GM_max_iterations>
  <NS_CG_max_iterations> 500 </NS_CG_max_iterations>
  <Tolerance> 0.4 </Tolerance>
  <NS_GM_tolerance> 0.01 </NS_GM_tolerance>
  <NS_CG_tolerance> 0.2 </NS_CG_tolerance>
  <Krylov_space_dimension> 300 </Krylov_space_dimension>
</LS>

The NS/Trilinos settings were:

<LS type="NS" >
  <Linear_algebra type="trilinos" >
    <GMRES_Preconditioner> trilinos-resistance </GMRES_Preconditioner>
    <CG_Preconditioner> trilinos-diagonal </CG_Preconditioner>
  </Linear_algebra>
  <Max_iterations> 10 </Max_iterations>
  <NS_GM_max_iterations> 200 </NS_GM_max_iterations>
  <NS_CG_max_iterations> 500 </NS_CG_max_iterations>
  <Tolerance> 0.4 </Tolerance>
  <NS_GM_tolerance> 0.01 </NS_GM_tolerance>
  <NS_CG_tolerance> 0.02 </NS_CG_tolerance>
  <Krylov_space_dimension> 300 </Krylov_space_dimension>
</LS>

I will attach the plots below. I am not diving too deeply into every detail here to keep the comment short, but there is definitely more information in these plots than what I summarize below.

plots

Key takeaways from the 5-step CFD runs:

• For cumulative wall time, the MGS and LS-MGS implementations are as good as or better than the native/classical GMRES runs among the cases tested.
• The mean percentage of time spent in the linear solver is comparable across the in-house GMRES implementations.
• The different residual histories for the pure GMRES runs are essentially identical and overlap with each other.
• The residual histories for the bipartition/NS runs are also comparable and in the same ballpark across implementations.
• The linear iteration counts are almost the same across processor counts and GMRES implementation types.

I also computed relative velocity L2 errors at each saved time step. These are computed as

||u_case - u_reference||_2 / ||u_reference||_2

For the first velocity-error figure, the reference is the native 4-rank case, and I compare native 8-rank and native 10-rank results against it. For the second velocity-error figure, each method is compared against the matching native case at the same processor count.

The larger velocity differences for the NS cases are likely a reflection of the coarser NS settings used in the XML. Those settings can be useful because the NS/bipartition runs are significantly faster, almost an order of magnitude faster than the pure GMRES solves in these short tests.

Next, I am running longer CFD comparisons, 4000 time steps, to see how these differences behave after the initial transients. I am also planning to run the same type of comparison for the FSI version of this AAA test case and compare flow variables over longer times.

More comments are very welcome. If there are specific quantities that would help decide the way forward, I can compute and add those as well.

PS - One thing I would like to clarify is the difference in the Max_iteration for the inhouse FSILS case and the Trilinos case. The inhouse GMRES counts the Max_iteration as Max_iteration x Krylov_space_dimension, hence for the Trilinos case, which treats them as they're passed in, I have kept the Max_iteration as 5000 x 300.