Caribou manufactured solution isn't converging

[jnbrunet]

The manufactured solution for the StVk material isn't converging anymore.

Time step                : 5
Context                  : /
Max iterations           : 10
Residual tolerance (abs) : 1e-15
Residual tolerance (rel) : 1e-10
Correction tolerance     : 1e-05
Linear solver            : /LDLTSolver

Newton iteration #1      |R| = 3.070686e+82  |R|/|R0| = 1.000000e+00  |du| / |U| = 1.000000e+00  Time = 6 ms
Newton iteration #2      |R| = 9.098330e+81  |R|/|R0| = 2.962963e-01  |du| / |U| = 4.000000e-01  Time = 3 ms
Newton iteration #3      |R| = 2.695802e+81  |R|/|R0| = 8.779150e-02  |du| / |U| = 2.105263e-01  Time = 3 ms
Newton iteration #4      |R| = 7.987560e+80  |R|/|R0| = 2.601229e-02  |du| / |U| = 1.230769e-01  Time = 3 ms
Newton iteration #5      |R| = 2.366684e+80  |R|/|R0| = 7.707347e-03  |du| / |U| = 7.582938e-02  Time = 3 ms
Newton iteration #6      |R| = 7.012399e+79  |R|/|R0| = 2.283658e-03  |du| / |U| = 4.812030e-02  Time = 3 ms
Newton iteration #7      |R| = 2.077748e+79  |R|/|R0| = 6.766395e-04  |du| / |U| = 3.108305e-02  Time = 3 ms
Newton iteration #8      |R| = 6.156290e+78  |R|/|R0| = 2.004858e-04  |du| / |U| = 2.030135e-02  Time = 3 ms
Newton iteration #9      |R| = 1.824086e+78  |R|/|R0| = 5.940319e-05  |du| / |U| = 1.335350e-02  Time = 3 ms
Newton iteration #10     |R| = 5.404699e+77  |R|/|R0| = 1.760095e-05  |du| / |U| = 8.823783e-03  Time = 3 ms
[DIVERGED] The number of Newton iterations reached the maximum of 10 iterations.

relative L2 error at 100% load: 125780184395449.69

[jnbrunet]

It looks like the p1 mesh is wrong. It is converging the p2 mesh. Using the following changes:

diff --git a/Validation/manufactured_solution.py b/Validation/manufactured_solution.py
index 8734ec1..c7cb6a2 100644
--- a/Validation/manufactured_solution.py
+++ b/Validation/manufactured_solution.py
@@ -22,7 +22,7 @@ import Sofa, SofaRuntime, SofaCaribou
 import meshio, numpy as np
 from manufactured_solution import assemble, integrate, compute_solution, ConstantForceField

-mesh = meshio.read('meshes/cylinder_p1.vtu')
+mesh = meshio.read('meshes/cylinder_p2.vtu')

 mu = 1.0
 l  = 1.25
@@ -41,9 +41,9 @@ def create_scene(root):
     root.addObject('StaticODESolver', newton_iterations=10, residual_tolerance_threshold=1e-10, printLog=True)
     root.addObject('LDLTSolver', backend='Pardiso')
     root.addObject('MechanicalObject', name='mo', position=mesh.points.tolist())
-    root.addObject('CaribouTopology', name='volume', template='Tetrahedron', indices=mesh.cells[1].data.tolist())
-    root.addObject('CaribouTopology', name='dirichlet_boundary', template='Triangle', indices=mesh.cells[0].data[np.ma.masked_equal(mesh.cell_data['gmsh:physical'][0], 1).mask].tolist())
-    root.addObject('CaribouTopology', name='neumann_boundary',   template='Triangle', indices=mesh.cells[0].data[np.ma.masked_inside(mesh.cell_data['gmsh:physical'][0], 2, 3).mask].tolist())
+    root.addObject('CaribouTopology', name='volume', template='Tetrahedron10', indices=mesh.cells[1].data.tolist())
+    root.addObject('CaribouTopology', name='dirichlet_boundary', template='Triangle6', indices=mesh.cells[0].data[np.ma.masked_equal(mesh.cell_data['gmsh:physical'][0], 1).mask].tolist())
+    root.addObject('CaribouTopology', name='neumann_boundary',   template='Triangle6', indices=mesh.cells[0].data[np.ma.masked_inside(mesh.cell_data['gmsh:physical'][0], 2, 3).mask].tolist())

     root.addObject('SaintVenantKirchhoffMaterial', young_modulus=young_modulus, poisson_ratio=poisson_ratio)
     root.addObject('HyperelasticForcefield', topology='@volume', printLog=True)

leads to convergence:

relative L2 error at 0% load: 0.9748987206941837
relative L2 error at 1% load: 0.895999603032051
relative L2 error at 10% load: 0.6415177176178076
relative L2 error at 15% load: 0.5652613135083865
relative L2 error at 50% load: 0.2536373549078224
relative L2 error at 100% load: 8.300875725328914e-06

[Ziemnono]

@jnbrunet, I quickly checked and if you open the 2 meshes (cylinder_p1.vtu and cylinder_p2.vtu) with Paraview, you will clearly see a difference in size in the 2 meshes. In wireframe, the p1 cylinder and the small one is the p2.

I can reproduce the same results for the p2 but when it comes to p1, I just adjust the length to the corresponding one (around 80 if I am not mistaken). Leading to the following output:

[INFO]    [StaticODESolver] ======= Starting static ODE solver =======
Time step                : 5
Context                  : /
Max iterations           : 10
Residual tolerance (abs) : 1e-15
Residual tolerance (rel) : 1e-10
Correction tolerance     : 1e-05
Linear solver            : /LDLTSolver

Newton iteration #1      |R| = 4.332215e+03  |R|/|R0| = 1.000000e+00  |du| / |U| = 1.000000e+00  Time = 7 ms
Newton iteration #2      |R| = 2.135447e+02  |R|/|R0| = 4.929226e-02  |du| / |U| = 2.521153e-01  Time = 3 ms
Newton iteration #3      |R| = 6.182438e-01  |R|/|R0| = 1.427085e-04  |du| / |U| = 1.338070e-02  Time = 3 ms
Newton iteration #4      |R| = 5.294556e-06  |R|/|R0| = 1.222136e-09  |du| / |U| = 3.618296e-05  Time = 2 ms
Newton iteration #5      |R| = 1.918934e-10  |R|/|R0| = 4.429452e-14  |du| / |U| = 2.697594e-10  Time = 2 ms
[CONVERGED] The correction's ratio |du|/|U| = 2.69759e-10 is smaller than the threshold of 1e-05.

relative L2 error at 100% load: 0.007460799780052525

What do you think?

[jnbrunet]

Nice find @Ziemnono !

The radius and length of the cylinder mesh must indeed match to one provided to the manufactured solution.

I guess the correct fix here is to update the mesh so that it matches to parameter of the manufactured solution. It would also be nice to add it to the CI unittests, I will have a look at it.

[Ziemnono]

Hi @jnbrunet, With @Sidaty1, we created 2 new meshes (p1_from_gmsh.msh and p2_from_gmsh.msh) that you can find here. We accordingly changed the length and radius in the manufactured_solution.py to end up with the following results for p1 and p2 respectively:

relative L2 error at 100% load: 0.007885064961087717
relative L2 error at 100% load: 9.582959168910188e-06

We deliberately chose coarse meshes to avoid excessive computational time but by refining, you definitely reach a much better accuracy. I will let you close the issue if this solution looks good to you. We can also provide you with finer meshes if you want to go below e-10 relative L2 error.