SimVascular / SimVascular/svMultiPhysics
Slow linear solver convergence for struct using trilinos preconditioner
Open
@sujaldave is already working on this.
Since Dec 16, 2025.
bug
- Dominant language
- C++
- Stars
- 45
- Forks
- 60
- Avg merge
- 5d 23h
- Merged PRs (30d)
- 11
Description
Description
While running a struct simulation using GMRES with Trilinos preconditioners, I am seeing:
- Very slow convergence compared to the default FSILS preconditioner
- 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
Contributor guide
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
Assessment
This issue has not been assessed yet.