Merge branch 'feature/newton-unit-test' into 'master'

Feature/newton unit test

See merge request !1836

dumux/nonlinear/newtonsolver.hh

@@ -205,20 +205,6 @@ public:
     void setMaxSteps(int maxSteps)
     { maxSteps_ = maxSteps; }
 
-    /*!
-     * \brief Run the Newton method to solve a non-linear system.
-     *        Does time step control when the Newton fails to converge
-     */
-    [[deprecated("Use attachConvergenceWriter(convWriter) and solve(x, *timeLoop) instead")]]
-    void solve(SolutionVector& uCurrentIter, TimeLoop& timeLoop,
-               std::shared_ptr<ConvergenceWriter> convWriter)
-    {
-        if (!convergenceWriter_)
-            attachConvergenceWriter(convWriter);
-
-        solve(uCurrentIter, timeLoop);
-    }
-
     /*!
      * \brief Run the Newton method to solve a non-linear system.
      *        Does time step control when the Newton fails to converge
@@ -684,20 +670,6 @@ public:
         return oldTimeStep*(1.0 + percent/1.2);
     }
 
-    /*!
-     * \brief Specifies if the Newton method ought to be chatty.
-     */
-    [[deprecated("Has been replaced by setVerbosity(int). Will be removed after 3.1 release!")]]
-    void setVerbose(bool val)
-    { verbosity_ = val; }
-
-    /*!
-     * \brief Returns true if the Newton method ought to be chatty.
-     */
-    [[deprecated("Has been replaced by int verbosity(). Will be removed after 3.1 release!")]]
-    bool verbose() const
-    { return verbosity_ ; }
-
     /*!
      * \brief Specifies the verbosity level
      */
