diff --git a/NumLib/ODESolver/NonlinearSolver-impl.h b/NumLib/ODESolver/NonlinearSolver-impl.h
index 95acaadd1a1dbd868095c00f9725fa3cec6e1cc9..303934b62cff7195310542910b4809f424584072 100644
--- a/NumLib/ODESolver/NonlinearSolver-impl.h
+++ b/NumLib/ODESolver/NonlinearSolver-impl.h
@@ -12,34 +12,31 @@
// for debugging
// #include <iostream>
+#include "BaseLib/ConfigTree.h"
#include "MathLib/LinAlg/BLAS.h"
#include "NonlinearSolver.h"
#include "MathLib/LinAlg/VectorNorms.h"
-// TODO: change
-#include "ODETypes.h" // for one shot linear solver
-
-
namespace NumLib
{
template<typename Matrix, typename Vector>
void
NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard>::
-assemble(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard> &sys,
- Vector const& x) const
+assemble(Vector const& x) const
{
- sys.assembleMatricesPicard(x);
+ _equation_system->assembleMatricesPicard(x);
}
template<typename Matrix, typename Vector>
bool
NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard>::
-solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard> &sys, Vector &x)
+solve(Vector &x)
{
namespace BLAS = MathLib::BLAS;
+ auto& sys = *_equation_system;
bool success = false;
@@ -57,11 +54,16 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard> &sys, Vector &
// std::cout << "A:\n" << Eigen::MatrixXd(A) << "\n";
// std::cout << "rhs:\n" << rhs << "\n\n";
- oneShotLinearSolve(_A, _rhs, _x_new);
+ if (!_linear_solver.solve(_A, _rhs, _x_new)) {
+ ERR("The linear solver failed.");
+ x = _x_new;
+ success = false;
+ break;
+ }
// x is used as delta_x in order to compute the error.
BLAS::aypx(x, -1.0, _x_new); // x = _x_new - x
- auto const error = norm(x);
+ auto const error = BLAS::norm2(x);
// INFO(" picard iteration %u error: %e", iteration, error);
// Update x s.t. in the next iteration we will compute the right delta x
@@ -86,18 +88,21 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard> &sys, Vector &
template<typename Matrix, typename Vector>
void
NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::
-assemble(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys,
- Vector const& x) const
+assemble(Vector const& x) const
{
- sys.assembleResidualNewton(x);
+ _equation_system->assembleResidualNewton(x);
+ // TODO if the equation system would be reset to nullptr after each
+ // assemble() or solve() call, the user would be forced to set the
+ // equation every time and could not forget it.
}
template<typename Matrix, typename Vector>
bool
NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::
-solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys, Vector &x)
+solve(Vector &x)
{
namespace BLAS = MathLib::BLAS;
+ auto& sys = *_equation_system;
bool success = false;
@@ -114,7 +119,7 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys, Vector &
// std::cout << " res:\n" << res << std::endl;
// TODO streamline that, make consistent with Picard.
- if (norm(_res) < _tol) {
+ if (BLAS::norm2(_res) < _tol) {
success = true;
break;
}
@@ -125,7 +130,12 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys, Vector &
// std::cout << " J:\n" << Eigen::MatrixXd(J) << std::endl;
- oneShotLinearSolve(_J, _res, _minus_delta_x);
+ if (!_linear_solver.solve(_J, _res, _minus_delta_x)) {
+ ERR("The linear solver failed.");
+ BLAS::axpy(x, -_alpha, _minus_delta_x);
+ success = false;
+ break;
+ }
// auto const dx_norm = _minus_delta_x.norm();
// INFO(" newton iteration %u, norm of delta x: %e", iteration, dx_norm);
@@ -142,4 +152,37 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys, Vector &
return success;
}
+
+template<typename Matrix, typename Vector>
+std::pair<
+ std::unique_ptr<NonlinearSolverBase<Matrix, Vector> >,
+ NonlinearSolverTag
+>
+createNonlinearSolver(MathLib::LinearSolver<Matrix, Vector>& linear_solver,
+ BaseLib::ConfigTree const& config)
+{
+ using AbstractNLS = NonlinearSolverBase<Matrix, Vector>;
+
+ auto const type = config.getConfParam<std::string>("type");
+ auto const tol = config.getConfParam<double>("tol");
+ auto const max_iter = config.getConfParam<unsigned>("max_iter");
+
+ if (type == "Picard")
+ {
+ auto const tag = NonlinearSolverTag::Picard;
+ using ConcreteNLS = NonlinearSolver<Matrix, Vector, tag>;
+ return std::make_pair(std::unique_ptr<AbstractNLS>(
+ new ConcreteNLS{linear_solver, tol, max_iter}), tag);
+ }
+ else if (type == "Newton")
+ {
+ auto const tag = NonlinearSolverTag::Newton;
+ using ConcreteNLS = NonlinearSolver<Matrix, Vector, tag>;
+ return std::make_pair(std::unique_ptr<AbstractNLS>(
+ new ConcreteNLS{linear_solver, tol, max_iter}), tag);
+ }
+ std::abort();
+}
+
+
}
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index 0077fc3cf8f8bc95339cc107749c4115c6dfa69f..cae9903d359479039e60d8cbbc568fc364853044 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -10,15 +10,58 @@
#ifndef NUMLIB_NONLINEARSOLVER_H
#define NUMLIB_NONLINEARSOLVER_H
+#include <memory>
+#include <utility>
#include <logog/include/logog.hpp>
+#include "MathLib/LinAlg/LinearSolver.h"
+
#include "Types.h"
#include "NonlinearSystem.h"
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+// TODO Document in the ODE solver lib, which matrices and vectors that are passed around
+// as method arguments are guaranteed to be of the right size (and zeroed out) and
+// which are not.
namespace NumLib
{
+/*! Common interface for nonlinear solvers.
+ *
+ * \tparam Matrix the type of matrices occuring in the linearization of the equation.
+ * \tparam Vector the type of the solution vector of the equation.
+ */
+template<typename Matrix, typename Vector>
+class NonlinearSolverBase
+{
+public:
+ /*! Only assemble the equation system.
+ *
+ * \note This method is needed to preload CrankNicolson time discretization scheme.
+ * It is not used for the general solver steps; in those only the solve() method
+ * is needed.
+ *
+ * \param x the state at which the equation system will be assembled.
+ */
+ virtual void assemble(Vector const& x) const = 0;
+
+ /*! Assemble and solve the equation system.
+ *
+ * \param x in: the initial guess, out: the solution.
+ *
+ * \retval true if the equation system could be solved
+ * \retval false otherwise
+ */
+ virtual bool solve(Vector& x) = 0;
+
+ virtual ~NonlinearSolverBase() = default;
+};
+
//! \addtogroup ODESolver
//! @{
@@ -39,44 +82,38 @@ class NonlinearSolver;
*/
template<typename Matrix, typename Vector>
class NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton> final
+ : public NonlinearSolverBase<Matrix, Vector>
{
public:
//! Type of the nonlinear equation system to be solved.
using System = NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton>;
+ using LinearSolver = MathLib::LinearSolver<Matrix, Vector>;
/*! Constructs a new instance.
*
+ * \param linear_solver the linear solver used by this nonlinear solver.
* \param tol the tolerance of the solver. \todo Be more specific about that!
* \param maxiter the maximum number of iterations used to solve the equation.
*/
explicit
- NonlinearSolver(double const tol, const unsigned maxiter)
- : _tol(tol)
+ NonlinearSolver(LinearSolver& linear_solver,
+ double const tol, const unsigned maxiter)
+ : _linear_solver(linear_solver)
+ , _tol(tol)
, _maxiter(maxiter)
{}
- /*! Only assemble the equation system.
- *
- * \note This method is needed to preload CrankNicolson time discretization scheme.
- * It is not used for the general solver steps; in those only the solve() method
- * is needed.
- *
- * \param sys the equation system to be assembled.
- * \param x the state at which the equation system will be assembled.
- */
- void assemble(System& sys, Vector const& x) const;
+ //! Set the nonlinear equation system that will be solved.
+ void setEquationSystem(System& eq) { _equation_system = &eq; }
- /*! Assemble and solve the equation system.
- *
- * \param sys the equation system to be solved.
- * \param x in: the initial guess, out: the solution.
- *
- * \retval true if the equation system could be solved
- * \retval false otherwise
- */
- bool solve(System& sys, Vector& x);
+ void assemble(Vector const& x) const override;
+
+ bool solve(Vector& x) override;
private:
+ LinearSolver& _linear_solver;
+ System* _equation_system = nullptr;
+
const double _tol; //!< tolerance of the solver
const unsigned _maxiter; //!< maximum number of iterations
@@ -95,44 +132,38 @@ private:
*/
template<typename Matrix, typename Vector>
class NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard> final
+ : public NonlinearSolverBase<Matrix, Vector>
{
public:
//! Type of the nonlinear equation system to be solved.
using System = NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard>;
+ using LinearSolver = MathLib::LinearSolver<Matrix, Vector>;
/*! Constructs a new instance.
*
+ * \param linear_solver the linear solver used by this nonlinear solver.
* \param tol the tolerance of the solver. \todo Be more specific about that!
* \param maxiter the maximum number of iterations used to solve the equation.
*/
explicit
- NonlinearSolver(double const tol, const unsigned maxiter)
- : _tol(tol)
+ NonlinearSolver(LinearSolver& linear_solver,
+ double const tol, const unsigned maxiter)
+ : _linear_solver(linear_solver)
+ , _tol(tol)
, _maxiter(maxiter)
{}
- /*! Only assemble the equation system.
- *
- * \note This method is needed to preload CrankNicolson time discretization scheme.
- * It is not used for the general solver steps; in those only the solve() method
- * is needed.
- *
- * \param sys the equation system to be assembled.
- * \param x the state at which the equation system will be assembled.
- */
- void assemble(System& sys, Vector const& x) const;
+ //! Set the nonlinear equation system that will be solved.
+ void setEquationSystem(System& eq) { _equation_system = &eq; }
- /*! Assemble and solve the equation system.
- *
- * \param sys the equation system to be solved.
- * \param x in: the initial guess, out: the solution.
- *
- * \retval true if the equation system could be solved
- * \retval false otherwise
- */
- bool solve(System& sys, Vector& x);
+ void assemble(Vector const& x) const override;
+
+ bool solve(Vector& x) override;
private:
+ LinearSolver& _linear_solver;
+ System* _equation_system = nullptr;
+
const double _tol; //!< tolerance of the solver
const unsigned _maxiter; //!< maximum number of iterations
@@ -141,9 +172,27 @@ private:
Vector _x_new; //!< \c Vector to store solutions of \f$ A \cdot x = \mathit{rhs} \f$.
};
+/*! Creates a new nonlinear solver from the given configuration.
+ *
+ * \param linear_solver the linear solver that will be used by the nonlinear
+ * solver
+ * \param config configuration settings
+ *
+ * \return a pair <tt>(nl_slv, tag)</tt> where \c nl_slv is the generated
+ * nonlinear solver instance and the \c tag indicates if it uses the
+ * Picard or Newton-Raphson method
+ */
+template<typename Matrix, typename Vector>
+std::pair<
+ std::unique_ptr<NonlinearSolverBase<Matrix, Vector> >,
+ NonlinearSolverTag
+>
+createNonlinearSolver(MathLib::LinearSolver<Matrix, Vector>& linear_solver,
+ BaseLib::ConfigTree const& config);
+
//! @}
-}
+} // namespace NumLib
#include "NonlinearSolver-impl.h"