Merge branch 'feature/linearsolver-multitype-matrix' into 'master'

Feature/linearsolver multitype matrix

Closes #435

See merge request !744

dumux/linear/linearsolveracceptsmultitypematrix.hh

@@ -29,54 +29,48 @@
 namespace Dumux {
 
 //! The default
-template<typename TypeTag, typename LinearSolver>
-struct LinearSolverAcceptsMultiTypeMatrixImpl
+template<typename LinearSolver>
+struct LinearSolverAcceptsMultiTypeMatrix
 { static constexpr bool value = true; };
 
-/*!
- * \ingroup Linear
- * \brief Trait checking if linear solvers accept Dune::MultiTypeBlockMatrix or we need to convert the matrix
- */
-template <typename TypeTag>
-using LinearSolverAcceptsMultiTypeMatrix = LinearSolverAcceptsMultiTypeMatrixImpl<TypeTag, typename GET_PROP_TYPE(TypeTag, LinearSolver)>;
 
 //! Solvers that don't accept multi-type matrices
 //! Those are all with ILU preconditioner that doesn't support the additional block level
 //! And the direct solvers that have BCRS Matrix hardcoded
 
-template<typename TypeTag>
-struct LinearSolverAcceptsMultiTypeMatrixImpl<TypeTag, ILUnBiCGSTABBackend>
+template<>
+struct LinearSolverAcceptsMultiTypeMatrix<ILUnBiCGSTABBackend>
 { static constexpr bool value = false; };
 
-template<typename TypeTag>
-struct LinearSolverAcceptsMultiTypeMatrixImpl<TypeTag, ILUnCGBackend>
+template<>
+struct LinearSolverAcceptsMultiTypeMatrix<ILUnCGBackend>
 { static constexpr bool value = false; };
 
-template<typename TypeTag>
-struct LinearSolverAcceptsMultiTypeMatrixImpl<TypeTag, ILU0BiCGSTABBackend>
+template<>
+struct LinearSolverAcceptsMultiTypeMatrix<ILU0BiCGSTABBackend>
 { static constexpr bool value = false; };
 
-template<typename TypeTag>
-struct LinearSolverAcceptsMultiTypeMatrixImpl<TypeTag, ILU0CGBackend>
+template<>
+struct LinearSolverAcceptsMultiTypeMatrix<ILU0CGBackend>
 { static constexpr bool value = false; };
 
-template<typename TypeTag>
-struct LinearSolverAcceptsMultiTypeMatrixImpl<TypeTag, ILU0RestartedGMResBackend>
+template<>
+struct LinearSolverAcceptsMultiTypeMatrix<ILU0RestartedGMResBackend>
 { static constexpr bool value = false; };
 
-template<typename TypeTag>
-struct LinearSolverAcceptsMultiTypeMatrixImpl<TypeTag, ILUnRestartedGMResBackend>
+template<>
+struct LinearSolverAcceptsMultiTypeMatrix<ILUnRestartedGMResBackend>
 { static constexpr bool value = false; };
 
 #if HAVE_SUPERLU
-template<typename TypeTag>
-struct LinearSolverAcceptsMultiTypeMatrixImpl<TypeTag, SuperLUBackend>
+template<>
+struct LinearSolverAcceptsMultiTypeMatrix<SuperLUBackend>
 { static constexpr bool value = false; };
 #endif // HAVE_SUPERLU
 
 #if HAVE_UMFPACK
-template<typename TypeTag>
-struct LinearSolverAcceptsMultiTypeMatrixImpl<TypeTag, UMFPackBackend>
+template<>
+struct LinearSolverAcceptsMultiTypeMatrix<UMFPackBackend>
 { static constexpr bool value = false; };
 #endif // HAVE_UMFPACK
 

dumux/nonlinear/staggerednewtoncontroller.hh

@@ -72,49 +72,54 @@ public:
                            SolutionVector& x,
                            SolutionVector& b)
     {
+        // check matrix sizes
+        assert(checkMatrix_(A) && "Sub blocks of MultiType matrix have wrong sizes!");
+
         try
         {
             if (this->numSteps_ == 0)
-                this->initialResidual_ = b.two_norm();
-
-            // check matrix sizes
-            using namespace Dune::Indices;
-            assert(A[_0][_0].N() == A[_0][_1].N());
-            assert(A[_1][_0].N() == A[_1][_1].N());
-
-            // create the bcrs matrix the IterativeSolver backend can handle
-            const auto M = MatrixConverter<JacobianMatrix>::multiTypeToBCRSMatrix(A);
-
-            // get the new matrix sizes
-            const std::size_t numRows = M.N();
-            assert(numRows == M.M());
-
-            // create the vector the IterativeSolver backend can handle
-            const auto bTmp = VectorConverter<SolutionVector>::multiTypeToBlockVector(b);
-            assert(bTmp.size() == numRows);
+            {
+                Scalar norm2 = b.two_norm2();
+                if (this->comm().size() > 1)
+                    norm2 = this->comm().sum(norm2);
 
-            // create a blockvector to which the linear solver writes the solution
-            using VectorBlock = typename Dune::FieldVector<Scalar, 1>;
-            using BlockVector = typename Dune::BlockVector<VectorBlock>;
-            BlockVector y(numRows);
+                using std::sqrt;
+                this->initialResidual_ = sqrt(norm2);
+            }
 
-            // solve
-            const bool converged = ls.template solve</*precondBlockLevel=*/2>(M, y, bTmp);
+            // solve by calling the appropriate implementation depending on whether the linear solver
+            // is capable of handling MultiType matrices or not
+            const bool converged = solveLinearSystem_(ls, A, x, b,
+                                                      std::integral_constant<bool, LinearSolverAcceptsMultiTypeMatrix<LinearSolver>::value>());
 
-            // copy back the result y into x
-            VectorConverter<SolutionVector>::retrieveValues(x, y);
+            // make sure all processes converged
+            int convergedRemote = converged;
+            if (this->comm().size() > 1)
+                convergedRemote = this->comm().min(converged);
 
-            if (!converged)
-                DUNE_THROW(NumericalProblem, "Linear solver did not converge");
+            if (!converged) {
+                DUNE_THROW(NumericalProblem,
+                           "Linear solver did not converge");
+            }
+            else if (!convergedRemote) {
+                DUNE_THROW(NumericalProblem,
+                           "Linear solver did not converge on a remote process");
+            }
         }
-        catch (const Dune::Exception &e)
-        {
-            Dumux::NumericalProblem p;
+        catch (const Dune::Exception &e) {
+            // make sure all processes converged
+            int converged = 0;
+            if (this->comm().size() > 1)
+                converged = this->comm().min(converged);
+
+            NumericalProblem p;
             p.message(e.what());
             throw p;
         }
     }
 
+
+
     /*!
     * \brief Update the current solution with a delta vector.
     *
@@ -199,6 +204,90 @@ public:
             this->shift_ = this->comm().max(this->shift_);
     }
 
+private:
+    /*!
+     * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
+     *
+     * Throws Dumux::NumericalProblem if the linear solver didn't
+     * converge.
+     *
+     * Specialization for linear solvers that can handle MultiType matrices.
+     *
+     */
+    template<class LinearSolver, class JacobianMatrix, class SolutionVector>
+    bool solveLinearSystem_(LinearSolver& ls,
+                            JacobianMatrix& A,
+                            SolutionVector& x,
+                            SolutionVector& b,
+                            std::true_type)
+    {
+        // TODO: automatically derive the precondBlockLevel
+        return ls.template solve</*precondBlockLevel=*/2>(A, x, b);
+    }
+
+    /*!
+     * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
+     *
+     * Throws Dumux::NumericalProblem if the linear solver didn't
+     * converge.
+     *
+     * Specialization for linear solvers that cannot handle MultiType matrices.
+     * We copy the matrix into a 1x1 block BCRS matrix before solving.
+     *
+     */
+    template<class LinearSolver, class JacobianMatrix, class SolutionVector>
+    bool solveLinearSystem_(LinearSolver& ls,
+                            JacobianMatrix& A,
+                            SolutionVector& x,
+                            SolutionVector& b,
+                            std::false_type)
+    {
+        // create the bcrs matrix the IterativeSolver backend can handle
+        const auto M = MatrixConverter<JacobianMatrix>::multiTypeToBCRSMatrix(A);
+
+        // get the new matrix sizes
+        const std::size_t numRows = M.N();
+        assert(numRows == M.M());
+
+        // create the vector the IterativeSolver backend can handle
+        const auto bTmp = VectorConverter<SolutionVector>::multiTypeToBlockVector(b);
+        assert(bTmp.size() == numRows);
+
+        // create a blockvector to which the linear solver writes the solution
+        using VectorBlock = typename Dune::FieldVector<Scalar, 1>;
+        using BlockVector = typename Dune::BlockVector<VectorBlock>;
+        BlockVector y(numRows);
+
+        // solve
+        const bool converged = ls.solve(M, y, bTmp);
+
+        // copy back the result y into x
+        if(converged)
+            VectorConverter<SolutionVector>::retrieveValues(x, y);
+
+        return converged;
+    }
+
+    //! helper method to assure the MultiType matrix's sub blocks have the correct sizes
+    template<class JacobianMatrix>
+    bool checkMatrix_(const JacobianMatrix& A)
+    {
+        bool matrixHasCorrectSize = true;
+        using namespace Dune::Hybrid;
+        using namespace Dune::Indices;
+        forEach(A, [&matrixHasCorrectSize](const auto& rowOfMultiTypeMatrix)
+        {
+            const auto numRowsLeftMostBlock = rowOfMultiTypeMatrix[_0].N();
+
+            forEach(rowOfMultiTypeMatrix, [&matrixHasCorrectSize, &numRowsLeftMostBlock](const auto& subBlock)
+            {
+                if (subBlock.N() != numRowsLeftMostBlock)
+                    matrixHasCorrectSize = false;
+            });
+        });
+        return matrixHasCorrectSize;
+    }
+
 };
 
 } // end namespace Dumux