[jjcasmar]

I have kind of solve this. I have reimplemented the CG solver using the one in Caribou. Based on that, I have kept only the implementation that doesn’t build the matrix, but added some calls to assemble the vector that I had to project after preconditioning according to the paper linked above.

I copy paste the implementation here, and when I have time, Ill try to push upstream.

void ConjugateGradientSolver::solve(sofa::core::behavior::MultiVecDeriv &b, sofa::core::behavior::MultiVecDeriv &x)
{
    sofa::simulation::common::VectorOperations vop(&m_mparams, this->getContext());
    sofa::simulation::common::MechanicalOperations mop(&m_mparams, this->getContext());

    // Create temporary vectors needed for the method
    sofa::core::behavior::MultiVecDeriv r(&vop);
    sofa::core::behavior::MultiVecDeriv p(&vop);
    sofa::core::behavior::MultiVecDeriv s(&vop);
    sofa::core::behavior::MultiVecDeriv h(&vop);

    // Get the method parameters
    const auto &maximum_number_of_iterations = m_iterations;
    const auto &residual_tolerance_threshold = m_toleranceThreshold;

    m_squaredResiduals.clear();
    m_squaredResiduals.reserve(m_iterations);

    // Get the matrices coefficient m, b and k : A = (mM + bB + kK)
    const auto m_coef = m_mparams.mFactor();
    const auto b_coef = m_mparams.bFactor();
    const auto k_coef = m_mparams.kFactor();

    // Declare the method variables
    VNCS::Real b_delta;
    VNCS::Real delta;            // RHS and residual squared norms
    VNCS::Real rho0, rho1 = 0.;  // Stores r*r as it is used two times per iterations
    VNCS::Real alpha, beta;      // Alpha and Beta coefficients
    VNCS::Real threshold;        // Residual threshold
    int currentIteration = 0;    // Current iteration number
    bool converged = false;

    sofa::component::linearsolver::DefaultMultiMatrixAccessor accessor;

    mop.getMatrixDimension(nullptr, nullptr, &accessor);
    const auto n = static_cast<Eigen::Index>(accessor.getGlobalDimension());

    if (m_preconditioner) {
        Eigen::VectorXd bEigen;
        bEigen.resize(n);
        sofa::component::linearsolver::EigenVectorWrapper<VNCS::Real> bEigenWrapped(bEigen);
        mop.multiVector2BaseVector(b.id(), &bEigenWrapped, &accessor);

        b_delta = bEigen.transpose() * m_preconditioner.value() * bEigen;

    } else {
        b_delta = b.dot(b);
    }

    threshold = residual_tolerance_threshold * residual_tolerance_threshold * b_delta;

    // INITIAL RESIDUAL
    // Do the A*x(0) with visitors since we did not construct the matrix A
    mop.propagateDxAndResetDf(x, s);                 // Set s = 0 and calls applyJ(x) on every mechanical mappings
    mop.addMBKdx(s, m_coef, b_coef, k_coef, false);  // s = (m M + b B + k K) x

    // Do the projection of the result in the constrained space
    mop.projectResponse(s);  // BaseProjectiveConstraintSet::projectResponse(q)

    // Finally, compute the initial residual r = b - A*x
    r.eq(b, s, -1.0);  // r = b - A*x

    Eigen::VectorXd rEigen;
    sofa::component::linearsolver::EigenVectorWrapper<VNCS::Real> rEigenWrapped(rEigen);
    if (m_preconditioner) {
        rEigen.resize(n);

        mop.multiVector2BaseVector(r.id(), &rEigenWrapped, &accessor);

        rEigen = m_preconditioner.value() * rEigen;  // It doesnt alias because m_preconditioner is diagonal

        mop.baseVector2MultiVector(&rEigenWrapped, p.id(), &accessor);  // Save it directly in p, (p = P * r)
        mop.projectResponse(p);                                         // p = S * P * r
    } else {
        p.eq(r);
    }

    delta = r.dot(p);

    // ITERATIONS
    while (delta > threshold && currentIteration < maximum_number_of_iterations) {
        // 1. Computes q(k+1) = A*p(k)
        mop.propagateDxAndResetDf(p, s);                 // Set s = 0 and calls applyJ(p) on every mechanical mappings
        mop.addMBKdx(s, m_coef, b_coef, k_coef, false);  // s = (m M + b B + k K) p 

        // We need to project the residual in the constrained space since the constraints haven't been added to the
        // matrix (the matrix is never constructed) and addMBKdx of the forcefields do not take constraints into
        // account.
        mop.projectResponse(s);  // BaseProjectiveConstraintSet::projectResponse(q)

        // 2. Computes x(k+1) and r(k+1)
        alpha = delta / p.dot(s);
        x.peq(p, alpha);   // x = x + alpha*p
        r.peq(s, -alpha);  // r = r - alpha*q

        if (m_preconditioner) {
            mop.multiVector2BaseVector(r.id(), &rEigenWrapped, &accessor);

            rEigen = m_preconditioner.value() * rEigen;

            mop.baseVector2MultiVector(&rEigenWrapped, h.id(), &accessor);  // Save it directly in h, (h = P * r)
        } else {
            h.eq(r);
        }

        // 3. Computes the new residual norm
        VNCS::Real deltaOld = delta;
        delta = r.dot(h);

        // 6. Compute the next search direction
        p.eq(h, p, delta / deltaOld);  // p = r + beta*p
        mop.projectResponse(p);

        ++currentIteration;
    }

    auto log = spdlog::get("VNCS");
    if (log) {
        log->info("Conjugate gradient finished");
        log->info("Iterations = {}/{}", currentIteration, maximum_number_of_iterations);
        log->info("Error = {:.4e}/{:.4e}", delta, threshold);
    }
}

