[NL] basic ODESolver types

Documentation/mainpage.dox

@@ -10,4 +10,9 @@
  * - Mailing list: https://groups.google.com/forum/#!forum/ogs-users
  * - Jenkins build server: https://svn.ufz.de/hudson
  *
+ *
+ * \section internal_modules Internal Modules
+ *
+ * * \subpage ODESolver
+ *
  */

NumLib/ODESolver/NonlinearSystem.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_NONLINEARSYSTEM_H
+#define NUMLIB_NONLINEARSYSTEM_H
+
+#include "Types.h"
+
+
+namespace NumLib
+{
+
+//! \addtogroup ODESolver
+//! @{
+
+/*! A System of nonlinear equations.
+ *
+ * \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 NonlinearSystem;
+
+
+/*! A System of nonlinear equations to be solved with the Newton-Raphson method.
+ *
+ * The Newton-Raphson method will iterate the linearized equation
+ * \f$ \mathtt{Jac} \cdot (-\Delta x_i) = \mathtt{res} \f$.
+ *
+ * \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 NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton>
+{
+public:
+    //! Assembles the residual at the point \c x.
+    virtual void assembleResidualNewton(Vector const& x) = 0;
+
+    //! Assembles the Jacobian of the residual at the point \c x.
+    virtual void assembleJacobian(Vector const& x) = 0;
+
+    /*! Writes the residual at point \c x to \c res.
+     *
+     * \pre assembleResidualNewton() must have been called before
+     *      with the same argument \c x.
+     *
+     * \todo Remove argument \c x.
+     */
+    virtual void getResidual(Vector const& x, Vector& res) const = 0;
+
+    /*! Writes the Jacobian of the residual to \c Jac.
+     *
+     * \pre assembleJacobian() must have been called before.
+     */
+    virtual void getJacobian(Matrix& Jac) const = 0;
+
+    //! Apply known solutions to the linearized equation system
+    //! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
+    virtual void applyKnownComponentsNewton(
+            Matrix& Jac, Vector& res, Vector& minus_delta_x) = 0;
+
+    /*! Check whether this is actually a linear equation system.
+     *
+     * \remark
+     * Depending on its parameters an in general nonlinear equation system
+     * can be linear in special cases. With this method it is possible to
+     * detect that at runtime and thus save an assembly call.
+     */
+    virtual bool isLinear() const = 0;
+
+    virtual ~NonlinearSystem() = default;
+};
+
+
+/*! A System of nonlinear equations to be solved with the Picard fixpoint
+ *  iteration method.
+ *
+ * The Picard method will iterate the linearized equation
+ * \f$ \mathtt{A} \cdot x_i = \mathtt{rhs} \f$.
+ *
+ * \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 NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard>
+{
+public:
+    //! Assembles the linearized eqation at point \c x.
+    virtual void assembleMatricesPicard(Vector const& x) = 0;
+
+    //! Writes the linearized equation system matrix to \c A.
+    virtual void getA(Matrix& A) const = 0;
+
+    //! Writes the linearized equation system right-hand side to \c rhs.
+    virtual void getRhs(Vector& rhs) const = 0;
+
+    //! Apply known solutions to the linearized equation system
+    //! \f$ A \cdot x = \mathit{rhs} \f$.
+    virtual void applyKnownComponentsPicard(
+            Matrix& A, Vector& rhs, Vector& x) = 0;
+
+    /*! Check whether this is actually a linear equation system.
+     *
+     * \remark
+     * Depending on its parameters an in general nonlinear equation system
+     * can be linear in special cases. With this method it is possible to
+     * detect that at runtime and thus save an assembly call.
+     */
+    virtual bool isLinear() const = 0;
+
+    virtual ~NonlinearSystem() = default;
+};
+
+//! @}
+
+}
+
+#endif // NUMLIB_NONLINEARSYSTEM_H

NumLib/ODESolver/ODESolver.dox

+/*! \defgroup ODESolver ODE Solver Library
+
+Put detailed ODE Solver Lib documentation here.
+
+*/

