Slow nonlinear convergence in ustruct/LV_Guccione_active

[mrp089]

Description

While adding missing results to test cases (#175), I noticed that the test ustruct/LV_Guccione_active has a very slow nonlinear convergence.

Reproduction

In ./tests, run pytest -vrP -k LV_Guccione_active or have a look at this pipeline.

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 ST 1-1  2.720e-01  [0 1.000e+00 1.000e+00 9.031e-13]  [96 -126 3]
 ST 1-2  5.520e-01  [-27 4.192e-02 4.192e-02 8.019e-13]  [95 -125 2]
 ST 1-3  8.300e-01  [-42 7.443e-03 7.443e-03 4.504e-12]  [88 -102 2]
 ST 1-4  1.108e+00  [-55 1.616e-03 1.616e-03 1.866e-11]  [82 -85 3]
 ST 1-5  1.399e+00  [-68 3.657e-04 3.657e-04 8.086e-11]  [77 -70 2]
 ST 1-6  1.678e+00  [-81 8.364e-05 8.364e-05 3.345e-10]  [72 -56 2]
 ST 1-7  1.968e+00  [-94 1.919e-05 1.919e-05 1.418e-09]  [67 -39 1]
 ST 1-8  2.254e+00  [-107 4.404e-06 4.404e-06 8.073e-09]  [61 -23 1]
 ST 1-9  2.529e+00  [-119 1.011e-06 1.011e-06 3.121e-08]  [56 -10 1]
 ST 1-10  2.793e+00  [-132 2.321e-07 2.321e-07 1.121e-07]  [50 -160 2]
 ST 1-11  3.063e+00  [-145 5.326e-08 5.326e-08 5.740e-07]  [45 -144 1]
 ST 1-12  3.331e+00  [-158 1.222e-08 1.222e-08 2.805e-06]  [40 -128 1]
 ST 1-13  3.597e+00  [-171 2.805e-09 2.805e-09 1.244e-05]  [35 -113 1]
 ST 1-14  3.865e+00  [-183 6.437e-10 6.437e-10 5.587e-05]  [30 -98 0]
 ST 1-15  4.134e+00  [-196 1.477e-10 1.477e-10 2.031e-04]  [26 -85 0]
 ST 1-16  4.400e+00  [-209 3.386e-11 3.386e-11 9.033e-04]  [21 -70 0]
 ST 1-17  4.659e+00  [-222 7.852e-12 7.852e-12 3.988e-03]  [16 -55 0]
 ST 1-18  4.929e+00  [-234 1.845e-12 1.845e-12 1.874e-02]  [11 -40 0]
 ST 1-19s 5.202e+00  [-248 3.939e-13 3.939e-13 8.439e-02]  [7 -25 0]

I increased Max_iterations and set Tolerance = 1e-12.

Expected behavior

Should converge to this tolerance in ~5 Newton iterations.

Additional context

@ktbolt, @aabrown100-git, similar to #189, was there anything changed recently in ustruct or in the Guccione material? Have you used a setup like this in a more "real world" example?

Code of Conduct

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

[ktbolt]

The original simulations I did for svFSI-Tests/06-ustruct/03-LV-Guccione-active using the parameters in svFSI.inp converged in three iterations

 ST 1-1  3.000e+00  [0 1.000e+00 1.000e+00 8.479e-05]  [27 -94 22]
 ST 1-2  7.000e+00  [-20 9.216e-02 9.216e-02 8.989e-05]  [21 -93 16]
 ST 1-3s 1.000e+01  [-81 8.256e-05 8.256e-05 9.155e-05]  [23 -93 21]

[mrp089]

That looks much better than the current convergence. However, it's still hard to say if the original simulation converged well. We would have to see how many iterations it takes the nonlinear solver to converge (i.e. not change anymore at something like 1e-11).

It's a different setup, but that's the convergence of struct/LV_Guccione_passive:

 ST 1-1  1.210e-01  [0 1.000e+00 1.000e+00 1.109e-13]  [154 -3 4]
 ST 1-2  2.420e-01  [-56 1.476e-03 1.476e-03 7.937e-11]  [126 -2 3]
 ST 1-3  3.600e-01  [-176 1.450e-09 1.450e-09 8.042e-05]  [51 -2 1]
 ST 1-4s 4.790e-01  [-209 3.494e-11 3.494e-11 3.582e-03]  [20 -1 0]

[ktbolt]

The results matched Fortran (if that is meaningful).

The Fortran tolerances are set to 1e-4.

[ktbolt]

Note that when using svSolver a residual of 1e-3 is considered to be a converged solution.

[mrp089]

What I find concerning with ustruct, cmm, and shell is not the residual tolerance in the tests but the rate of convergence. I made this convergence plot from our current test cases. You can reproduce it by setting

• Number_of_time_steps=1
• Max_iterations=35
• Tolerance=1e-16

There's a clear difference between physics that converge quadratically (struct, fluid, heats) and those that converge linearly with different rates (ustruct, cmm, shell). The former ones are linearized consistently, the latter ones not. The physics might still work in many cases; they will just be slow and not robust (this is the same problem we had with svOneDSolver before I changed it to use a finite difference linearization).

You can also nicely see when the simulations converge: The residual becomes asymptotic and just oscillates through noise. For testing, we need a converged solution to compare it across different platforms, library versions, number of processors, ... . The specific tolerance for that varies from test case to test case because it's a relative norm. That's why I went through all test cases in #175 and adapted that value for each test case.

From the discussions this week, it seems like cmm and shell converged similarily in svFSI, so we will have to live with that until someone wants to fix the linearization. For ustruct, it would be interesting to hear from @aabrown100-git what his experience was in svFSI and svFSIplus. Maybe there's a bug in svFSIplus that can be fixed.

[vvedula22]

@mrp089 ustruct/LV_Guccione_active should also show quadratic accuracy. Here is the result I have for a tolerance set to 1e-9:

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 New time step size: 5.000E-3
 ST     1-1  9.8E-1  [   0 1.00000 1.00000 9.1E-10]  [  67  -43  73]
 ST     1-2  1.9010  [ -20 9.05E-2 9.05E-2 9.3E-10]  [  59  -20  74]
 ST     1-3  2.9220  [ -81 8.05E-5 8.05E-5 8.7E-10]  [  53   -2  69]
 ST     1-4  3.6919  [-159 1.07E-8 1.07E-8 9.6E-10]  [  50 -208  65]
 ST     1-5  4.4909  [-231 2.7E-12 2.7E-12 7.4E-10]  [  51 -210  63]
 ST     1-6  5.4749  [-281 8.7E-15 8.7E-15 7.7E-10]  [  59  -19  66]
 ST     2-1  6.3929  [   0 1.00000 6.12E-1 8.8E-10]  [  68  -45  74]
 ST     2-2  7.2779  [ -18 1.17E-1 7.17E-2 8.3E-10]  [  59  -24  74]
 ST     2-3  8.1209  [ -80 9.24E-5 5.66E-5 9.1E-10]  [  53   -2  73]
 ST     2-4  8.7820  [-155 1.72E-8 1.05E-8 9.5E-10]  [  49 -208  64]
 ST     2-5  9.4340  [-219 1.1E-11 6.7E-12 9.9E-10]  [  48 -207  63]
 ST     2-6  1.03E1  [-273 2.0E-14 1.2E-14 9.9E-10]  [  59  -20  70]

Something is probably not right in the shell linearization even in the Fortran code. This is a working case I have:

 SH     1-1  8.0E-3  [   0 1.00000 1.00000 5.18E-7]  [  17 -145  13]
 SH     1-2  1.4E-2  [ -69 3.16E-4 3.16E-4 8.74E-7]  [  15 -139  17]
 SH     1-3  1.9E-2  [-150 3.12E-8 3.12E-8 1.03E-3]  [   8  -69   0]
 SH     2-1  2.5E-2  [   0 1.00000 5.35E-1 5.17E-7]  [  17 -145  14]
 SH     2-2  3.1E-2  [ -63 6.81E-4 3.65E-4 5.49E-7]  [  16 -144  17]
 SH     2-3  3.8E-2  [-129 3.49E-7 1.87E-7 1.10E-4]  [  10  -91  14]

[aabrown100-git]

Sorry for the delay. Just adding another data point to this discussion. I just ran my 0D-coupled biventricular case with ustruct and the HO-ma material model, using https://github.com/vvedula22/svFSI. Here is the convergence history

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 New time step size: 1.000E-3
 ST     1-1  4.1229  [   0 1.00000 1.00000 4.0E-12]  [ 596  -14  35]
 ST     1-2  7.6210  [ -49 3.20E-3 3.20E-3 1.28E-9]  [ 436  -10  29]
 ST     1-3  1.07E1  [-111 2.56E-6 2.56E-6 1.59E-6]  [ 278   -8  21]
 ST     1-4  1.36E1  [-172 2.43E-9 2.43E-9 1.61E-3]  [ 105   -1  10]
 ST     1-5  1.64E1  [-227 4.2E-12 4.2E-12 9.15E-1]  [   3   -1   0]
 ST     1-6  1.89E1  [-228 3.9E-12 3.9E-12 1.00000]  !   0    0   0!
 ST     1-7  2.16E1  [-227 4.3E-12 4.3E-12 8.16E-1]  [   3   -2   0]
 ST     1-8  2.43E1  [-229 3.5E-12 3.5E-12 1.00000]  !   0    0   0!
 ST     1-9  2.68E1  [-229 3.5E-12 3.5E-12 1.00000]  !   0    0   0!
 ST     1-10  2.94E1  [-225 5.2E-12 5.2E-12 6.25E-1]  [   3   -5   0]
 ST     1-11  3.18E1  [-229 3.2E-12 3.2E-12 1.00000]  !   0    0   0!
 ST     1-12  3.44E1  [-227 4.2E-12 4.2E-12 9.49E-1]  [   2   -1   0]
 ST     1-13  3.76E1  [-227 4.0E-12 4.0E-12 1.00000]  !   0    0   0!

(The not-quite-quadratic convergence may be related to the 0D-coupling)

[aabrown100-git]

I found some likely typos/discrepancies in the 1) ustruct Guccione model and 2) ustruct follower pressure load. These changes can be found in this commit: 751b966

