[ODESolver] Introduction of generic Newton-Raphson algorithm

Sofa/Component/ODESolver/Backward/CMakeLists.txt

@@ -10,6 +10,7 @@ set(HEADER_FILES
     ${SOFACOMPONENTODESOLVERBACKWARD_SOURCE_DIR}/StaticSolver.h
     ${SOFACOMPONENTODESOLVERBACKWARD_SOURCE_DIR}/VariationalSymplecticSolver.h
     ${SOFACOMPONENTODESOLVERBACKWARD_SOURCE_DIR}/NewmarkImplicitSolver.h
+    ${SOFACOMPONENTODESOLVERBACKWARD_SOURCE_DIR}/NewtonRaphsonSolver.h
 )
 
 set(SOURCE_FILES
@@ -18,6 +19,7 @@ set(SOURCE_FILES
     ${SOFACOMPONENTODESOLVERBACKWARD_SOURCE_DIR}/StaticSolver.cpp
     ${SOFACOMPONENTODESOLVERBACKWARD_SOURCE_DIR}/VariationalSymplecticSolver.cpp
     ${SOFACOMPONENTODESOLVERBACKWARD_SOURCE_DIR}/NewmarkImplicitSolver.cpp
+    ${SOFACOMPONENTODESOLVERBACKWARD_SOURCE_DIR}/NewtonRaphsonSolver.cpp
 )
 
 sofa_find_package(Sofa.Simulation.Core REQUIRED)

Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/NewtonRaphsonSolver.h

@@ -0,0 +1,87 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2006 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* This program is free software; you can redistribute it and/or modify it     *
+* under the terms of the GNU Lesser General Public License as published by    *
+* the Free Software Foundation; either version 2.1 of the License, or (at     *
+* your option) any later version.                                             *
+*                                                                             *
+* This program is distributed in the hope that it will be useful, but WITHOUT *
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
+* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
+* for more details.                                                           *
+*                                                                             *
+* You should have received a copy of the GNU Lesser General Public License    *
+* along with this program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: [email protected]                             *
+******************************************************************************/
+#pragma once
+#include <sofa/component/odesolver/backward/config.h>
+
+#include <sofa/core/behavior/LinearSolverAccessor.h>
+#include <sofa/core/behavior/OdeSolver.h>
+#include <sofa/core/behavior/BaseIntegrationMethod.h>
+#include <sofa/simulation/MechanicalOperations.h>
+
+
+namespace sofa::component::odesolver::backward
+{
+
+class SOFA_COMPONENT_ODESOLVER_BACKWARD_API NewtonRaphsonSolver
+    : public sofa::core::behavior::OdeSolver
+    , public sofa::core::behavior::LinearSolverAccessor
+{
+public:
+    SOFA_CLASS2(NewtonRaphsonSolver, sofa::core::behavior::OdeSolver, sofa::core::behavior::LinearSolverAccessor);
+
+    void init() override;
+    void reset() override;
+    void solve(const core::ExecParams* params, SReal dt, sofa::core::MultiVecCoordId xResult, sofa::core::MultiVecDerivId vResult) override;
+
+    SingleLink<NewtonRaphsonSolver, core::behavior::BaseIntegrationMethod, BaseLink::FLAG_STRONGLINK> l_integrationMethod;
+
+    Data<unsigned int> d_maxNbIterationsNewton;
+    Data<SReal> d_relativeSuccessiveStoppingThreshold;
+    Data<SReal> d_relativeInitialStoppingThreshold;
+    Data<SReal> d_absoluteResidualStoppingThreshold;
+    Data<SReal> d_relativeEstimateDifferenceThreshold;
+    Data<SReal> d_absoluteEstimateDifferenceThreshold;
+    Data<unsigned int> d_maxNbIterationsLineSearch;
+    Data<SReal> d_lineSearchCoefficient;
+
+    MAKE_SELECTABLE_ITEMS(NewtonStatus,
+        sofa::helper::Item{"Undefined", "The solver has not been called yet"},
+        sofa::helper::Item{"Running", "The solver is still running and/or did not finish"},
+        sofa::helper::Item{"ConvergedEquilibrium", "Converged: the iterations did not start because the system is already at equilibrium"},
+        sofa::helper::Item{"DivergedLineSearch", "Diverged: line search failed"},
+        sofa::helper::Item{"DivergedMaxIterations", "Diverged: Reached the maximum number of iterations"},
+        sofa::helper::Item{"ConvergedResidualSuccessiveRatio", "Converged: Residual successive ratio is smaller than the threshold"},
+        sofa::helper::Item{"ConvergedResidualInitialRatio", "Converged: Residual initial ratio is smaller than the threshold"},
+        sofa::helper::Item{"ConvergedAbsoluteResidual", "Converged: Absolute residual is smaller than the threshold"},
+        sofa::helper::Item{"ConvergedRelativeEstimateDifference", "Converged: Relative estimate difference is smaller than the threshold"},
+        sofa::helper::Item{"ConvergedAbsoluteEstimateDifference", "Converged: Absolute estimate difference is smaller than the threshold"});
+
+    Data<NewtonStatus> d_status;
+
+protected:
+    NewtonRaphsonSolver();
+
+    core::behavior::MultiVecDeriv m_linearSystemSolution;
+
+    void computeRightHandSide(const core::ExecParams* params, SReal dt,
+                      core::MultiVecDerivId force,
+                      core::MultiVecDerivId b,
+                      core::MultiVecDerivId velocity_i,
+                      core::MultiVecCoordId position_i) const;
+
+
+    SReal computeResidual(const core::ExecParams* params, sofa::simulation::common::MechanicalOperations& mop, SReal dt,
+        core::MultiVecDerivId force,
+        core::MultiVecDerivId oldVelocity, core::MultiVecDerivId newVelocity);
+};
+
+}

Sofa/Component/ODESolver/Backward/src/sofa/component/odesolver/backward/init.cpp

@@ -30,6 +30,7 @@ extern void registerEulerImplicitSolver(sofa::core::ObjectFactory* factory);
 extern void registerNewmarkImplicitSolver(sofa::core::ObjectFactory* factory);
 extern void registerStaticSolver(sofa::core::ObjectFactory* factory);
 extern void registerVariationalSymplecticSolver(sofa::core::ObjectFactory* factory);
+extern void registerNewtonRaphsonSolver(sofa::core::ObjectFactory* factory);
 
 extern "C" {
     SOFA_EXPORT_DYNAMIC_LIBRARY void initExternalModule();
@@ -59,6 +60,7 @@ void registerObjects(sofa::core::ObjectFactory* factory)
     registerNewmarkImplicitSolver(factory);
     registerStaticSolver(factory);
     registerVariationalSymplecticSolver(factory);
+    registerNewtonRaphsonSolver(factory);
 }
 
 void init()

Sofa/Component/ODESolver/Backward/tests/CMakeLists.txt

@@ -7,6 +7,7 @@ set(SOURCE_FILES
     EulerImplicitSolverStatic_test.cpp
     EulerImplicitSolver_withDamping_test.cpp
     NewmarkImplicitSolverDynamic_test.cpp
+    NewtonRaphsonSolver_test.cpp
     StaticSolver_test.cpp
     SpringSolverDynamic_test.cpp
     VariationalSymplecticExplicitSolverDynamic_test.cpp

Sofa/Component/ODESolver/Backward/tests/NewtonRaphsonSolver_test.cpp

@@ -0,0 +1,212 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2006 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* This program is free software; you can redistribute it and/or modify it     *
+* under the terms of the GNU Lesser General Public License as published by    *
+* the Free Software Foundation; either version 2.1 of the License, or (at     *
+* your option) any later version.                                             *
+*                                                                             *
+* This program is distributed in the hope that it will be useful, but WITHOUT *
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
+* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
+* for more details.                                                           *
+*                                                                             *
+* You should have received a copy of the GNU Lesser General Public License    *
+* along with this program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: [email protected]                             *
+******************************************************************************/
+#include <sofa/component/odesolver/backward/NewtonRaphsonSolver.h>
+#include <sofa/component/statecontainer/MechanicalObject.h>
+#include <sofa/simpleapi/SimpleApi.h>
+#include <sofa/testing/BaseSimulationTest.h>
+#include <sofa/testing/NumericTest.h>
+
+
+namespace sofa
+{
+
+struct NewtonRaphsonParameters
+{
+    SReal relativeSuccessiveStoppingThreshold {};
+    SReal relativeInitialStoppingThreshold {};
+    SReal absoluteResidualStoppingThreshold {};
+    SReal relativeEstimateDifferenceThreshold {};
+    SReal absoluteEstimateDifferenceThreshold {};
+    unsigned int maxNbIterationsLineSearch = 1;
+    unsigned int maxNbIterationsNewton = 1;
+
+    void apply(component::odesolver::backward::NewtonRaphsonSolver* solver) const
+    {
+        solver->d_relativeSuccessiveStoppingThreshold.setValue(relativeSuccessiveStoppingThreshold);
+        solver->d_relativeInitialStoppingThreshold.setValue(relativeInitialStoppingThreshold);
+        solver->d_absoluteResidualStoppingThreshold.setValue(absoluteResidualStoppingThreshold);
+        solver->d_relativeEstimateDifferenceThreshold.setValue(relativeEstimateDifferenceThreshold);
+        solver->d_absoluteEstimateDifferenceThreshold.setValue(absoluteEstimateDifferenceThreshold);
+        solver->d_maxNbIterationsLineSearch.setValue(maxNbIterationsLineSearch);
+        solver->d_maxNbIterationsNewton.setValue(maxNbIterationsNewton);
+    }
+};
+
+struct NewtonRaphsonTest : public testing::BaseSimulationTest
+{
+    SceneInstance m_scene{};
+
+    component::odesolver::backward::NewtonRaphsonSolver* m_solver { nullptr };
+    component::statecontainer::MechanicalObject<defaulttype::Vec1Types>* m_state { nullptr };
+
+    void onSetUp() override
+    {
+        sofa::simpleapi::importPlugin("Sofa.Component.LinearSolver.Direct");
+        sofa::simpleapi::importPlugin("Sofa.Component.ODESolver.Backward");
+        sofa::simpleapi::importPlugin("Sofa.Component.ODESolver.Integration");
+        sofa::simpleapi::importPlugin("Sofa.Component.StateContainer");
+        sofa::simpleapi::importPlugin("Sofa.Component.Topology.Container.Dynamic");
+
+        sofa::simpleapi::createObject(m_scene.root, "DefaultAnimationLoop");
+        sofa::simpleapi::createObject(m_scene.root, "BDF1");
+        auto solverObject = sofa::simpleapi::createObject(m_scene.root, "NewtonRaphsonSolver", {{"printLog", "true"}});
+        m_solver = dynamic_cast<component::odesolver::backward::NewtonRaphsonSolver*>(solverObject.get());
+        sofa::simpleapi::createObject(m_scene.root, "EigenSimplicialLDLT", {{"template", "CompressedRowSparseMatrix"}});
+        auto stateObject = sofa::simpleapi::createObject(m_scene.root, "MechanicalObject", {{"template", "Vec1"}, {"position", "1"}, {"rest_position", "0"}});
+        m_state = dynamic_cast<component::statecontainer::MechanicalObject<defaulttype::Vec1Types>*>(stateObject.get());
+        sofa::simpleapi::createObject(m_scene.root, "PointSetTopologyContainer", {{"name", "topology"}, {"position", "1"}});
+
+        ASSERT_NE(nullptr, m_solver);
+        ASSERT_NE(nullptr, m_state);
+    }
+
+    void noForce()
+    {
+        m_scene.initScene();
+        EXPECT_EQ(m_solver->d_status.getValue(), component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus("Undefined"));
+
+        m_scene.simulate(0.01_sreal);
+        EXPECT_EQ(m_solver->d_status.getValue(), component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus("ConvergedEquilibrium"));
+    }
+
+    void gravity(const NewtonRaphsonParameters& params, const component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus& expectedStatus)
+    {
+        static constexpr SReal gravity = -1;
+        m_scene.root->setGravity({gravity, 0, 0});
+
+        sofa::simpleapi::importPlugin("Sofa.Component.Mass");
+        sofa::simpleapi::createObject(m_scene.root, "UniformMass", {{"totalMass", "1"}, {"topology", "@topology"}, {"printLog", "true"}});
+        params.apply(m_solver);
+
+        m_scene.initScene();
+        EXPECT_EQ(m_solver->d_status.getValue(), component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus("Undefined"));
+
+        auto previousVelocity = m_state->read(sofa::core::vec_id::read_access::velocity)->getValue();
+        auto previousPosition = m_state->read(sofa::core::vec_id::read_access::position)->getValue();
+
+        static constexpr SReal dt = 0.1_sreal;
+        m_scene.simulate(dt);
+        EXPECT_EQ(m_solver->d_status.getValue(), expectedStatus);
+
+        const auto& velocity = m_state->read(sofa::core::vec_id::read_access::velocity)->getValue();
+        const auto& position = m_state->read(sofa::core::vec_id::read_access::position)->getValue();
+
+        EXPECT_FLOATINGPOINT_EQ(velocity[0].x(), gravity * dt)
+        EXPECT_FLOATINGPOINT_EQ(position[0].x(), previousPosition[0].x() + velocity[0].x() * dt);
+
+        for (unsigned int i = 1; i < 10; ++i)
+        {
+            previousVelocity = m_state->read(sofa::core::vec_id::read_access::velocity)->getValue();
+            previousPosition = m_state->read(sofa::core::vec_id::read_access::position)->getValue();
+
+            m_scene.simulate(dt);
+            EXPECT_FLOATINGPOINT_EQ(velocity[0].x(), previousVelocity[0].x() + gravity * dt);
+
+            EXPECT_FLOATINGPOINT_EQ(position[0].x(), previousPosition[0].x() + velocity[0].x() * dt);
+        }
+    }
+
+    void spring()
+    {
+        sofa::simpleapi::importPlugin("Sofa.Component.SolidMechanics.Spring");
+        sofa::simpleapi::createObject(m_scene.root, "RestShapeSpringsForceField", {{"points", "0"}});
+
+        m_scene.initScene();
+        EXPECT_EQ(m_solver->d_status.getValue(), component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus("Undefined"));
+
+        m_scene.simulate(0.01_sreal);
+    }
+
+};
+
+TEST_F(NewtonRaphsonTest, noForce)
+{
+    this->noForce();
+}
+
+TEST_F(NewtonRaphsonTest, gravity_noStopping)
+{
+    // no stopping criteria is defined so it cannot converge
+    this->gravity(
+        NewtonRaphsonParameters{},
+        component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus("DivergedMaxIterations"));
+}
+
+TEST_F(NewtonRaphsonTest, gravity_relativeSuccessiveStopping)
+{
+    NewtonRaphsonParameters params;
+    params.relativeSuccessiveStoppingThreshold = 1e-15_sreal;
+
+    this->gravity(
+        params,
+        component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus("ConvergedResidualSuccessiveRatio"));
+}
+
+TEST_F(NewtonRaphsonTest, gravity_relativeInitialStopping)
+{
+    NewtonRaphsonParameters params;
+    params.relativeInitialStoppingThreshold = 1e-15_sreal;
+
+    this->gravity(
+        params,
+        component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus("ConvergedResidualInitialRatio"));
+}
+
+TEST_F(NewtonRaphsonTest, gravity_absoluteResidualStopping)
+{
+    NewtonRaphsonParameters params;
+    params.absoluteResidualStoppingThreshold = 1e-15_sreal;
+
+    this->gravity(
+        params,
+        component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus("ConvergedAbsoluteResidual"));
+}
+
+TEST_F(NewtonRaphsonTest, gravity_relativeEstimateDifferenceStopping)
+{
+    NewtonRaphsonParameters params;
+    params.relativeEstimateDifferenceThreshold = 1e-15_sreal;
+
+    // considered as diverged because this threshold requires more than one iteration
+    this->gravity(
+        params,
+        component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus("DivergedMaxIterations"));
+}
+
+TEST_F(NewtonRaphsonTest, gravity_absoluteEstimateDifferenceStopping)
+{
+    NewtonRaphsonParameters params;
+    params.absoluteEstimateDifferenceThreshold = 1e-15_sreal;
+
+    // considered as diverged because this threshold requires more than one iteration
+    this->gravity(
+        params,
+        component::odesolver::backward::NewtonRaphsonSolver::NewtonStatus("DivergedMaxIterations"));
+}
+
+TEST_F(NewtonRaphsonTest, spring)
+{
+    this->spring();
+}
+
+
+}

Sofa/Component/ODESolver/CMakeLists.txt

@@ -4,7 +4,7 @@ project(Sofa.Component.ODESolver LANGUAGES CXX)
 set(SOFACOMPONENTODESOLVER_SOURCE_DIR "src/sofa/component/odesolver")
 
 sofa_add_subdirectory_modules(SOFACOMPONENTODESOLVER_TARGETS
-    DIRECTORIES Forward Backward
+    DIRECTORIES Forward Backward Integration
 )
 
 set(HEADER_FILES

Sofa/Component/ODESolver/Integration/CMakeLists.txt

@@ -0,0 +1,35 @@
+cmake_minimum_required(VERSION 3.22)
+project(Sofa.Component.ODESolver.Integration LANGUAGES CXX)
+
+set(SOFACOMPONENTODESOLVERINTEGRATION_SOURCE_DIR "src/sofa/component/odesolver/integration")
+
+set(HEADER_FILES
+    ${SOFACOMPONENTODESOLVERINTEGRATION_SOURCE_DIR}/config.h.in
+    ${SOFACOMPONENTODESOLVERINTEGRATION_SOURCE_DIR}/init.h
+
+    ${SOFACOMPONENTODESOLVERINTEGRATION_SOURCE_DIR}/BDF1.h
+)
+
+set(SOURCE_FILES
+    ${SOFACOMPONENTODESOLVERINTEGRATION_SOURCE_DIR}/init.cpp
+
+    ${SOFACOMPONENTODESOLVERINTEGRATION_SOURCE_DIR}/BDF1.cpp
+)
+
+sofa_find_package(Sofa.Simulation.Core REQUIRED)
+
+add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES})
+target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Simulation.Core)
+
+sofa_create_package_with_targets(
+    PACKAGE_NAME ${PROJECT_NAME}
+    PACKAGE_VERSION ${Sofa_VERSION}
+    TARGETS ${PROJECT_NAME} AUTO_SET_TARGET_PROPERTIES
+    INCLUDE_SOURCE_DIR "src"
+    INCLUDE_INSTALL_DIR "${PROJECT_NAME}"
+)
+
+#cmake_dependent_option(SOFA_COMPONENT_ODESOLVER_INTEGRATION_BUILD_TESTS "Compile the automatic tests" ON "SOFA_BUILD_TESTS OR NOT DEFINED SOFA_BUILD_TESTS" OFF)
+#if(SOFA_COMPONENT_ODESOLVER_INTEGRATION_BUILD_TESTS)
+#    add_subdirectory(tests)
+#endif()

Sofa/Component/ODESolver/Integration/Sofa.Component.ODESolver.IntegrationConfig.cmake.in

@@ -0,0 +1,12 @@
+# CMake package configuration file for the @PROJECT_NAME@ module
+
+@PACKAGE_GUARD@
+@PACKAGE_INIT@
+
+find_package(Sofa.Simulation.Core QUIET REQUIRED)
+
+if(NOT TARGET @PROJECT_NAME@)
+    include("${CMAKE_CURRENT_LIST_DIR}/@[email protected]")
+endif()
+
+check_required_components(@PROJECT_NAME@)