Slow linear solver convergence for struct using trilinos preconditioner

[dseyler]

Description

While running a struct simulation using GMRES with Trilinos preconditioners, I am seeing:

1. Very slow convergence compared to the default FSILS preconditioner
2. A runtime crash with MueLu, which throws the following exception:

'std::length_error'
  what():  vector::_M_fill_insert

Posted below are the convergence histories for several preconditioner options. For the preconditioners that fail to converge, I tried increasing the number of linear solver iterations to 50,000 and observed similar results.

FSILS:

---------------------------------------------------------------------
Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 ST 1-1  3.679e+01  [0 1.000e+00 1.000e+00 9.979e-04]  [883 -3 94]
 ST 1-2  8.305e+01  [-30 2.890e-02 2.890e-02 9.989e-04]  [1111 -3 95]
 ST 1-3s 1.634e+02  [-90 3.074e-05 3.074e-05 9.982e-04]  [2186 -2 97]

MueLu:

std::length_error'
  what():  vector::_M_fill_insert

RILUK0:

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 ST 1-1  2.418e+02  [0 1.000e+00 1.000e+00 7.052e-01]  !2500 -2 98!  WARNING: The linear system solution has not converged
 ST 1-2  4.794e+02  [-3 7.053e-01 7.053e-01 9.499e-01]  !2500 0 98!  WARNING: The linear system solution has not converged
 ST 1-3  5.047e+02  [-3 6.700e-01 6.700e-01 9.775e-01]  !2500 0 92!  WARNING: The linear system solution has not converged

ILUT:

 ST 1-1  4.705e+01  [0 1.000e+00 1.000e+00 7.444e-01]  !2500 -1 90!  WARNING: The linear system solution has not converged
 ST 1-2  9.424e+01  [-2 7.445e-01 7.445e-01 9.937e-01]  !2500 0 90!  WARNING: The linear system solution has not converged
 ST 1-3  1.412e+02  [-2 7.398e-01 7.398e-01 1.000e+00]  !2500 0 90!  WARNING: The linear system solution has not converged

blockjacobi:

 ST 1-1  1.417e+01  [0 1.000e+00 1.000e+00 9.975e-04]  [874 -30 69]
 ST 1-2  3.113e+01  [-30 2.890e-02 2.890e-02 9.986e-04]  [1100 -30 73]
 ST 1-3s 5.944e+01  [-90 3.074e-05 3.074e-05 9.981e-04]  [2164 -30 84]

Has anyone had success using Trilinos for struct problems recently? Do the preconditioner settings need to be tuned differently for structural mechanics problems than for fluid or FSI? I know @aabrown100-git mentioned he had used trilinos-ilut for struct problems in the past, so this may be a bug related to recent changes or user error on my part.

@sujaldave @dcodoni

Reproduction

I have been using the following test case with a ~500,000 element mesh

<?xml version="1.0" encoding="UTF-8" ?>
<svMultiPhysicsFile version="0.1">

<GeneralSimulationParameters>
  <Continue_previous_simulation> 0 </Continue_previous_simulation>
  <Number_of_spatial_dimensions> 3 </Number_of_spatial_dimensions> 
  <Number_of_time_steps> 3 </Number_of_time_steps>
  <Time_step_size> 0.01 </Time_step_size>
  <Spectral_radius_of_infinite_time_step> 0.50 </Spectral_radius_of_infinite_time_step> 
  <Searched_file_name_to_trigger_stop> STOP_SIM </Searched_file_name_to_trigger_stop> 

  <Save_results_to_VTK_format> 1 </Save_results_to_VTK_format> 
  <Name_prefix_of_saved_VTK_files> result </Name_prefix_of_saved_VTK_files> 
  <Save_results_in_folder> results_trilinos </Save_results_in_folder>
  <Increment_in_saving_VTK_files> 1 </Increment_in_saving_VTK_files>
  <Start_saving_after_time_step> 1 </Start_saving_after_time_step> 

  <Increment_in_saving_restart_files> 1000 </Increment_in_saving_restart_files> 
  <Convert_BIN_to_VTK_format> 0 </Convert_BIN_to_VTK_format> 

  <Verbose> 1 </Verbose> 
  <Warning> 1 </Warning> 
  <Debug> 0 </Debug> 