NumLib/ODESolver/ODESystem.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_ODESYSTEM_H
+#define NUMLIB_ODESYSTEM_H
+
+#include "Types.h"
+
+// TODO move to other namespace
+namespace ProcessLib
+{
+template <typename IndexType>
+struct DirichletBc;
+}
+
+namespace NumLib
+{
+
+//! \addtogroup ODESolver
+//! @{
+
+/*! ODE system interface.
+ *
+ * This is the interface an ODE has to implement in order to be solved with this
+ * ODE solver library.
+ *
+ * \tparam Matrix the type of matrices occuring in the linearization of the ODE.
+ * \tparam Vector the type of the solution vector of the ODE.
+ * \tparam ODETag a tag indicating the type of ODE.
+ * \tparam NLTag  a tag indicating the method used for resolving nonlinearities.
+ */
+template<typename Matrix, typename Vector, ODESystemTag ODETag, NonlinearSolverTag NLTag>
+class ODESystem;
+
+
+/*! Interface for a first-order implicit quasi-linear ODE.
+ *
+ * \tparam Matrix the type of matrices occuring in the linearization of the ODE.
+ * \tparam Vector the type of the solution vector of the ODE.
+ *
+ * \see ODESystemTag::FirstOrderImplicitQuasilinear
+ */
+template<typename Matrix, typename Vector>
+class ODESystem<Matrix, Vector,
+                ODESystemTag::FirstOrderImplicitQuasilinear,
+                NonlinearSolverTag::Picard>
+{
+public:
+    //! A tag indicating the type of ODE.
+    static const ODESystemTag ODETag = ODESystemTag::FirstOrderImplicitQuasilinear;
+
+    /*! Check whether this is actually a linear equation system.
+     *
+     * \remark
+     * Depending on its parameters an in general nonlinear ODE
+     * can be linear in special cases. With this method it is possible to
+     * detect that at runtime and thus save an assembly call.
+     */
+    virtual bool isLinear() const = 0;
+
+    //! Get the number of equations.
+    virtual std::size_t getNumEquations() const = 0;
+
+    //! Assemble \c M, \c K and \c b at the state (\c t, \c x).
+    virtual void assemble(const double t, Vector const& x,
+                          Matrix& M, Matrix& K, Vector& b) = 0;
+
+    using Index = typename MatrixTraits<Matrix>::Index;
+
+    // TODO doc
+    virtual std::vector<ProcessLib::DirichletBc<Index> > const* getKnownComponents() const
+    {
+        return nullptr; // by default there are no known components
+    }
+
+    virtual ~ODESystem() = default;
+};
+
+
+/*! Interface for a first-order implicit quasi-linear ODE.
+ *
+ * ODEs using this interface also provide a Jacobian.
+ */
+template<typename Matrix, typename Vector>
+class ODESystem<Matrix, Vector,
+                ODESystemTag::FirstOrderImplicitQuasilinear,
+                NonlinearSolverTag::Newton>
+        : public ODESystem<Matrix, Vector,
+                           ODESystemTag::FirstOrderImplicitQuasilinear,
+                           NonlinearSolverTag::Picard>
+{
+public:
+    /*! Assemble \f$ \mathtt{Jac} := \partial r/\partial x \f$ at the state (\c t, \c x).
+     *
+     * For the meaning of the other parameters refer to the the introductory remarks on
+     * \ref concept_time_discretization "time discretization".
+     *
+     * TODO document how to assemble the Jacobian!
+     */
+    virtual void assembleJacobian(const double t, Vector const& x, Vector const& xdot,
+                                  const double dxdot_dx, Matrix const& M,
+                                  const double dx_dx, Matrix const& K,
+                                  Matrix& Jac) = 0;
+};
+
+//! @}
+
+}
+
+#endif // NUMLIB_ODESYSTEM_H

NumLib/ODESolver/Types.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_TYPES_H
+#define NUMLIB_TYPES_H
+
+namespace NumLib
+{
+
+//! \addtogroup ODESolver
+//! @{
+
+//! Tag used to specify which nonlinear solver will be used.
+enum class NonlinearSolverTag : bool {
+    Picard /*!< Picard fixpoint iteration scheme */,
+    Newton /*!< Newton-Raphson iteration scheme */
+};
+
+//! Tag used to specify the type of ODE.
+enum class ODESystemTag : char
+{
+    /*! First order implicit quasi-linear ODE
+     *
+     * This is an ODE of the form
+     * \f$ M(x,t)\cdot \dot x + K(x,t) \cdot x - b(x,t)
+     *  =: r(\dot x, x, t) \stackrel{!}{=} 0 \f$
+     */
+    FirstOrderImplicitQuasilinear,
+    NeumannBC // Sure, that's misuse of this enum, so sue me!
+};
+
+//! @}
+
+
+// TODO move MatrixTraits to NumericsConfig.h?
+template<typename Matrix>
+struct MatrixTraits
+/*
+{
+    using Index = int;
+}
+// */
+;
+
+} // namespace NumLib
+
+
+
+#ifdef OGS_USE_EIGEN
+
+#include<Eigen/Core>
+
+namespace NumLib
+{
+template<>
+struct MatrixTraits<Eigen::MatrixXd>
+{
+    using Index = Eigen::MatrixXd::Index;
+};
+}
+
+#endif
+
+
+#ifdef USE_PETSC
+
+namespace MathLib { class PETScMatrix; }
+
+namespace NumLib
+{
+template<>
+struct MatrixTraits<MathLib::PETScMatrix>
+{
+    using Index = MathLib::PETScMatrix::IndexType;
+};
+}
+
+#elif defined(OGS_USE_EIGEN)
+
+namespace MathLib { class EigenMatrix; }
+
+namespace NumLib
+{
+template<>
+struct MatrixTraits<MathLib::EigenMatrix>
+{
+    using Index = MathLib::EigenMatrix::IndexType;
+};
+}
+
+#endif
+
+
+#endif // NUMLIB_TYPES_H