[SofaSimulationCore] Add tons of details in Euler solver

[alxbilger]

• Lots of comments
• Split into multiple functions (solve function easier to read)

Next steps:

• Create a visitor to detect automatically if the mass matrix is trivially invertible (diagonal matrix), in order to replace d_optimizedForDiagonalMatrix.
• Benchmark with/without SOFA_NO_VMULTIOP
• More tests, in particular for the standard Euler method (not symplectic)
• Update the doc

---

By submitting this pull request, I acknowledge that
I have read, understand, and agree SOFA Developer Certificate of Origin (DCO).

---

Reviewers will merge this pull-request only if

• it builds with SUCCESS for all platforms on the CI.
• it does not generate new warnings.
• it does not generate new unit test failures.
• it does not generate new scene test failures.
• it does not break API compatibility.
• it is more than 1 week old (or has fast-merge label).

[jnbrunet]

Hey @alxbilger ,

Awesome PR, thanks a lot. Reading through this class will now be much clearer and approachable.

Small general notes (not necessarily for this PR):

1. We could probably rename the "EulerSolver.h" and "EulerSolver.cpp" files to "EulerExplicitSolver.h" and "EulerExplicitSolver.cpp" to follow the name of the class.
2. Maybe we could factor the sympletic option as a new time integrator class? I think it would be much clearer to have, for example,

<EulerExplicitSolver />
<EulerImplicitSolver />
<EulerSemiImplicitSolver />

and might speak to more people.

[jnbrunet]

I know you didn't introduce this and it was like that before, but, don't you think this "single-operation optimization" introduces a different behavior with respect to constraints? i.e.:

Without optimization:

v = v + a*t
solve_constraints(v)
x = x + v*t
solve_constraints(x)

With optimization:

v = v + a*t
x = x + v*t
solve_constraints(v)
solve_constraints(x)

[alxbilger]

I noticed that too, but I don't know the answer. @ChristianDuriez can you help us on this question?

[alxbilger]

Short answer: it depends on how the constraint is implemented.

In the non-optimized version, constraints are applied on the velocity, and then the position is updated. Therefore, the new position takes into account the changes due to the constraints on the velocity. Finally, constraints on the position are applied.
In the optimized version, the position is updated first, and then the constraints on the velocity are applied. Therefore, the new position does not take into account the constraints on the velocity. However, it may have the same effect, if the constraint applied on the position has the same effect than the constraints applied on the velocity.
To be completely rigorous, I think the non-optimized version is the good one. However, in practice the optimized version works similarly to the non-optimized version, because constraints on velocity and position are consistent.
Let's see if it worths keeping the "optimized" version after a benchmark.

[alxbilger]

1. We could probably rename the "**EulerSolver.h**" and "**EulerSolver.cpp**" files to "**EulerExplicitSolver.h**" and "**EulerExplicitSolver.cpp**" to follow the name of the class.

2. Maybe we could factor the sympletic option as a new time integrator class? I think it would be much clearer to have, for example,

<EulerExplicitSolver />
<EulerImplicitSolver />
<EulerSemiImplicitSolver />

and might speak to more people.

I definitively agree with both your suggestions.
I also thought about your second suggestion, and it was my plan to talk about it during the next dev meeting. The problem I see is that EulerExplicitSolver is by default symplectic. So, applying your suggestion would change the ODE solver when the user write <EulerExplicitSolver />. Let's discuss it on Wednesday.

[alxbilger]

One thing I changed compared to the original code is:
from

ops[op_pos].second.push_back(std::make_pair(newVel.id(),dt));

to

ops_pos.second.emplace_back(d_symplectic.getValue() ? newVel.id() : vel.id(), dt);

Instead of providing newVel as the first argument, I check d_symplectic, and if it's false, I provide vel. I do that because I think the non-symplectic variant works on the velocity of the previous time step. Also, I think newVel is not yet computed for the non-symplectic variant. I would like your opinion @jnbrunet

[jnbrunet]

1. We could probably rename the "**EulerSolver.h**" and "**EulerSolver.cpp**" files to "**EulerExplicitSolver.h**" and "**EulerExplicitSolver.cpp**" to follow the name of the class.

2. Maybe we could factor the sympletic option as a new time integrator class? I think it would be much clearer to have, for example,

<EulerExplicitSolver />
<EulerImplicitSolver />
<EulerSemiImplicitSolver />

and might speak to more people.

I definitively agree with both your suggestions.
I also thought about your second suggestion, and it was my plan to talk about it during the next dev meeting. The problem I see is that EulerExplicitSolver is by default symplectic. So, applying your suggestion would change the ODE solver when the user write <EulerExplicitSolver />. Let's discuss it on Wednesday.

Unfortunately, I'm not sure I will be able to join this Wednesday and the next.

Could we do something like this?