I made these changes to match the following lines in svFSI (Fortran) code
https://github.com/SimVascular/svFSI/blob/7bdb74cdd517310d6bf32c0bb7d42072aeac8658/Code/Source/svFSI/MATMODELS.f#L809

https://github.com/SimVascular/svFSI/blob/7bdb74cdd517310d6bf32c0bb7d42072aeac8658/Code/Source/svFSI/MATMODELS.f#L812

https://github.com/SimVascular/svFSI/blob/7bdb74cdd517310d6bf32c0bb7d42072aeac8658/Code/Source/svFSI/USTRUCT.f#L1320

Unfortunately, after making these changes, the convergence for ustruct/LV_Guccione_active is still not good. There must be another bug in the tangent somewhere, but I've gone through the Guccione model in get_pk2cc_dev() several times and can't find anything.

Still, these appear to be typos so I'll submit a pull request for this.

[aabrown100-git]

I just added some material model tests to test the correct implementation of the 2nd PK stress and it's linearization, the material elasticity tensor (#225, #229). I'm wondering if it's possible and worthwhile to try a similar approach to test the linearization of various implementations in a general way. I'm thinking, for example, we could use finite difference to testing whether struct_3d_carray(..., lR, lK) produces a local tangent lK that is consistent with the local residual lR. By varying the inputs to this function, we could test the various physical models implemented in struct, which includes not just material models but also things like solid viscosity and inertia. I guess this would be testing the tangent consistency at the Gauss integration point level.