[jjcasmar]

I have been looking at the code from Caribou, as well as the code from ShewchukPCGLinearSolver and I have posted a possible bug on the github tracker. If I understand the code correctly, there are some differences with the algorithm described in https://www.cs.ubc.ca/~ascher/papers/ab.pdf

When I have that clear, I will have a better look on how to implement a custom LinearSolver for my preconditioner.

One thing that is a bit unclear right now is that the EulerImplicitSolver requires a LinearSolver and it uses the first one it finds in the tree. If that one is a ShewchukPCGLinearSolver, it is possible there is another LinearSolver on the same node, no? Therefore, it is important to keep the one use for the preconditioner the second one in the node so the EulerImplicitSolver doesn’t find it, no?

[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?

[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?

[jjcasmar]

So I still need to update the mapping to propagate the positions, right?

[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.

[jjcasmar]

the mop.computeForce takes care of updating the mappings in both directions?

[jjcasmar]

Sorry, cant share the scene description right now. Ill do it ASAP.

One thing I normally dont have clear if is I have to add the kFactor or dt on the computations and also the sign of the result.

For example, if used in a EulerImplicit solver, the force is multiplied by the dt and the hessian is multiplied by dt+beta value (or something like that, I dont have the details right infront of me). It would be nice to add in the documentation if we need to multiply our force and hessian by these quantities.

Also, the force is the -dE/du, but I always have to double check if the solver is assuming Im returning the force or the derivative of the energy. It would be nice to specify all this in the documentation clearly.

[jjcasmar]

@Froy can you provide a bit more details or documentation on how to do this? I am trying to mimic what I see in other types, but I fail to actually make it work. I am having linking errors when I try to load my plugin.

[ERROR]   [PluginManager] Plugin loading failed (/home/jjcasmar/projects/VNCSSofaPlugin/build/Debug/lib/libVNCSSofaPlugin.so): /home/jjcasmar/projects/VNCSSofaPlugin/build/Debug/lib/libVNCSSofaPlugin.so: undefined symbol: _ZTv0_n88_N4sofa4core5StateINS_11defaulttype14StdVectorTypesINS2_3VecILi9EdEES5_dEEE9baseWriteENS0_6TVecIdILNS0_7VecTypeE0ELNS0_9VecAccessE1EEE
[ERROR]   [RequiredPlugin(/home/jjcasmar/projects/VNCSSofaPlugin/build/Debug/lib/libVNCSSofaPlugin.so)] Plugin loading failed (/home/jjcasmar/projects/VNCSSofaPlugin/build/Debug/lib/libVNCSSofaPlugin.so): /home/jjcasmar/projects/VNCSSofaPlugin/build/Debug/lib/libVNCSSofaPlugin.so: undefined symbol: _ZTv0_n88_N4sofa4core5StateINS_11defaulttype14StdVectorTypesINS2_3VecILi9EdEES5_dEEE9baseWriteENS0_6TVecIdILNS0_7VecTypeE0ELNS0_9VecAccessE1EEE

I have tested two approaches:
1) Creating my own type which exposes Eigen::Matrix3d as Coord and Deriv types. Since this looks complicated, I decided to go to case 2
2) Using StdVectorTypes with a Vec9 for Coord and Deriv. This is the case I am trying to make it work right now.

Some documentation on this will be really useful. I think that being able to use our own types in MechanicalObject gives a lot of freedom. What I am trying here is not possible wihout a custom templated MechanicalObject…

[jjcasmar]

Thanks for the answer. I guess I will wait until that is ready and code that part in C++ for now.

[jjcasmar]

Well, that depends on your preconditioner. For my particular use case, the preconditioned can be represented as a matrix, but its easier to represent it as some functions, kind of matrix-free preconditioner. I think saying that a its based on a external linear direct solver may arise to some misconceptions later. Better to just say that is has to be a linear solver, and maybe then cite a few and add “or any user implementation” or something like that.

Ill close the question and mark it as solved, as the original question is solved.

Thanks!

[jjcasmar]

Thanks Hugo. I have read the document you proposed and I think that this: The ShewchukPCGLinearSolver allows to choose the preconditioner of our choice based on an external direct linear solver: LULinearSolver, SparseLDLSolver, etc Is not exactly correct.

From reading the source code, I think any linear solver can be used, not only direct linear solver. This also includes user defined implementations.

Maybe at some point SOFA will be interested in deprecating CGLinearSolver in favour of Shewchuk?

[jjcasmar]

So ShewchukPCGSolver is a preconditioned conjugate gradient implementation. Why is not the default CG implementation? Seems like it is a bit hidden…

If I understand correctly, I just need to provide an implementation of my preconditioner, deriving it from sofa::core::behavior::LinearSolver and link it to the preconditioners property, right?

[jjcasmar]

Okay, was my problem. Didnt realize the the example scene I was looking to had the fixed collider outside the scope of the ODESolver and the CGSolver. I have put my shere outside the scope of my ODESolver and now is fixed.

[jjcasmar]

I have compiled in debug mode to check what is happening and I confirm that masses are mapped through the mapping, at least when using a matrix-free solver like CG

[jjcasmar]

Okay, I see the point. If the DataEngine initializing the data is after the MO, the MO init function will be called before the init function of my DataEngine so the MO will see nothing when executhing the link, no?

[jjcasmar]

Not sure I can share my scene. It uses some plugin code I can’t share atm.
I am trying to initialize my degrees of freedom using a DataEngine. Still, its a bit unclear what is the order of initialization of the different components in a scene. Any hint about this?

[jjcasmar]

Apparently the plugin was compiled without debug symbols, not sure why. I have wiped the plugin build directory and rerun cmake and now the compiler is adding the correct flag for building the debug symbols.

With that, in case anyone need it in the future: I am using linux and qtcreator as IDE. I just need to ‘Start and Debug an external application’ and run ‘runSofa’ with the scene you want. Breakpoints work fine for the plugin. SOFA itself can be compiled without debug symbols and with -O2 or -O3.

[jjcasmar]

Wow, awesome work! Ill check ASAP!

[jjcasmar]

Hi @Hugo, you can find all the code in the gitlab (its easier to share it that way). Anyway, @jnbrunet solution works fine.
Shouldnt the .inl be already included in sofa/core/MultiMapping.h?

[jjcasmar]

Would it make sense to provide some default integration schemes in which you can specify the number of Newton steps done per time step?

[jjcasmar]

Please, finish the “Main principles” documentation!

I think that it will be nice to have some kind of tool to look for a component and get feedback about the available templates and the data names for the .scn file.

I know that with the Modeler you can do this, but I think that a lighter application to just search for a component and getting the information about that component will be nice. There’s no need for all the extra tools available in the Modeler (launch sofa, tutorials…)

I dont know if SOFA already has multithread support in the core, but a task scheduler mechanism would be great. Just schedule tasks for updating maps, computing forces, whatever, with their dependencies during graph traverse and run afterwards in parallel

Thanks for building sucha a great tool. Prototyping new algorithms is easier with SOFA 😉

[jjcasmar]

Hi Hugo,

in fact it have some relationship as I’m simply trying to understand how constraints work.

The main issue I’m facing in this example is that the LMConstraints are not being applied, but the solver is added in the scene :-s

Thanks!

[jjcasmar]

Thanks Hugo!

I will have a look at StateMask and try to compile the plugin.

[jjcasmar]

Finally Ive been able to compile v16.12. The problem was that FileRepository.cpp wasn’t added to the project so some functions weren’t found. This problem was due to some missing boost libraries (in particular, FileSystem).

if(Boost_FILESYSTEM_FOUND AND Boost_LOCALE_FOUND)
    list(APPEND HEADER_FILES system/FileRepository.h)
    list(APPEND SOURCE_FILES system/FileRepository.cpp)
endif()

Maybe there should be some kind of check to verify that the needed boost libs are installed?

Thanks

[Author]

Posts