1. If EulerExplicitSolver is used with default parameters or with explicit sympletic=True, then warn the user that he should use EulerSemiImplicitSolver instead, and create this new component automatically for him (this is the part where I'm not sure it's feasible..)
2. Else, the user manually specified sympletic=False, warn him that this option no longer exists and that he can remove it from its scene file.

[alxbilger]

Could we do something like this?

1. If `EulerExplicitSolver` is used with default parameters or with explicit `sympletic=True`, then warn the user that he should use `EulerSemiImplicitSolver` instead, and create this new component automatically for him (this is the part where I'm not sure it's feasible..)

2. Else, the user manually specified `sympletic=False`, warn him that this option no longer exists and that he can remove it from its scene file.

In your first suggestion, it would not be possible to create EulerExplicitSolver.

[jnbrunet]

One thing I changed compared to the original code is:
from

ops[op_pos].second.push_back(std::make_pair(newVel.id(),dt));

to

ops_pos.second.emplace_back(d_symplectic.getValue() ? newVel.id() : vel.id(), dt);

Instead of providing newVel as the first argument, I check d_symplectic, and if it's false, I provide vel. I do that because I think the non-symplectic variant works on the velocity of the previous time step. Also, I think newVel is not yet computed for the non-symplectic variant. I would like your opinion @jnbrunet

Well, your question took me some time haha. I'm not sure about this... This is what I have at the moment when I try to derive them manually:

First option

Begin time step
Solve a in : (M + hC) a = P -  C v - R(x)
x = x + h v
v = v + h a
solve_constraints(x)
solve_constraints(v)
End time step

Second option

Begin time step
x = x + h v + h^2 a
v = v + h a
solve_constraints(x)
solve_constraints(v)
Solve a in : M a = P -  C v - R(x)
End time step

I've put my maths from the middle of the night here:

https://www.overleaf.com/read/mdbcmmnqxzkk

But I might have made some mistakes... I don't get how SOFA doesn't have the same thing...

Edit: To clarify, I didn't really answer your question because I just can't reproduce neither of those equations from the ones found in the theory...

[alxbilger]

@jnbrunet If I understand correctly, your problem with this solver is not the approximation equations on which I asked your opinion, but with the way a_{n+1} is computed. From your document, a_{n+1} should be computed differently depending on if it symplectic or not. In EulerExplicitSolver, we just compute a in Ma=f for both symplectic and non-symplectic, right?

[alxbilger]

[ci-build][with-all-tests]

[jnbrunet]

Hey @alxbilger ,

I'm very sorry for the late reply.

If I understand correctly, your problem with this solver is not the approximation equations on which I asked your opinion

Well kinda. I am unable to answer your initial question, since when I try to derive the equation to validate your changes, I'm unable to get to the same formulation as you have. But, in addition, I'm also unable to derive the old one (before your change)....

Anyway, if you can validate that the multiop "optimized" version gives the same result as the "unoptimized" one for both explicit and semi-implicit, and for both before and after this PR, you have my +1 to merge!

[alxbilger]

@jnbrunet To be clear, there are 2 optimizations in this solver: 1) optimization for diagonal masses, 2) multiop optimization.
I did not touch anything in the non-optimized multiop code. I changed the velocity input in the position update in the optimized multiop code.

Non-optimized multiop code

I did not change this code

Symplectic

newVel = vec + acc * dt
newPos = pos + newVel * dt

Non-Symplectic

newPos = pos + vel * dt
newVel = vel + acc * dt

Optimized multiop code

ops_vel.second.emplace_back(vel.id(), 1.0);
ops_vel.second.emplace_back(acc.id(), dt);
ops_pos.second.emplace_back(pos.id(), 1.0);
ops_pos.second.emplace_back(d_symplectic.getValue() ? newVel.id() : vel.id(), dt);

ops_vel is executed first if symplectic. ops_pos is executed first if non-symplectic.

If I translate this code to formulas:

Symplectic

newVel = vel + acc * dt
newPos = pos + newVel * dt

Non-symplectic

newPos = pos + vel * dt
newVel = vel + acc * dt

This is exactly the same formulas than for the non-optimized code!
My change was on the line:
ops_pos.second.emplace_back(d_symplectic.getValue() ? newVel.id() : vel.id(), dt);
Before it was:
ops_pos.second.emplace_back(vel.id(), dt);

It would translate to:

Symplectic

newVel = vel + acc * dt
newPos = pos + vel * dt  <-- difference here with the non-optimized code

Non-symplectic

newPos = pos + vel * dt
newVel = vel + acc * dt

Conclusion

I believe that I fixed the multiop optimized symplectic version. Both optimized and non-optimized codes are now consistent. Note that this makes a difference only for the FreeMotion animation loops.

@jnbrunet Regarding your math problem, I don't understand it. I don't see where you see a difference with SOFA.

If you agree, I think we can merge this PR because it does not change the behavior (except for the fix). If there are some maths problems, it can probably be tackled later. What do you think?

[jnbrunet]

Hey @alxbilger , yes go ahead, you can absolutely merge this. The issue I think I see here isn't related to this PR.