Skip to content
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.

Sign up for GitHub

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

Add callback mechanism to Newton-Raphson solvers #107

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add callback mechanism to Newton-Raphson solvers #107

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
cmake_minimum_required(VERSION 3.12)
project(Caribou VERSION 21.06.00)
project(Caribou VERSION 22.06.99)

# Policies
cmake_policy(SET CMP0072 NEW)
Expand Down
26 changes: 19 additions & 7 deletions src/SofaCaribou/Ode/NewtonRaphsonSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,6 @@ void NewtonRaphsonSolver::solve(const ExecParams *params, SReal dt, MultiVecCoor
}

// Local variables used for the iterations
unsigned n_it=0;
double dx_squared_norm, du_squared_norm, R_squared_norm = 0;
const auto squared_residual_threshold = residual_tolerance_threshold*residual_tolerance_threshold;
const auto squared_correction_threshold = correction_tolerance_threshold*correction_tolerance_threshold;
Expand Down Expand Up @@ -248,16 +247,20 @@ void NewtonRaphsonSolver::solve(const ExecParams *params, SReal dt, MultiVecCoor
// # Newton iterations #
// ###########################################################################

while (not converged and n_it < newton_iterations) {
p_current_iteration = 0;
while (not converged and p_current_iteration < newton_iterations) {
sofa::helper::ScopedAdvancedTimer step_timer ("NewtonStep");
t = steady_clock::now();

trigger_event(Event::ITERATION_BEGIN);

// Part 1. Assemble the system matrix.
{
sofa::helper::ScopedAdvancedTimer _t_("MBKBuild");
p_A->clear();
this->assemble_system_matrix(mechanical_parameters, accessor, p_A.get());
linear_solver->set_system_matrix(p_A.get());
trigger_event(Event::MATRIX_ASSEMBLED);
}

// Part 2. Analyze the pattern of the matrix in order to compute a permutation matrix.
Expand All @@ -281,6 +284,8 @@ void NewtonRaphsonSolver::solve(const ExecParams *params, SReal dt, MultiVecCoor
break;
}

trigger_event(Event::MATRIX_ANALYZED);

p_has_already_analyzed_the_pattern = true;
}
}
Expand All @@ -293,6 +298,8 @@ void NewtonRaphsonSolver::solve(const ExecParams *params, SReal dt, MultiVecCoor
diverged = true;
break;
}

trigger_event(Event::MATRIX_FACTORIZED);
}

// Part 4. Solve the unknown increment.
Expand All @@ -303,12 +310,14 @@ void NewtonRaphsonSolver::solve(const ExecParams *params, SReal dt, MultiVecCoor
diverged = true;
break;
}
trigger_event(Event::INCREMENT_SOLVED);
}

// Part 5. Propagating the solution increment and update geometry.
{
sofa::helper::ScopedAdvancedTimer _t_("PropagateDx");
this->propagate_solution_increment(mechanical_parameters, accessor, p_DX.get(), x_id, v_id, dx_id);
trigger_event(Event::INCREMENT_PROPAGATED);
}

// The next two parts are only necessary when doing more than one Newton iteration
Expand All @@ -323,6 +332,7 @@ void NewtonRaphsonSolver::solve(const ExecParams *params, SReal dt, MultiVecCoor
sofa::helper::AdvancedTimer::stepBegin("UpdateResidual");
R_squared_norm = SofaCaribou::Algebra::dot(p_F.get(), p_F.get());
sofa::helper::AdvancedTimer::stepEnd("UpdateResidual");
trigger_event(Event::RESIDUAL_UPDATED);
}