Alternatively, maybe it would be more generalizable to test higher level functions, like global_eq_assem() or other functions in iterate_solution() inner loop, in which case we could test with simple 1-element meshes (a single tet, a single hex, etc...). I guess this would basically be setting up the linear system for one Newton iteration using a 1-element mesh and testing if com_mod.Val is consistent with com_mod.R. By turning on/off various physics (boundary conditions, body forces, viscosity, contact), we could test whether their linearizations are consistent with their residual contributions.

I think by extending the mock objects currently used for material model testing, these things should be possible with a bit of work. Let me know any thought, potential issues, or other ideas!

[mrp089]

I think this is a great idea! I used an FD tangent on an element and global level in an abandoned G&R approach (maybe there's something useful in there).

You need to make an additional modification to the assembly, which currently does residual and tangent simultaneously. If you calculate the residual repeatedly as in FD, it adds to your tangent every time. I thus split up the assembly of residual and tangent (there might be a more elegant way to do this).

[aabrown100-git]

I'm making changes in this branch.

I'm not sure I need to split up the assembly of residual and tangent. So far, I've encapsulated the residual and tangent calculation into a single function, called construct_residual_and_tangent(). . My plan is to call this function within the Google testing framework to get global R and K. Then I'll create a random perturbation to the state vector (I guess Yn?) to get a perturbed R, and check if K * $\Delta$ Yn = $\Delta$ R

[ktbolt]

@aabrown100-git Like you mentioned above I think it is important to ask if this is important and worthwhile at this point in the code's evolution when lots of the code is poorly understood, not internally well-documented and should be reorganized into C++ classes.

And this Issue seems to be moving away from its initial purpose. Open another Issue so your testing ideas are not buried in an unrelated Issue.

[aabrown100-git]

@ktbolt Thanks for the thoughts. I think testing if the tangent is consistent with the residual would be somewhat useful. To provide the argument against doing this, we can more or less tell if the tangent is consistent with the residual by just observing the convergence behavior, essentially as Martin did earlier in this issue. @mrp089 Any additional thoughts on this?

[ktbolt]

@aabrown100-git If this testing helps understand the code then have at it!

@alisonmarsden @dcodoni It would be good to identify tasks and set priorities that will help evolve the code towards a real C++ implementation.