</GeneralSimulationParameters>


<Add_mesh name="msh" > 

  <Mesh_file_path> ./Meshes/in_vivo_002_1d5mm-mesh-complete/mesh-complete.mesh.vtu </Mesh_file_path>

  <Add_face name="base">
      <Face_file_path> ./Meshes/in_vivo_002_1d5mm-mesh-complete/mesh-surfaces/top.vtp </Face_file_path>
  </Add_face>

  <Add_face name="epi">
      <Face_file_path> ./Meshes/in_vivo_002_1d5mm-mesh-complete/mesh-surfaces/epi.vtp </Face_file_path>
  </Add_face>

  <Add_face name="endo_lv">
      <Face_file_path> ./Meshes/in_vivo_002_1d5mm-mesh-complete/mesh-surfaces/endo_lv.vtp </Face_file_path>
  </Add_face>

  <Add_face name="endo_rv">
      <Face_file_path> ./Meshes/in_vivo_002_1d5mm-mesh-complete/mesh-surfaces/endo_rv.vtp </Face_file_path>
  </Add_face>

  <Fiber_direction_file_path> ./Meshes/in_vivo_002_1d5mm-mesh-complete/fibersLong.vtu </Fiber_direction_file_path>
  <Fiber_direction_file_path> ./Meshes/in_vivo_002_1d5mm-mesh-complete/fibersSheet.vtu </Fiber_direction_file_path>

  <Domain_file_path> ./Meshes/in_vivo_002_1d5mm-mesh-complete/mesh-complete_DOMAIN_ID.dat </Domain_file_path>

  <Mesh_scale_factor> 0.001 </Mesh_scale_factor> <!-- Mesh is in mm. Convert to SI (m) -->

</Add_mesh>