@@ -992,8 +964,9 @@ private:
     void newtonUpdateShiftImpl_(const SolVec &uLastIter,
                                 const SolVec &deltaU)
     {
-        for (int i = 0; i < int(uLastIter.size()); ++i) {
-            typename SolVec::block_type uNewI = uLastIter[i];
+        for (int i = 0; i < int(uLastIter.size()); ++i)
+        {
+            auto uNewI = uLastIter[i];
             uNewI -= deltaU[i];
 
             Scalar shiftAtDof = relativeShiftAtDof_(uLastIter[i], uNewI);
@@ -1023,7 +996,6 @@ private:
                                    const SolutionVector &deltaU)
     {
         Scalar lambda = 1.0;
-        SolutionVector tmp(uLastIter);
 
         while (true)
         {
@@ -1082,7 +1054,7 @@ private:
         //! to this field vector type in Dune ISTL
         //! Could be avoided for vectors that already have the right type using SFINAE
         //! but it shouldn't impact performance too much
-        constexpr auto blockSize = JacobianMatrix::block_type::rows;
+        constexpr auto blockSize = std::decay_t<decltype(b[0])>::dimension;
         using BlockType = Dune::FieldVector<Scalar, blockSize>;
         Dune::BlockVector<BlockType> xTmp; xTmp.resize(b.size());
         Dune::BlockVector<BlockType> bTmp(xTmp);

test/nonlinear/CMakeLists.txt

 add_subdirectory(findscalarroot)
+add_subdirectory(newton)

test/nonlinear/newton/CMakeLists.txt

+dumux_add_test(SOURCES test_newton.cc
+               LABELS unit nonlinear)
+dumux_add_test(NAME test_newton_linesearch
+               TARGET test_newton
+               COMMAND test_newton
+               CMD_ARGS "-Newton.UseLineSearch" "true"
+               LABELS unit nonlinear)

test/nonlinear/newton/test_newton.cc

+#include <config.h>
+
+#include <iostream>
+#include <cmath>
+#include <iomanip>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/float_cmp.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+/*
+
+  This test currently solves a scalar non-linear equation using the
+  Dumux::NewtonSolver. The Mock classes expose which dependencies the
+  current implementation has on different other classes it interacts with.
+  Several dependencies seem unnecessary. In particular, the current
+  implementation is basically hard-coded to double indexable residual vectors
+  with several math operators. A good idea would seem to somehow delegate
+  this dependency to something like a linear algebra backend or at least
+  the assembler. The assembler requires a lot of interface functions which
+  are not always needed. The linear solver interdependency is much better (small).
+
+  This test is to ensure that the dependencies do not grow more in the future.
+
+ */
+
+namespace Dumux {
+
+class MockScalarVariables
+{
+public:
+    static constexpr int dimension = 1;
+    MockScalarVariables() : var_(0.0) {}
+    explicit MockScalarVariables(double v) : var_(v) {}
+    MockScalarVariables& operator-=(const MockScalarVariables& other) { var_ -= other.var_; return *this; }
+    double& operator[] (int i) { return var_; }
+    const double& operator[] (int i) const { return var_; }
+private:
+    double var_;
+};
+
+class MockScalarResidual
+{
+public:
+    MockScalarResidual() : res_(0.0) {}
+    explicit MockScalarResidual(double r) : res_(r) {}
+    MockScalarResidual& operator=(double r) { res_[0] = r; return *this; }
+    MockScalarResidual& operator*=(double a) { res_[0] *= a; return *this; }
+    MockScalarResidual& operator+=(const MockScalarResidual& other) { res_[0] += other.res_[0]; return *this; }
+    MockScalarResidual& operator-=(const MockScalarResidual& other) { res_[0] -= other.res_[0]; return *this; }
+    constexpr std::size_t size() const { return 1; }
+    MockScalarVariables& operator[] (int i) { return res_; }
+    const MockScalarVariables& operator[] (int i) const { return res_; }
+    double two_norm2() const { return res_[0]*res_[0]; }
+private:
+    MockScalarVariables res_;
+};
+
+class MockScalarAssembler
+{
+public:
+    using ResidualType = MockScalarResidual;
+    using Scalar = double;
+    using JacobianMatrix = double;
+
+    void setLinearSystem() {}
+
+    bool isStationaryProblem() { return false; }
+
+    ResidualType prevSol() { return ResidualType(0.0); }
+
+    void resetTimeStep(const ResidualType& sol) {}
+
+    void assembleResidual(const ResidualType& sol)
+    {
+        res_[0][0] = sol[0][0]*sol[0][0] - 5.0;
+    }
+
+    void assembleJacobianAndResidual (const ResidualType& sol)
+    {
+        assembleResidual(sol);
+        jac_ = 2.0*sol[0][0];
+    }
+
+    JacobianMatrix& jacobian() { return jac_; }
+
+    ResidualType& residual() { return res_; }
+
+    double residualNorm(const ResidualType& sol)
+    {
+        assembleResidual(sol);
+        return res_[0][0];
+    }
+
+    void updateGridVariables(const ResidualType& sol) {}
+
+private:
+    JacobianMatrix jac_;
+    ResidualType res_;
+};
+
+class MockScalarLinearSolver
+{
+public:
+    void setResidualReduction(double residualReduction) {}
+
+    template<class Vector>
+    bool solve (const double& A, Vector& x, const Vector& b) const
+    {
+        x[0][0] = b[0][0]/A;
+        return true;
+    }
+};
+
+} // end namespace Dumux
+
+int main(int argc, char* argv[]) try
+{
+    using namespace Dumux;
+
+    // maybe initialize MPI
+    Dune::MPIHelper::instance(argc, argv);
+
+    // initialize  parameters
+    // TODO this is necessary because there are some global default used in the Newton solver
+    // Do we really need them to be global defaults???
+    Parameters::init(argc, argv);
+
+    // use the Newton solver to find a solution to a scalar equation
+    using Assembler = MockScalarAssembler;
+    using LinearSolver = MockScalarLinearSolver;
+    using Solver = NewtonSolver<Assembler, LinearSolver, DefaultPartialReassembler>;
+
+    auto assembler = std::make_shared<Assembler>();
+    auto linearSolver = std::make_shared<LinearSolver>();
+    auto solver = std::make_shared<Solver>(assembler, linearSolver);
+
+    double initialGuess = 0.1;
+    MockScalarResidual x(initialGuess);
+
+    std::cout << "Solving: x^2 - 5 = 0" << std::endl;
+    solver->solve(x);
+    std::cout << "Solution: " << std::setprecision(15) << x[0][0]
+              << ", exact: " << std::sqrt(5.0)
+              << ", error: " << std::abs(x[0][0]-std::sqrt(5.0))/std::sqrt(5.0)*100 << "%" << std::endl;
+
+    if (Dune::FloatCmp::ne(x[0][0], std::sqrt(5.0), 1e-13))
+        DUNE_THROW(Dune::Exception, "Didn't find correct root: " << std::setprecision(15) << x[0][0] << ", exact: " << std::sqrt(5.0));
+
+    return 0;
+
+}
+catch (const Dune::Exception& e)
+{
+    std::cout << e << std::endl;
+    return 1;
+}
+catch (...)
+{
+    std::cout << "Unknown exception thrown!" << std::endl;
+    return 1;
+}