-
Notifications
You must be signed in to change notification settings - Fork 311
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
[SofaSimulationCore] Add tons of details in Euler solver #2163
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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):
- We could probably rename the "EulerSolver.h" and "EulerSolver.cpp" files to "EulerExplicitSolver.h" and "EulerExplicitSolver.cpp" to follow the name of the class.
- 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.
|
||
// Calls "solveConstraint" on every ConstraintSolver objects found in the current context tree. | ||
mop->solveConstraint(newVel,core::ConstraintParams::VEL); | ||
mop->solveConstraint(newPos,core::ConstraintParams::POS); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I noticed that too, but I don't know the answer. @ChristianDuriez can you help us on this question?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
SofaKernel/modules/SofaExplicitOdeSolver/src/SofaExplicitOdeSolver/EulerSolver.cpp
Outdated
Show resolved
Hide resolved
SofaKernel/modules/SofaExplicitOdeSolver/src/SofaExplicitOdeSolver/EulerSolver.cpp
Outdated
Show resolved
Hide resolved
So it is more similar to the other signature functions. Suggestion by jnbrunet PR#2163
I definitively agree with both your suggestions. |
One thing I changed compared to the original code is: 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 |
Unfortunately, I'm not sure I will be able to join this Wednesday and the next. Could we do something like this?
|
In your first suggestion, it would not be possible to create |
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
Second option
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... |
@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 |
[ci-build][with-all-tests] |
Hey @alxbilger , I'm very sorry for the late reply.
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! |
@jnbrunet To be clear, there are 2 optimizations in this solver: 1) optimization for diagonal masses, 2) multiop optimization. Non-optimized multiop codeI did not change this code Symplectic
Non-Symplectic
Optimized multiop codeops_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);
If I translate this code to formulas: Symplectic
Non-symplectic
This is exactly the same formulas than for the non-optimized code! It would translate to: Symplectic
Non-symplectic
ConclusionI 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? |
Hey @alxbilger , yes go ahead, you can absolutely merge this. The issue I think I see here isn't related to this PR. |
solve
function easier to read)Next steps:
d_optimizedForDiagonalMatrix
.SOFA_NO_VMULTIOP
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