<!-- Using SI units-->
<Add_equation type="struct" > 

   <Coupled> true </Coupled>
   <Min_iterations> 3 </Min_iterations>  
   <Max_iterations> 20 </Max_iterations> 
   <Tolerance> 1e-4 </Tolerance> 

   <!-- Myocardium (Domain 1) -->
   <Domain id="1">
      <Equation> struct </Equation>
      <Density> 1055.0 </Density>                           <!-- kg/m^3 -->
      <Elasticity_modulus> 100000 </Elasticity_modulus>   <!-- Pa. Used to compute material wave speeds for struct stab parameters -->
      <Poisson_ratio> 0.48333 </Poisson_ratio>

      <Viscosity model="Newtonian">
         <Value> 100 </Value> <!-- Pa.s -->
      </Viscosity>

      <Momentum_stabilization_coefficient> 1e-5 </Momentum_stabilization_coefficient>
      <Continuity_stabilization_coefficient> 0.0 </Continuity_stabilization_coefficient>

      <!-- Holzapfel-Ogden constitutive model (parameters optimized by iFEA) -->
      <Constitutive_model type="HolzapfelOgden">
         <a> 590 </a>                                   <!-- Pa -->
         <b> 8.023 </b>                                   <!-- no units -->
         <a4f> 184720 </a4f>                             <!-- Pa -->
         <b4f> 16.026 </b4f>                             <!-- no units -->
         <a4s> 0 </a4s>                             <!-- Pa (fixed at 0) -->
         <b4s> 11.12 </b4s>                             <!-- no units (fixed from Arostica et al. 2025) -->
         <afs> 0 </afs>                             <!-- Pa (fixed at 0) -->
         <bfs> 11.436 </bfs>                             <!-- no units (fixed from Arostica et al. 2025) -->
         <k> 100.0 </k>                                  <!-- Heaviside smoothing parameter for HO model -->
      </Constitutive_model>

      <Dilational_penalty_model> ST91 </Dilational_penalty_model>
      <Penalty_parameter> 1.0e6 </Penalty_parameter>

   </Domain>

   <!-- Valve rings (Domain 2) -->
   <Domain id="2">
      <Equation> struct </Equation>
      <Density> 1055.0 </Density>                           <!-- kg/m^3 -->
      <Elasticity_modulus> 100000.0 </Elasticity_modulus>   <!-- Pa -->
      <Poisson_ratio> 0.48333 </Poisson_ratio>

      <Viscosity model="Newtonian">
         <Value> 100.0 </Value> <!-- Pa.s -->
      </Viscosity>

      <Momentum_stabilization_coefficient> 1e-5 </Momentum_stabilization_coefficient>
      <Continuity_stabilization_coefficient> 1e-5 </Continuity_stabilization_coefficient>

      <!-- Isotropic and stiff for valve rings -->
      <Constitutive_model type="HolzapfelOgden">
         <a> 1.0e6 </a>        <!-- Pa -->
         <b> 5.0 </b>           <!-- no units -->
         <a4f> 0 </a4f>         <!-- Pa -->
         <b4f> 6.0 </b4f>       <!-- no units -->
         <a4s> 0 </a4s>         <!-- Pa -->
         <b4s> 11.1 </b4s>      <!-- no units -->
         <afs> 0 </afs>         <!-- Pa -->
         <bfs> 11.436 </bfs>    <!-- no units -->
         <k> 100.0 </k>         <!-- Heaviside smoothing parameter for HO model -->
      </Constitutive_model>

      <Dilational_penalty_model> ST91 </Dilational_penalty_model>
      <Penalty_parameter> 1.0e6 </Penalty_parameter>
   </Domain>

   <Output type="Spatial" >
      <Displacement> true </Displacement>
      <Velocity> true </Velocity>
      <Acceleration> true </Acceleration>
      <Jacobian> true </Jacobian>
      <Stress> true </Stress>
      <Strain> true </Strain>
      <Cauchy_stress> true </Cauchy_stress>
      <Def_grad> true </Def_grad>
      <VonMises_stress> true </VonMises_stress>
      <Fiber_shortening> true </Fiber_shortening>
   </Output>

   <LS type="GMRES" >
      <Linear_algebra type="trilinos" >
         <Preconditioner> trilinos-riluk0 </Preconditioner>
      </Linear_algebra> 
      <Tolerance> 1e-3 </Tolerance>
      <Max_iterations> 2500 </Max_iterations> 
      <Krylov_space_dimension> 100 </Krylov_space_dimension>
   </LS>

   <Output type="Boundary_integral" >
      <Velocity> true </Velocity>
   </Output>

   <!-- LV endocardial boundary condition: pressure from file -->
   <Add_BC name="endo_lv" > 
      <Type> Neu </Type> 
      <Time_dependence> Unsteady </Time_dependence> 
      <Temporal_values_file_path> ./0D_pressure_data/p_LV_0D_Pascal.dat </Temporal_values_file_path>
      <Ramp_function> false </Ramp_function>
      <Follower_pressure_load> true </Follower_pressure_load>
   </Add_BC>

   <!-- RV endocardial boundary condition: pressure from file -->
   <Add_BC name="endo_rv" > 
      <Type> Neu </Type> 
      <Time_dependence> Unsteady </Time_dependence> 
      <Temporal_values_file_path> ./0D_pressure_data/p_RV_0D_Pascal.dat </Temporal_values_file_path>
      <Ramp_function> false </Ramp_function>
      <Follower_pressure_load> true </Follower_pressure_load>
   </Add_BC>

   <Add_BC name="epi" > 
      <Type> Robin </Type> 
      <Time_dependence> Steady </Time_dependence> 
      <Stiffness> 1.0e6 </Stiffness>  <!-- Pa/m -->
      <Damping> 5.0e3 </Damping>      <!-- Pa.s/m -->
   </Add_BC>

   <!-- Base boundary condition: Robin BC -->
   <Add_BC name="base" > 
      <Type> Robin </Type> 
      <Time_dependence> Steady </Time_dependence> 
      <Stiffness> 1.0e5 </Stiffness>  <!-- Pa/m -->
      <Damping> 5.0e3 </Damping>      <!-- Pa.s/m -->
   </Add_BC>
   
</Add_equation>

</svMultiPhysicsFile>

Expected behavior

