Trilinos MueLu default number of equations is incorrect for struct

[dseyler]

Description

In the default Trilinos ML options, the "number of equations" is hardcoded to 4 (Trilinos_impl.cpp, line 767), causing struct simulations to throw an error. Users can currently avoid this by including a custom Trilinos options file, but I think a more permanent default solution would be preferable.

Changing the "number of equations" parameter value to the number of degrees of freedom fixes this problem:

mueluParams.set("number of equations", 4);   // dof

should be

mueluParams.set("number of equations", dof); 

This is a simple one-line fix, which I implemented and tested in this branch.

I think it would also be a good idea to include a struct Trilinos ML test case to catch this bug (and potentially others) in the future. The branch I linked includes one that is a variant of LV_HolzapfelOgden_active.

Reproduction

Running a sufficiently large struct simulation (above coarse: max_size = 2000) with the default Trilinos ML settings results in the following error:

std::length_error from vector::_M_fill_insert

Expected behavior

Struct simulations should run without error when using the default Trilinos ML settings.

Additional context

No response

Code of Conduct

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

[michelebucelli]

Thanks @dseyler! I would also be interested in trying out MueLu for both structural mechanics and electrophysiology, in the context of electromechanics.

I took a look at the code, and one thing that took me a moment to figure out is from where the value of dof is coming from. What I could understand is the following:

1. dof is a global variable defined in trilinos_impl.cpp;
2. its value is set by calling trilinos_lhs_create, which contains the assignment dof = Dof, with Dof being one of the function arguments;
3. trilinos_lhs_create is called from TrilinosLinearAlgebra::TrilinosImpl::alloc(), which sets Dof = com_mod.dof;
4. the documentation for ComMod::dof says Current equation degrees of freedom, but I could not figure out whether this variable is updated when iterating between different equations.

All this is to say:

1. can this be expected to work if I were to try and use MueLU for multiple equations within the same coupled problem?
2. if not, and I were to use it for one equation but not the other, is there any implicit assumption on how those equations need to be ordered?

[dseyler]

@michelebucelli Thanks for asking this, since it made me dig through the code a bit to verify that everything is working correctly.

1. dof is a global variable, but it is updated with each Newton iteration in Integrator::initiator() where the current equation's dof is copied into ComMod.
2. It then reaches Trilinos at TrilinosImpl::alloc as dof = com_mod.dof (still the dof of the current equation).
3. Then it is passed into trilinos_lhs_create as you already noted.
4. The Sparsity pattern and parallel distribution map is then built using whatever value of dof was passed into trilinos_lhs_create for that Newton iteration.

Since everything on the Trilinos side is rebuilt per equation, I believe this should work for multiple equations and ordering should not matter. It's still probably a good idea to include a coupled test case to verify this. Do you have one you could use to check it? As long as the mesh has more than ~2000 elements, it should trigger the error if there is one. Electro-mechanics would work since they have different dof's.

Related, I think we shouldn't recompute the Tpetra parallel distribution map and sparsity graph at each Newton iteration since these are fairly large problems and can be cached for each equation. Only the values of the matrix need to be recomputed at each iteration. For the structural mechanics simulation I've been working with, this accounts for 2-3% of the total runtime (10% of the linear solve).

[ktbolt]

@dseyler trilinos_impl.cpp has a bunch of variables declared at the file level so they have global namespace scope.

[dseyler]

@ktbolt Gotcha. Do you agree with how I traced dof through the code?

For the larger change regarding recomputing Tpetra, I think we'd have to move some of those global variables from the .ccp file into the Trilinos struct in Trilinos_impl.h

[ktbolt]

@dseyler Your trace of dof is correct I think.

The Trillinos interface was copied from the svFSI version; has since been modified but it still can only be used in a single equation.

The Trillinos interface probably does need to be refactored to support using it in multiple equations. It is not clear how well Trillinos might support this. @sujaldave Can you comment on this ?

[sujaldave]

I agree with @dseyler 's trace of dof.
For the full-system MueLu path, number of equations should be the current equation block size, i.e. com_mod.dof, not hard-coded to 4. The value 4 is only appropriate for a 3D fluid/FSI velocity-pressure system. For structural mechanics it should be 3, and for EP it should be 1.

Trilinos/MueLu can support these cases; the main limitation is in our current svMultiPhysics Trilinos adapter. It still has file-scope state in trilinos_impl.cpp, so it supports one active equation at a time. Since com_mod.dof is updated from the current equation in Integrator::initiator(), and TrilinosImpl::alloc() rebuilds the Tpetra map/graph using that current value, I would expect the current implementation to work for multiple equations as long as each equation allocates/rebuilds its Trilinos objects before solving. In that mode, equation ordering should not matter.

That said, I agree we should verify this with a coupled test where different equations use different dof, for example electromechanics with EP (dof=1) and mechanics (dof=3). I'll start working on this and create a PR soon.

[dseyler]

@sujaldave do the Trilinos objects need to be rebuilt before each solve or can they be computed once per equation and then reused at each subsequent iteration? I think this should be fine as long as there's no remeshing? For my test case, the Trilinos object rebuild takes longer than the linear solve itself

