From 5d8da526797ab0bdea70e75a7491fc8028b8f5ed Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Thu, 1 Feb 2018 09:15:48 +0100
Subject: [PATCH] Use decoupling strategy of Scheichl/Masson in the CPR
 preconditioner.

Often it is claimed by CPR-AMG evangelists that this will make the pressure
system more elliptic. That may be the case. But even more important it also
decouples the pressure from the saturations.

Without this approach the pressure correction influences the smoothing of the
full system too much and e.g. for Norne CPR is worse than simple block ILU0.
---
 opm/autodiff/BlackoilAmg.hpp | 92 ++++++++++++++++++++++++++++++++++--
 1 file changed, 89 insertions(+), 3 deletions(-)

diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp
index 751599637..11d9989fa 100644
--- a/opm/autodiff/BlackoilAmg.hpp
+++ b/opm/autodiff/BlackoilAmg.hpp
@@ -27,6 +27,7 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/preconditioners.hh>
 #include <dune/istl/schwarz.hh>
+#include <dune/istl/operators.hh>
 #include <dune/istl/scalarproducts.hh>
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
@@ -53,6 +54,84 @@ namespace Opm
 namespace Detail
 {
 
+template<class M, class X, class Y, class T>
+Dune::MatrixAdapter<M,X,Y> createOperator(const Dune::MatrixAdapter<M,X,Y>&, const M& matrix, const T&)
+{
+    return Dune::MatrixAdapter<M,X,Y>(matrix);
+}
+
+template<class M, class X, class Y, class T>
+Dune::OverlappingSchwarzOperator<M,X,Y,T> createOperator(const Dune::OverlappingSchwarzOperator<M,X,Y,T>&,
+                                                         const M& matrix, const T& comm)
+{
+    return Dune::OverlappingSchwarzOperator<M,X,Y,T>(matrix, comm);
+}
+
+//! \brief Applies diagonal scaling to the discretization Matrix (Scheichl, 2003)
+//!
+//! See section 3.2.3 of Scheichl, Masson: Decoupling and Block Preconditioning for
+//! Sedimentary Basin Simulations, 2003.
+//! \param op The operator that stems from the discretization.
+//! \param comm The communication objecte describing the data distribution.
+//! \param pressureIndex The index of the pressure in the matrix block
+//! \retun A pair of the scaled matrix and the associated operator-
+template<class Operator, class Communication>
+std::tuple<std::unique_ptr<typename Operator::matrix_type>, Operator>
+scaleMatrixQuasiImpes(const Operator& op, const Communication& comm,
+                      std::size_t pressureIndex)
+{
+    using Matrix = typename Operator::matrix_type;
+    using Block = typename Matrix::block_type;
+    std::unique_ptr<Matrix> matrix(new Matrix(op.getmat()));
+
+    for ( auto& row : *matrix )
+    {
+        for ( auto& block : row )
+        {
+            for ( std::size_t i = 0; i < Block::rows; i++ )
+            {
+                if ( i != pressureIndex )
+                {
+                    for(std::size_t j=0; j < Block::cols; j++)
+                    {
+                        block[pressureIndex][j] += block[i][j];
+                    }
+                }
+            }
+        }
+    }
+    return std::make_tuple(std::move(matrix), createOperator(op, *matrix, comm));
+}
+
+//! \brief Applies diagonal scaling to the discretization Matrix (Scheichl, 2003)
+//!
+//! See section 3.2.3 of Scheichl, Masson: Decoupling and Block Preconditioning for
+//! Sedimentary Basin Simulations, 2003.
+//! \param vector The vector to scale
+//! \param pressureIndex The index of the pressure in the matrix block
+template<class Vector>
+void scaleVectorQuasiImpes(Vector& vector, std::size_t pressureIndex)
+{
+    using Block = typename Vector::block_type;
+
+    for ( auto& block: vector)
+    {
+        for ( std::size_t i = 0; i < Block::dimension; i++ )
+        {
+            if ( i != pressureIndex )
+            {
+                block[pressureIndex] += block[i];
+            }
+        }
+    }
+}
+
+//! \brief TMP to create the scalar pendant to a real block matrix, vector, smoother, etc.
+//!
+//! \code
+//! using Scalar = ScalarType<BlockType>::value;
+//! \endcode
+//! will get the corresponding scalar type.
 template<typename NonScalarType>
 struct ScalarType
 {
@@ -656,6 +735,7 @@ class BlackoilAmg
 {
 public:
     using Operator = O;
+    using Matrix = typename Operator::matrix_type;
     using Criterion = C;
     using Communication = P;
     using Smoother = S;
@@ -698,10 +778,13 @@ public:
 
     BlackoilAmg(const Operator& fineOperator, const Criterion& criterion,
                 const SmootherArgs& smargs, const Communication& comm)
-        : smoother_(Detail::constructSmoother<Smoother>(fineOperator,smargs,comm)),
+        : scaledMatrixOperator_(Detail::scaleMatrixQuasiImpes(fineOperator, comm,
+                                                              COMPONENT_INDEX)),
+          smoother_(Detail::constructSmoother<Smoother>(std::get<1>(scaledMatrixOperator_),
+                                                        smargs, comm)),
           levelTransferPolicy_(criterion, comm),
           coarseSolverPolicy_(smargs, criterion),
-          twoLevelMethod_(fineOperator, smoother_,
+          twoLevelMethod_(std::get<1>(scaledMatrixOperator_), smoother_,
                           levelTransferPolicy_,
                           coarseSolverPolicy_, 0, 1)
     {}
@@ -720,9 +803,12 @@ public:
     void apply(typename TwoLevelMethod::FineDomainType& v,
                const typename TwoLevelMethod::FineRangeType& d)
     {
-        twoLevelMethod_.apply(v, d);
+        auto scaledD = d;
+        Detail::scaleVectorQuasiImpes(scaledD, COMPONENT_INDEX);
+        twoLevelMethod_.apply(v, scaledD);
     }
 private:
+    std::tuple<std::unique_ptr<Matrix>, Operator> scaledMatrixOperator_;
     std::shared_ptr<Smoother> smoother_;
     LevelTransferPolicy levelTransferPolicy_;
     CoarseSolverPolicy coarseSolverPolicy_;