The linear solver should converge, ideally faster than FSILS when using Trilinos preconditioners.

Additional context

No response

Code of Conduct

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

[aabrown100-git]

Related to this is I think there are no struct or ustruct test cases that use Trilinos linear solvers. We should add a struct/ustruct Trilinos test case, probably just a modification of an existing case.

[dcodoni]

@dseyler not all preconditioners are good for all cases, they are problem specific. I see that blockjacobi works fine and it is faster than fsils.
Also it is important to remember that the max iterations for fsils means max number of restarts so 100x2500 total number of iterations. For trilinos you need to specify the actual number of iterations. Moreover, I noticed that usually the linear solver performance in Trilinos deteriorates for low number of Krylov iterations.
I would suggest:

• krylov number: 300
• max iterations: 5000 or more

Ifpack2 does not implement ILUT in the same sense that it was in the previous Ifpack (older version), so it is suggested to use riluk0 which was working well for ustruct, but I don't think I tested struct.

[dseyler]

@dcodoni Thanks for looking into this. I just tried using riluk0 with a krylov number of 300 and 10,000 max iterations for struct, but it still doesn't converge. I tried it previously with 50,000 max iterations (but 100 krylov dimensions) and still didn't see any convergence.

Here's the output for riluk0 with 300 krylov dimenstions and 10,000 iterations:

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 ST 1-1  7.570e+02  [0 1.000e+00 1.000e+00 4.234e-01]  !10000 -4 99!  WARNING: The linear system solution has not converged
 ST 1-2  9.168e+02  [-7 4.249e-01 4.249e-01 5.921e-01]  !10000 -2 99!  WARNING: The linear system solution has not converged
 ST 1-3  1.069e+03  [-11 2.516e-01 2.516e-01 9.429e-01]  !10000 0 99!  WARNING: The linear system solution has not converged

When I tried the same for ustruct, I received the following error:

terminate called after throwing an instance of 'std::runtime_error'
  what():  [svMultiPhysics] Equations with ustruct physics must use fsils for assembly.

Edit: after explicitly changing the asssembly type to fsils for ustruct, it produces a segfault:

[sh03-08n18:06360] *** Process received signal ***
[sh03-08n18:06360] Signal: Segmentation fault (11)
[sh03-08n18:06360] Signal code: Address not mapped (1)
[sh03-08n18:06360] Failing at address: (nil)
[sh03-08n18:06360] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x45330)[0x7f3c11fd7330]
[sh03-08n18:06360] [ 1] /home/groups/amarsden/svMultiPhysics/svMultiPhysics_aabrown100-git/build/svMultiPhysics-build/bin/svmultiphysics(+0x4f34d9)[0x55b7360844d9]
[sh03-08n18:06360] [ 2] /home/groups/amarsden/svMultiPhysics/svMultiPhysics_aabrown100-git/build/svMultiPhysics-build/bin/svmultiphysics(+0x4f3fad)[0x55b736084fad]
[sh03-08n18:06360] [ 3] /home/groups/amarsden/svMultiPhysics/svMultiPhysics_aabrown100-git/build/svMultiPhysics-build/bin/svmultiphysics(+0x5e50f7)[0x55b7361760f7]
[sh03-08n18:06360] [ 4] /home/groups/amarsden/svMultiPhysics/svMultiPhysics_aabrown100-git/build/svMultiPhysics-build/bin/svmultiphysics(+0x49b422)[0x55b73602c422]
[sh03-08n18:06360] [ 5] /lib/x86_64-linux-gnu/libc.so.6(+0x2a1ca)[0x7f3c11fbc1ca]
[sh03-08n18:06360] [ 6] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0x8b)[0x7f3c11fbc28b]
[sh03-08n18:06360] [ 7] /home/groups/amarsden/svMultiPhysics/svMultiPhysics_aabrown100-git/build/svMultiPhysics-build/bin/svmultiphysics(+0x4e1cc5)[0x55b736072cc5]
[sh03-08n18:06360] *** End of error message ***

Which case were you using to test riluk0 for ustruct?

Also, do you have any thoughts on why ML is failing to run entirely?