[NL] added nonlinear solver classes

NumLib/ODESolver/NonlinearSolver-impl.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <logog/include/logog.hpp>
+
+// for debugging
+// #include <iostream>
+
+#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
+{
+    sys.assembleMatricesPicard(x);
+}
+
+template<typename Matrix, typename Vector>
+bool
+NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard>::
+solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard> &sys, Vector &x)
+{
+    namespace BLAS = MathLib::BLAS;
+
+    bool success = false;
+
+    for (unsigned iteration=1; iteration<_maxiter; ++iteration)
+    {
+        sys.assembleMatricesPicard(x);
+        sys.getA(_A);
+        sys.getRhs(_rhs);
+        // TODO _x_new might not have been initialized
+        sys.applyKnownComponentsPicard(_A, _rhs, _x_new);
+
+        // std::cout << "A:\n" << Eigen::MatrixXd(A) << "\n";
+        // std::cout << "rhs:\n" << rhs << "\n\n";
+
+        oneShotLinearSolve(_A, _rhs, _x_new);
+
+        BLAS::aypx(x, -1.0, _x_new); // x = _x_new - x
+        auto const error = norm(x);
+        // INFO("  picard iteration %u error: %e", iteration, error);
+
+        x = _x_new;
+
+        if (error < _tol) {
+            success = true;
+            break;
+        }
+
+        if (sys.isLinear()) {
+            // INFO("  picard linear system. not looping");
+            success = true;
+            break;
+        }
+    }
+
+    return success;
+}
+
+
+template<typename Matrix, typename Vector>
+void
+NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::
+assemble(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys,
+         Vector const& x) const
+{
+    sys.assembleResidualNewton(x);
+}
+
+template<typename Matrix, typename Vector>
+bool
+NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::
+solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys, Vector &x)
+{
+    namespace BLAS = MathLib::BLAS;
+
+    bool success = false;
+
+    for (unsigned iteration=1; iteration<_maxiter; ++iteration)
+    {
+        sys.assembleResidualNewton(x);
+        sys.getResidual(x, _res);
+
+        // std::cout << "  res:\n" << res << std::endl;
+
+        if (norm(_res) < _tol) {
+            success = true;
+            break;
+        }
+
+        sys.assembleJacobian(x);
+        sys.getJacobian(_J);
+        // TODO _minus_delta_x might not have been initialized
+        sys.applyKnownComponentsNewton(_J, _res, _minus_delta_x);
+
+        // std::cout << "  J:\n" << Eigen::MatrixXd(J) << std::endl;
+
+        oneShotLinearSolve(_J, _res, _minus_delta_x);
+
+        // auto const dx_norm = _minus_delta_x.norm();
+        // INFO("  newton iteration %u, norm of delta x: %e", iteration, dx_norm);
+
+        BLAS::axpy(x, -1.0, _minus_delta_x);
+
+        if (sys.isLinear()) {
+            // INFO("  newton linear system. not looping");
+            success = true;
+            break;
+        }
+    }
+
+    return success;
+}
+
+}

NumLib/ODESolver/NonlinearSolver.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef NUMLIB_NONLINEARSOLVER_H
+#define NUMLIB_NONLINEARSOLVER_H
+
+#include <logog/include/logog.hpp>
+
+#include "Types.h"
+#include "NonlinearSystem.h"
+
+
+namespace NumLib
+{
+
+//! \addtogroup ODESolver
+//! @{
+
+/*! Find a solution to a nonlinear equation.
+ *
+ * \tparam Matrix the type of matrices occuring in the linearization of the equation.
+ * \tparam Vector the type of the solution vector of the equation.
+ * \tparam NLTag  a tag indicating the method used for solving the equation.
+ */
+template<typename Matrix, typename Vector, NonlinearSolverTag NLTag>
+class NonlinearSolver;
+
+
+/*! Find a solution to a nonlinear equation using the Newton-Raphson method.
+ *
+ * \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 NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton> final
+{
+public:
+    //! Type of the nonlinear equation system to be solved.
+    using System = NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton>;
+
+    /*! Constructs a new instance.
+     *
+     * \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)
+        , _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;
+
+    /*! 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);
+
+private:
+    const double _tol;       //!< tolerance of the solver
+    const unsigned _maxiter; //!< maximum number of iterations
+
+    Vector _res;           //!< The residual.
+    Matrix _J;             //!< The Jacobian of the residual.
+    Vector _minus_delta_x; //!< The Newton-Raphson method solves the linearized equation
+                           //!< \f$ J \cdot (-\Delta x) = r \f$ repeatedly.
+};
+
+
+/*! Find a solution to a nonlinear equation using the Picard fixpoint iteration method.
+ *
+ * \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 NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard> final
+{
+public:
+    //! Type of the nonlinear equation system to be solved.
+    using System = NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard>;
+
+    /*! Constructs a new instance.
+     *
+     * \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)
+        , _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;
+
+    /*! 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);
+
+private:
+    const double _tol;       //!< tolerance of the solver
+    const unsigned _maxiter; //!< maximum number of iterations
+
+    Matrix _A;     //!< \c Matrix describing the linearized system.
+    Vector _rhs;   //!< \c Vector describing the linearized system.
+    Vector _x_new; //!< \c Vector to store solutions of \f$ A \cdot x = \mathit{rhs} \f$.
+};
+
+//! @}
+
+}
+
+#include "NonlinearSolver-impl.h"
+
+#endif // NUMLIB_NONLINEARSOLVER_H