Non linear implicit solver

[jjcasmar]

I have looked at the EulerImplicitSolver and it only does a Newton step integration. I am trying to implement a version of it doing several Newton steps for integrating more complex energies and systems but I am stuck on the design of it.

Normally, when I do this kind of solvers outside of SOFA, I implement a solve very similar to what EulerImplicitSolver has, do a lineSearch using the direction provided by the solve and take the step when the error is smaller than the error I had.

Then I repeat this linearizing in the new state until the system error is smaller than a cerain value.

In SOFA is not possible to implement this directly as I need to update the state of the system, which is not done directly from the solver from what I have seen.

Could someone guide me on how to implement a non linear newton solver?

[jnbrunet]

Hey @jjcasmar,

Have a look at the StaticSolver which does a series of Newton-Raphson iterations. It should give you hints on how to achieve such state updates within the same time step.

Let us know if you hit specific difficulties and don’t hesitate to show your work to the community, this is a component that could be very useful to others.

Jean-Nicolas

[jjcasmar]

Great, that is what I was looking for.

So if I understand correctly,
1) Compute forces and current error to bootstrap the solver
2) Compute a dx
3) Compute x as x0 + dx
4) Propagate this changes through the graph using MechanicalPropagateOnlyPositionAndVelocityVisitor
5) Recompute forces, now in the new state
6) Repea

Somewhere between 4 and 5, the mappings have to be updated. Who is in charge of this? MechanicalPropagateOnlyPositionAndVelocityVisitor or mop.computeForces()

I will hack something fast later and push it to github so you can have a look.

[jnbrunet]

Hey @jjcasmar,

You understood correctly.

mop.computeForces() is only responsible to propagate forces from the mapped mechanical object to their parent.

MechanicalPropagateOnlyPositionAndVelocityVisitor is responsible to propagate the position and velocity vectors from parent (top level, or system) vectors to mapped vectors.

You can have a look at this StaticSolver from the Caribou plugin. It is a little bit more advanced compared to the SOFA one and its documentation is also richer.

[jjcasmar]

I think I still have one doubt.

From Caribou

        // Updating the geometry
        x.peq(dx); // x := x + dx

        // Solving constraints
        // Calls "solveConstraint" method of every ConstraintSolver objects found in the current context tree.
        // todo(jnbrunet): Shouldn't this be done AFTER the position propagation of the mapped nodes?
        mop.solveConstraint(x, sofa::core::ConstraintParams::POS);

        // Propagate positions to mapped mechanical objects, for example, identity mappings, barycentric mappings, ...
        // This will call the methods apply and applyJ on every mechanical mappings.
        sofa::simulation::MechanicalPropagateOnlyPositionAndVelocityVisitor(&mparams).execute(this->getContext());

You set x as x+dx. Is that operation already triggering to store the new values in the mechanical objects or when is that done? I understand that the mappings are update using the visitor, but the DoF values should be set somewhere, no?

[jnbrunet]

Yeah this part isn’t very straightforward to understand.

Here, x.peq(dx); will store the new positions inside the mechanical object’s vector pointed by x. Now, which one is x? If you follow its definition, it is set to the solve() method’s parameter sofa::core::MultiVecCoordId xResult which is defaulted to VecCoordId::position(), hence the MechanicalObject.position datafield. Now, the x.peq(dx); operation starts a visitor that stops at every mapped mechanical objects, i.e. any mechanical object where a mechanical mapping is found.

Next, MechanicalPropagateOnlyPositionAndVelocityVisitor will go through all mapped mechanical objects and compute their new position (MechanicalObject.position datafield) using mapped_x = mapping.J*parent.x.

Let’s say you have this scene:

<Node name="root">
  <StaticSolver />
  <LinearSolver />
  <Node name="object1">
    <MechanicalObject name="mo" />
    <Forcefield name="volume_ff" />
    <Node name="surface">
      <MechanicalObject name="mo" />
      <Forcefield name="surface_ff" />
      <Mapping name='mapping'/>
    </Node>
  </Node>

  <Node name="object2">
    <MechanicalObject name="mo" />
    <Forcefield name="volume_ff" />
  </Node>
</Node>

It goes like this:

1. x.peq(dx);
root.object1.mo.position += root.object1.mo.dx
root.object2.mo.position += root.object2.mo.dx

2. MechanicalPropagateOnlyPositionAndVelocityVisitor
root.object1.surface.position = root.object1.surface.mapping.J * root.object1.mo.position

Hope that cleared things a little bit.
J-N

[jjcasmar]

I think this solves all the issue.

I think you did a small typo in
root.object2.mo.position += root.object1.mo.dx
right?

[jnbrunet]

Yes sorry, I’ve corrected it.

[Author]

Posts