[sujaldave]

@dseyler I think the Tpetra map, importer, graph, and vector objects should be reusable once per equation as long as there is no remeshing and the equation’s dof, local/global node numbering, ownership, and sparsity pattern remain unchanged. The matrix and RHS values still need to be reset and reassembled each nonlinear iteration. I shall look further into it.

[dseyler]

Thanks @sujaldave! In the meantime, do you think the dof fix is safe to implement?

[sujaldave]

@dseyler Yes, based on my testing I think the dof fix is safe to implement.

In addition to the new LV_HolzapfelOgden_active_trilinos_ml case that you added, I tested this on an abdominal aortic aneurysm FSI case. The fluid mesh has NumberOfPoints="98523" and NumberOfCells="570927", and the solid mesh has NumberOfPoints="24938" and NumberOfCells="82852".

I compared the current main branch against your branch

The FSI equation uses trilinos-ml in all tests. I then tested two mesh-equation preconditioner choices -trilinos-diagonal and trilinos-ml for the mesh equation

For the trilinos-diagonal mesh solve, main and your branch produce matching nonlinear/linear iteration histories. For example, the FSI and mesh residual histories and iteration counts are essentially identical across the two runs, with only small timing differences. This suggests that changing MueLu’s "number of equations" from hard-coded 4 to the active equation dof does not perturb the FSI equation path, where dof = 4.

main branch

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 FS 1-1  5.258e+00  [0 1.000e+00 1.000e+00 8.781e-05]  [39 -41 32]
 MS 1-1  7.277e+00  [0 1.000e+00 1.000e+00 9.947e-05]  [230 -40 16]
 FS 1-2  1.244e+01  [-31 2.818e-02 2.818e-02 2.731e-05]  [38 -46 32]
 MS 1-2  1.444e+01  [-53 2.137e-03 2.137e-03 9.883e-05]  [282 -40 20]
 FS 1-3  1.934e+01  [-75 1.652e-04 1.652e-04 1.896e-05]  [34 -47 29]
 MS 1-3  2.141e+01  [-109 3.176e-06 3.176e-06 9.996e-05]  [310 -40 23]
 FS 1-4s 2.638e+01  [-119 1.061e-06 1.061e-06 1.710e-05]  [34 -48 28]

fix/muelu-number-of-equations-dof branch

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 FS 1-1  5.342e+00  [0 1.000e+00 1.000e+00 8.781e-05]  [39 -41 32]
 MS 1-1  7.301e+00  [0 1.000e+00 1.000e+00 9.947e-05]  [230 -40 19]
 FS 1-2  1.248e+01  [-31 2.818e-02 2.818e-02 2.731e-05]  [38 -46 31]
 MS 1-2  1.448e+01  [-53 2.137e-03 2.137e-03 9.883e-05]  [282 -40 20]
 FS 1-3  1.946e+01  [-75 1.652e-04 1.652e-04 1.896e-05]  [34 -47 27]
 MS 1-3  2.154e+01  [-109 3.176e-06 3.176e-06 9.996e-05]  [310 -40 23]
 FS 1-4s 2.662e+01  [-119 1.061e-06 1.061e-06 1.710e-05]  [34 -48 27]

For the trilinos-ml mesh solve, main fails immediately after the first FSI solve:

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
  FS 1-1  5.270e+00  [0 1.000e+00 1.000e+00 8.781e-05]  [39 -41 32]

With your branch, the same trilinos-ml mesh-equation case runs. The mesh CG solve shows poor convergence with the current default MueLu parameters, frequently reaching the 1000-iteration limit

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 FS 1-1  5.146e+00  [0 1.000e+00 1.000e+00 8.781e-05]  [39 -41 31]
 MS 1-1  2.457e+01  [0 1.000e+00 1.000e+00 2.758e-03]  !1000 -26 90!  WARNING: The linear system solution has not converged
 FS 1-2  2.963e+01  [-30 2.820e-02 2.820e-02 2.729e-05]  [38 -46 31]
 MS 1-2  4.910e+01  [-49 3.496e-03 3.496e-03 4.205e-02]  !1000 -14 90!  WARNING: The linear system solution has not converged
 FS 1-3  5.381e+01  [-65 5.334e-04 5.334e-04 2.422e-05]  [30 -46 25]
 MS 1-3  7.338e+01  [-76 1.471e-04 1.471e-04 6.685e-02]  !1000 -12 90!  WARNING: The linear system solution has not converged
 FS 1-4  7.786e+01  [-91 2.698e-05 2.698e-05 2.152e-05]  [28 -47 22]
 FS 1-4s 9.739e+01  [-91 2.698e-05 2.698e-05 2.152e-05]  [28 -47 5]

but this appears to be a separate tuning issue. The default MueLu parameters seem better suited to the monolithic FSI/GMRES system than to the mesh CG system. Importantly, the dof fix allows the mesh MueLu hierarchy to be constructed correctly instead of crashing.