// Part 8. Compute the updated displacement residual.
Expand All @@ -341,11 +351,13 @@ void NewtonRaphsonSolver::solve(const ExecParams *params, SReal dt, MultiVecCoor

p_squared_residuals.emplace_back(R_squared_norm);

trigger_event(Event::ITERATION_END);

// We completed one iteration, increment the counter
n_it++;
p_current_iteration++;

if( print_log ) {
info << "Newton iteration #" << std::left << std::setw(5) << n_it
info << "Newton iteration #" << std::left << std::setw(5) << p_current_iteration
<< std::scientific
<< " |R| = " << std::setw(12) << sqrt(R_squared_norm)
<< " |R|/|R0| = " << std::setw(12) << sqrt(R_squared_norm / p_squared_residuals[0])
Expand Down Expand Up @@ -405,9 +417,9 @@ void NewtonRaphsonSolver::solve(const ExecParams *params, SReal dt, MultiVecCoor
vop.v_clear(dx_id);
} // End while (not converged and not diverged and n_it < newton_iterations)

n_it--; // Reset to the actual index of the last iteration completed
p_current_iteration--; // Reset to the actual index of the last iteration completed

if (not converged and not diverged and n_it == (newton_iterations-1)) {
if (not converged and not diverged and p_current_iteration == (newton_iterations-1)) {
if (print_log) {
info << "[DIVERGED] The number of Newton iterations reached the maximum of " << newton_iterations << " iterations" << ".\n";
}
Expand All @@ -416,7 +428,7 @@ void NewtonRaphsonSolver::solve(const ExecParams *params, SReal dt, MultiVecCoor
d_converged.setValue(converged);

sofa::helper::AdvancedTimer::valSet("has_converged", converged ? 1 : 0);
sofa::helper::AdvancedTimer::valSet("nb_iterations", n_it+1);
sofa::helper::AdvancedTimer::valSet("nb_iterations", p_current_iteration+1);
}

void NewtonRaphsonSolver::init() {
Expand Down
83 changes: 83 additions & 0 deletions src/SofaCaribou/Ode/NewtonRaphsonSolver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@ DISABLE_ALL_WARNINGS_BEGIN
DISABLE_ALL_WARNINGS_END

#include <memory>
#include <functional>
#include <array>

namespace SofaCaribou::ode {

Expand Down Expand Up @@ -63,6 +65,35 @@ class NewtonRaphsonSolver : public sofa::core::behavior::OdeSolver {
ALWAYS
};

/**
* Events associated to registered callback function for different time in the NR iteration process.
*/
enum class Event : unsigned int {
/** Event triggered at the very beginning of the Newton iteration. */
ITERATION_BEGIN = 0,

/** Event triggered after the system matrix has been assembled. */
MATRIX_ASSEMBLED,

/** Event triggered after the system matrix has been analyzed. */
MATRIX_ANALYZED,

/** Event triggered after the system matrix has been factorized. */
MATRIX_FACTORIZED,

/** Event triggered after the system solution vector has been solved. */
INCREMENT_SOLVED,

/** Event triggered after the system solution vector has been propagated through mappings. */
INCREMENT_PROPAGATED,

/** Event triggered after the system force vector has been updated. */
RESIDUAL_UPDATED,

/** Event triggered at the very end of the Newton iteration. */
ITERATION_END,
};


NewtonRaphsonSolver();

Expand Down Expand Up @@ -92,6 +123,45 @@ class NewtonRaphsonSolver : public sofa::core::behavior::OdeSolver {

void set_pattern_analysis_strategy(const PatternAnalysisStrategy & strategy);

/** Get the complete force vector of the system */
auto F() const -> const sofa::defaulttype::BaseVector * {
return p_F.get();
}

/** Get the complete increment vector (solution) of the system */
auto dx() const -> const sofa::defaulttype::BaseVector * {
return p_DX.get();
}

/** Get the complete matrix of the linearized part of system */
auto A() const -> const sofa::defaulttype::BaseMatrix * {
return p_A.get();
}

/** Get the current newton iteration within the time step solve call. First iteration is 1 (index starts from 1). */
auto current_iteration() const -> unsigned int {
return p_current_iteration+1;
}

/**
* Register a callback function to be called at the specific NR event (see NewtonRaphsonSolver::Event for the list
* of events).
*
* @example
* @code
* using NewtonRaphson::Event;
* solver->register_callback(MATRIX_ASSEMBLED, [](const NewtonRaphsonSolver & solver){
* std::cout << "The system tangent matrix has been assembled in NR iteration #" << solver.current_iteration << '\n';
* });
* @endcode
* @param event
*/
inline void register_callback(const Event & event, const std::function<void(const NewtonRaphsonSolver &)>& fct) {
const auto event_id = static_cast<unsigned int>(event);
caribou_assert(event_id >= 0 and event_id <= static_cast<unsigned int>(Event::ITERATION_END));
p_iteration_callbacks[event_id].push_back(fct);
}

private:

/**
Expand Down Expand Up @@ -164,6 +234,13 @@ class NewtonRaphsonSolver : public sofa::core::behavior::OdeSolver {
sofa::core::MultiVecDerivId & v_id,
sofa::core::MultiVecDerivId & dx_id) = 0;

void trigger_event(const Event & event) const {
const auto event_id = static_cast<unsigned int>(event);
for (const auto & cb : p_iteration_callbacks[event_id]) {
cb(*this);
}
}

protected:
/** Check that the linked linear solver is not null and that it implements the SofaCaribou::solver::LinearSolver interface */

Expand Down Expand Up @@ -207,5 +284,11 @@ class NewtonRaphsonSolver : public sofa::core::behavior::OdeSolver {

/// Either or not the pattern of the system matrix was analyzed at the beginning of the simulation
bool p_has_already_analyzed_the_pattern = false;

/// Current newton iteration within a time step solve call
unsigned int p_current_iteration;

/// List of callback function to be called at each events of a NR iteration
std::array<std::vector<std::function<void(const NewtonRaphsonSolver &)>>, static_cast<unsigned int>(Event::ITERATION_END)+1> p_iteration_callbacks;
};
}
25 changes: 21 additions & 4 deletions src/SofaCaribou/Python/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ set(HEADER_FILES
Forcefield/HyperelasticForcefield.h
Mass/CaribouMass.h
Ode/LegacyStaticODESolver.h
Ode/NewtonRaphsonSolver.h
Ode/StaticODESolver.h
Solver/ConjugateGradientSolver.h
Solver/LDLTSolver.h
Expand All @@ -22,6 +23,7 @@ set(SOURCE_FILES
Forcefield/HyperelasticForcefield.cpp
Mass/CaribouMass.cpp
Ode/LegacyStaticODESolver.cpp
Ode/NewtonRaphsonSolver.cpp
Ode/StaticODESolver.cpp
Solver/ConjugateGradientSolver.cpp
Solver/LDLTSolver.cpp
Expand All @@ -38,6 +40,7 @@ set(PYTHON_FILES
set(PYTHON_TEST_FILES
pytest/SofaCaribou_Forcefield_HyperelasticForcefield.py
pytest/SofaCaribou_Mass_CaribouMass.py
pytest/SofaCaribou_StaticSolver.py
)

find_package(SofaPython3 REQUIRED)
Expand All @@ -64,17 +67,31 @@ caribou_add_python_module(SofaCaribou
)

list(APPEND target_rpath
"$ORIGIN/../../../../../../lib"
"$loader_path/../../../../../../lib"
"$ORIGIN/../../../../lib"
"$loader_path/../../../../lib"
"$ORIGIN/../../../../../../lib"
"$$ORIGIN/../../../../../../lib"
"$loader_path/../../../../../../lib"
"@loader_path/../../../../../../lib"
"$executable_path/../../../../../../lib"
"@executable_path/../../../../../../lib"
"$ORIGIN/../../../../lib"
"$$ORIGIN/../../../../lib"
"$loader_path/../../../../lib"
"@loader_path/../../../../lib"
"$executable_path/../../../../lib"
"@executable_path/../../../../lib"
"$loader_path/../Caribou"
"@loader_path/../Caribou"
)
set_target_properties(
${PROJECT_NAME}
PROPERTIES
INSTALL_RPATH "${target_rpath}"
)


# Make the target relocatable
set_target_properties(${PROJECT_NAME} PROPERTIES RELOCATABLE_INSTALL_DIR "plugins/SofaCaribou")
set_target_properties(${PROJECT_NAME} PROPERTIES EXPORT_PROPERTIES "RELOCATABLE_INSTALL_DIR")
if (SOFA_VERSION VERSION_GREATER_EQUAL "20.12.99")
# Set RPATH towards relocatable dependencies
sofa_auto_set_target_rpath(TARGETS ${PROJECT_NAME} RELOCATABLE "plugins")
Expand Down
Loading