Commit 6aefe4a3 authored by Ricardo Ortiz's avatar Ricardo Ortiz
Browse files

ENH: Streamline and refactor solvers

and their tests. Streamline and refactor
time integrators and their tests.
parent 1c01c312
Pipeline #4819 passed with stage
......@@ -106,7 +106,7 @@ int main()
auto planeObject = std::make_shared<StaticSceneObject>();
planeObject->setModel(planeModel);
scene->addSceneObject(planeObject);
// scene->addSceneObject(planeObject);
// Setup Scene lighting
auto light1 = Light::getDefaultLighting();
......
......@@ -8,7 +8,7 @@ simmedtk_add_library(${MODULE_NAME}
DirectLinearSolver.cpp
ForwardGaussSeidel.cpp
ForwardSOR.cpp
InexactNewton.cpp
NewtonMethod.cpp
IterativeLinearSolver.cpp
NonLinearSolver.cpp
PUBLIC_HEADERS
......@@ -19,7 +19,7 @@ simmedtk_add_library(${MODULE_NAME}
DirectLinearSolver.h
ForwardGaussSeidel.h
ForwardSOR.h
InexactNewton.h
NewtonMethod.h
IterativeLinearSolver.h
LinearSolver.h
LinearSolver.hpp
......
......@@ -22,7 +22,6 @@
//---------------------------------------------------------------------------
#include "Solvers/ConjugateGradient.h"
#include <iostream>
ConjugateGradient::ConjugateGradient(const core::SparseMatrixd &A,
const core::Vectord &rhs) : solver(A)
......@@ -46,7 +45,6 @@ void ConjugateGradient::solve(core::Vectord &x)
{
if(this->linearSystem)
{
this->solver.compute(this->linearSystem->getMatrix());
x = this->solver.solve(this->linearSystem->getRHSVector());
}
......@@ -57,7 +55,6 @@ void ConjugateGradient::solve(core::Vectord &x)
void ConjugateGradient::solve(core::Vectord &x, double tolerance)
{
this->setTolerance(tolerance);
this->solver.setTolerance(tolerance);
this->solve(x);
}
......@@ -80,3 +77,10 @@ void ConjugateGradient::setMaximumIterations(const size_t maxIter)
this->maxIterations = maxIter;
this->solver.setMaxIterations(maxIter);
}
//---------------------------------------------------------------------------
void ConjugateGradient::setSystem(std::shared_ptr<LinearSystem<core::SparseMatrixd>> newSystem)
{
LinearSolver<core::SparseMatrixd>::setSystem(newSystem);
this->solver.compute(this->linearSystem->getMatrix());
}
......@@ -96,6 +96,13 @@ public:
///
void setMaximumIterations(const size_t maxIter) override;
///
/// \brief Sets the system. System of linear equations.
///
/// \param newSystem Linear system replacement.
///
void setSystem(std::shared_ptr<LinearSystemType> newSystem) override;
private:
///> Pointer to the Eigen's Conjugate gradient solver
Eigen::ConjugateGradient<core::SparseMatrixd> solver;
......
......@@ -91,7 +91,7 @@ public:
///
/// \brief Default constructor/destructor
///
DirectLinearSolver() = delete;
DirectLinearSolver() = default;
~DirectLinearSolver() = default;
///
......
// This file is part of the SimMedTK project.
// Copyright (c) Center for Modeling, Simulation, and Imaging in Medicine,
// Rensselaer Polytechnic Institute
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
//---------------------------------------------------------------------------
//
// Authors:
//
// Contact:
//---------------------------------------------------------------------------
#include "InexactNewton.h"
#include "Solvers/IterativeLinearSolver.h"
#include "Solvers/ConjugateGradient.h"
#include <iostream>
InexactNewton::InexactNewton():
linearSolver(std::make_shared<ConjugateGradient>()),
forcingTerm(0.9),
absoluteTolerance(1e-3),
relativeTolerance(1e-6),
gamma(0.9),
etaMax(0.9),
maxIterations(50),
useArmijo(true) {}
//---------------------------------------------------------------------------
void InexactNewton::solve(core::Vectord &x)
{
if(!this->linearSolver)
{
return;
}
// Compute norms, set tolerances and other temporaries
double fnorm = this->nonLinearSystem->eval(x).norm();
double stopTolerance = this->absoluteTolerance + this->relativeTolerance * fnorm;
this->linearSolver->setTolerance(stopTolerance);
core::Vectord dx = x;
for(size_t i = 0; i < this->maxIterations; ++i)
{
if(fnorm < stopTolerance)
{
return;
}
this->updateJacobian(x);
this->linearSolver->solve(dx);
this->updateIterate(-dx,x);
double newNorm = fnorm;
if(this->useArmijo)
{
newNorm = this->armijo(dx, x);
}
if(this->forcingTerm > 0 && newNorm > stopTolerance)
{
double ratio = newNorm / fnorm; // Ratio of successive residual norms
this->updateForcingTerm(ratio, stopTolerance, fnorm);
// Reset tolerance in the linear solver according to the new forcing term
// to avoid over solving of the system.
if(this->linearSolver)
{
this->linearSolver->setTolerance(this->forcingTerm);
}
}
fnorm = newNorm;
}
}
//---------------------------------------------------------------------------
void InexactNewton::updateForcingTerm(const double ratio,
const double stopTolerance,
const double fnorm)
{
double eta = this->gamma * ratio * ratio;
double forcingTermSqr = this->forcingTerm * this->forcingTerm;
// Save guard to prevent the forcing term to become too small for far away iterates
if(this->gamma * forcingTermSqr > .1)
{
eta = std::max(eta, this->gamma * forcingTermSqr);
}
this->forcingTerm = std::max(std::min(eta, this->etaMax), 0.5 * stopTolerance / fnorm);
}
//---------------------------------------------------------------------------
void InexactNewton::setLinearSolver(std::shared_ptr< InexactNewton::LinearSolverType > newLinearSolver)
{
this->linearSolver = newLinearSolver;
}
//---------------------------------------------------------------------------
std::shared_ptr< InexactNewton::LinearSolverType > InexactNewton::getLinearSolver() const
{
return this->linearSolver;
}
//---------------------------------------------------------------------------
void InexactNewton::updateJacobian(const core::Vectord &x)
{
// Evaluate the Jacobian and sets the matrix
if(!this->jacobian)
{// TODO: Print message or log error.
return;
}
const core::SparseMatrixd &jacobianMatrix = this->jacobian(x);
const core::Vectord &f = this->nonLinearSystem->getFunctionValue();
if(!(jacobianMatrix.innerSize() > 0))
{// TODO: Print message and/or log error.
return;
}
auto linearSystem = std::make_shared<LinearSolverType::LinearSystemType>(
jacobianMatrix, f);
this->linearSolver->setSystem(linearSystem);
}
//---------------------------------------------------------------------------
void InexactNewton::setJacobian(const NonLinearSolver::JacobianType &newJacobian)
{
this->jacobian = newJacobian;
}
//---------------------------------------------------------------------------
const NonLinearSolver::JacobianType &InexactNewton::getJacobian() const
{
return this->jacobian;
}
//---------------------------------------------------------------------------
void InexactNewton::setAbsoluteTolerance(const double aTolerance)
{
this->absoluteTolerance = aTolerance;
}
//---------------------------------------------------------------------------
double InexactNewton::getAbsoluteTolerance() const
{
return this->absoluteTolerance;
}
// This file is part of the SimMedTK project.
// Copyright (c) Center for Modeling, Simulation, and Imaging in Medicine,
// Rensselaer Polytechnic Institute
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
//---------------------------------------------------------------------------
//
// Authors:
//
// Contact:
//---------------------------------------------------------------------------
#ifndef INEXACT_NEWTON_H
#define INEXACT_NEWTON_H
#include "Solvers/NonLinearSolver.h"
#include "Solvers/LinearSolver.h"
///
/// @brief Inexact Newton method. This version of the newton method is based on the work
/// Tim Kelly and others at NC State University. The algorithm is globally convergent
/// in the sense that for any initial iterate the iteration either converge to a
/// root of F or fails. Global convergence is achieved using a line search sub-process
/// and the Armijo rule.
///
class InexactNewton : public NonLinearSolver
{
public:
using LinearSolverType = LinearSolver<core::SparseMatrixd>;
public:
///
/// \brief Constructors/Destructors.
///
InexactNewton();
~InexactNewton() = default;
InexactNewton(const InexactNewton &other) = delete;
InexactNewton &operator=(const InexactNewton &other) = delete;
///
/// \brief Solve the non linear system of equations for a contant Jacobian matrix.
///
/// \param x Current iterate.
///
void solve(core::Vectord &x) override;
///
/// \brief Update forcing term according to Eisenstat-Walker criteria.
/// \see http://softlib.rice.edu/pub/CRPC-TRs/reports/CRPC-TR94463.pdf.
///
/// \param ratio Ratio of consecutive function norms.
/// \param stopTolerance Linear combinarion of relative and absolute tolerances.
/// \param fnorm Current function norm.
///
void updateForcingTerm(const double ratio,
const double stopTolerance,
const double fnorm);
///
/// \brief Set LinearSolver.
///
/// \param newLinearSolver Linear solver pointer.
///
void setLinearSolver(std::shared_ptr<LinearSolverType> newLinearSolver);
///
/// \brief Get LinearSolver.
///
std::shared_ptr<LinearSolverType> getLinearSolver() const;
///
/// \brief Update jacobians.
///
/// \param x Current iterate.
///
void updateJacobian(const core::Vectord &x);
///
/// \brief Set Jacobian. This sets the jacobian function, used to compute the jacobian matrix.
///
/// \param newJacobian Jacobian function.
///
void setJacobian(const JacobianType &newJacobian);
///
/// \brief Get Jacobian. Returns the Jacobian function.
///
const JacobianType &getJacobian() const;
///
/// \brief Set JacobianMatrix.
///
/// \param newJacobianMatrix Replaces the stored jacobian matrix.
///
void setJacobianMatrix(const core::SparseMatrixd &newJacobianMatrix);
///
/// \brief Get JacobianMatrix. Returns jacobian matrix.
///
const core::SparseMatrixd &getJacobianMatrix() const;
///
/// \brief Get JacobianMatrix. Returns jacobian matrix.
///
core::SparseMatrixd &getJacobianMatrix();
///
/// \brief Set AbsoluteTolerance.
///
/// \param aTolerance New absolute tolerance value.
///
void setAbsoluteTolerance(const double aTolerance);
///
/// \brief Get AbsoluteTolerance. Returns current tolerance value.
///
double getAbsoluteTolerance() const;
///
/// \brief Set RelativeTolerance.
///
/// \param newRelativeTolerance New relative tolerance value
///
void setRelativeTolerance(const double newRelativeTolerance)
{
this->relativeTolerance = newRelativeTolerance;
}
///
/// \brief Get RelativeTolerance. Returns current relative tolerance value.
///
double getRelativeTolerance() const
{
return this->relativeTolerance;
}
///
/// \brief Set Gamma. Sets the gamma parameter used to update the forcing term.
///
/// \param newGamma New gamma value.
///
void setGamma(const double newGamma)
{
this->gamma = newGamma;
}
///
/// \brief Get Gamma. Returns current gamma value.
///
double getGamma() const
{
return this->gamma;
}
///
/// \brief Set EtaMax. Maximum error tolerance for residual in inner
/// iteration. The inner iteration terminates when the relative linear residual is
/// smaller than eta*| F(x_c) |. eta is determined by the modified
/// Eisenstat-Walker formula if etamax > 0. If etamax < 0, then eta = |etamax|
/// for the entire iteration.
/// \see updateForcingTerm()
///
void setEtaMax(const size_t newEtaMax)
{
this->etaMax = newEtaMax;
}
///
/// \brief Get EtaMax. Returns current etaMax value.
///
size_t getEtaMax() const
{
return this->etaMax;
}
///
/// \brief Set MaxIterations. The maximum number of nonlinear iterations.
///
void setMaxIterations(const size_t newMaxIterations)
{
this->maxIterations = newMaxIterations;
}
///
/// \brief Get MaxIterations. Returns current maximum nonlinear iterations.
///
size_t getMaxIterations() const
{
return this->maxIterations;
}
///
/// \brief Set the useArmijo flag. If useArmijo is true a line search is performed
/// using the Armijo-Goldstein condition.
///
void setUseArmijo(const bool value)
{
this->useArmijo = value;
}
///
/// \brief Get MaxIterations. Returns current maximum nonlinear iterations.
///
bool getUseArmijo() const
{
return this->useArmijo;
}
private:
std::shared_ptr<LinearSolverType> linearSolver; ///> Linear solver to use. Default: Conjugate gradient.
JacobianType jacobian; ///> Jacobian matrix function.
double forcingTerm; ///> Method's forcing term (Default: 0.9).
double absoluteTolerance; ///> Tolerance for the method (Default: 1.0e-3).
double relativeTolerance; ///> Relative (to the rhs) tolerance (Default: 1.0e-6).
double gamma; ///> Internal parameter used to update the forcing term (Default: 0.9).
double etaMax; ///> Maximum tolerance for the linear solver (Default: 0.9).
size_t maxIterations; ///> Maximum number of nonlinear iterations (Default: 50).
bool useArmijo; ///> True if Armijo liner search is desired (Default: true).
};
#endif // INEXACT_NEWTON_H
......@@ -81,7 +81,7 @@ double NonLinearSolver::armijo(const core::Vectord &dx, core::Vectord &x, const
if(i == this->armijoMax)
{
// TODO: Add to logger
std::cout << "Maximum number of Armijo iterations reached." << std::endl;
// std::cout << "Maximum number of Armijo iterations reached." << std::endl;
}
......
......@@ -62,23 +62,11 @@ public:
///
inline const core::Vectord &eval(const core::Vectord &x)
{
this->f = this->F(x);
return this->f;
}
///
/// \brief Get the value of the function F
///
/// \return Function value.
///
inline core::Vectord &getFunctionValue()
{
return this->f;
return this->F(x);
}
protected:
FunctionType F; ///> Function associated with the system of equation to solve.
core::Vectord f; ///> Storage for function values
};
///
......@@ -190,9 +178,20 @@ public:
return this->A.template triangularView<Eigen::StrictlyUpper>();
}
///
/// \brief Get the value of the function F
///
/// \return Function value.
///
inline core::Vectord &getFunctionValue()
{
return this->f;
}
private:
const MatrixType &A;
const core::Vectord &rhs;
core::Vectord f; ///> Storage for matrix vector product
};
#endif // SYSTEM_OF_EQUATIONS
......@@ -32,7 +32,7 @@
#include "Solvers/BackwardSOR.h"
#include "Solvers/DirectLinearSolver.h"
#include "Solvers/ConjugateGradient.h"
#include "Solvers/InexactNewton.h"
#include "Solvers/NewtonMethod.h"
#include "Core/Matrix.h"
#include "Testing/ReadSparseMatrix.h"
......@@ -190,7 +190,7 @@ go_bandit([]()
describe("Newton non-linear solver",[&]()
{
auto NewtonSolver = std::make_shared<InexactNewton>();
auto NewtonSolver = std::make_shared<NewtonMethod>();
it("constructs ", [&]()
{
AssertThat(NewtonSolver != nullptr, IsTrue());
......
......@@ -22,6 +22,7 @@
//---------------------------------------------------------------------------
#include "BackwarEuler.h"
#include "Solvers/DirectLinearSolver.h"
BackwardEuler::BackwardEuler(OdeSystem *odeSystem): TimeIntegrator(odeSystem)
{}
......@@ -36,71 +37,31 @@ void BackwardEuler::solve(const OdeSystemState &state,
return;
}
auto G = [this](const core::Vectord &) -> core::Vectord&
// This function updates the state.
auto updateIterate = [&](const core::Vectord &dv, core::Vectord &v)
{
return this->rhs;
v += dv;
newState.getPositions() = state.getPositions() + timeStep*v;
};
const auto &cState = state;
auto &cNewState = newState;
auto DG = [&,this](const core::Vectord &)
// Function to evaluate the nonlinear objective function.
auto G = [&,this](const core::Vectord &) -> const core::Vectord&
{
this->computeSystemMatrix(cState,cNewState,timeStep);
return this->systemMatrix;
this->system->computeImplicitSystemRHS(state,newState,timeStep);
return this->system->getRHS();
};
auto updateIterate = [&](const core::Vectord &dv, core::Vectord &v)
// Jacobian of the objective function.
auto DG = [&,this](const core::Vectord &) -> const core::SparseMatrixd&
{
v += dv;
cNewState.getPositions() = cState.getPositions() + timeStep*v;
this->system->computeImplicitSystemLHS(state,newState,timeStep,false);
return this->system->getSystemMatrix();
};
auto NewtonSolver = std::make_shared<InexactNewton>();
newState = state;
this->computeSystemRHS(state,newState,timeStep);
NewtonSolver->setUpdateIterate(updateIterate);
NewtonSolver->setSystem(G);
NewtonSolver->setJacobian(DG);
NewtonSolver->solve(newState.getVelocities());
}
//---------------------------------------------------------------------------
void BackwardEuler::computeSystemMatrix(const OdeSystemState &state,
OdeSystemState &newState,
const double timeStep,
bool computeRHS)
{
auto &M = this->system->evalMass(newState);
auto &K = this->system->evalDFx(newState);
auto &C = this->system->evalDFv(newState);
this->systemMatrix = (1.0 / timeStep) * M;
this->systemMatrix += C;
this->systemMatrix += timeStep * K;