From 908b2bdf81ea1ce669827bb474015a5fe90958af Mon Sep 17 00:00:00 2001
From: Kai Bao <kai.bao@sintef.no>
Date: Wed, 7 Jun 2017 14:23:43 +0200
Subject: [PATCH] making wellMod_ in WellModelMatrixAdapter is actually a well
 model

instead of the BlackoilModelEbos, which reduces some confusion and also
removes a few unnecessary wrapping functions.
---
 opm/autodiff/BlackoilModelEbos.hpp       | 29 +++++++-----------------
 opm/autodiff/StandardWellsDense.hpp      |  4 ++--
 opm/autodiff/StandardWellsDense_impl.hpp |  4 ++--
 3 files changed, 12 insertions(+), 25 deletions(-)

diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp
index 797d67f81..6b7d3383c 100644
--- a/opm/autodiff/BlackoilModelEbos.hpp
+++ b/opm/autodiff/BlackoilModelEbos.hpp
@@ -100,7 +100,6 @@ namespace Opm {
     /// upwind weighting of mobilities.
     class BlackoilModelEbos
     {
-        typedef BlackoilModelEbos ThisType;
     public:
         // ---------  Types and enums  ---------
         typedef BlackoilState ReservoirState;
@@ -437,18 +436,6 @@ namespace Opm {
             return istlSolver().iterations();
         }
 
-        template <class X, class Y>
-        void applyWellModelAdd(const X& x, Y& y )
-        {
-            wellModel().apply(x, y);
-        }
-
-        template <class X, class Y>
-        void applyWellModelScaleAdd(const Scalar alpha, const X& x, Y& y )
-        {
-            wellModel().applyScaleAdd(alpha, x, y);
-        }
-
         /// Solve the Jacobian system Jx = r where J is the Jacobian and
         /// r is the residual.
         void solveJacobianSystem(BVector& x, BVector& xw) const
@@ -462,15 +449,15 @@ namespace Opm {
             // Solve system.
             if( isParallel() )
             {
-                typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, true > Operator;
-                Operator opA(ebosJac, const_cast< ThisType& > (*this), istlSolver().parallelInformation() );
+                typedef WellModelMatrixAdapter< Mat, BVector, BVector, StandardWellsDense<TypeTag>, true > Operator;
+                Operator opA(ebosJac, well_model_, istlSolver().parallelInformation() );
                 assert( opA.comm() );
                 istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
             }
             else
             {
-                typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, false > Operator;
-                Operator opA(ebosJac, const_cast< ThisType& > (*this) );
+                typedef WellModelMatrixAdapter< Mat, BVector, BVector, StandardWellsDense<TypeTag>, false > Operator;
+                Operator opA(ebosJac, well_model_);
                 istlSolver().solve( opA, x, ebosResid );
             }
 
@@ -516,7 +503,7 @@ namespace Opm {
           };
 
           //! constructor: just store a reference to a matrix
-          WellModelMatrixAdapter (const M& A, WellModel& wellMod, const boost::any& parallelInformation = boost::any() )
+          WellModelMatrixAdapter (const M& A, const WellModel& wellMod, const boost::any& parallelInformation = boost::any() )
               : A_( A ), wellMod_( wellMod ), comm_()
           {
 #if HAVE_MPI
@@ -533,7 +520,7 @@ namespace Opm {
           {
             A_.mv( x, y );
             // add well model modification to y
-            wellMod_.applyWellModelAdd(x, y );
+            wellMod_.apply(x, y );
 
 #if HAVE_MPI
             if( comm_ )
@@ -546,7 +533,7 @@ namespace Opm {
           {
             A_.usmv(alpha,x,y);
             // add scaled well model modification to y
-            wellMod_.applyWellModelScaleAdd( alpha, x, y );
+            wellMod_.applyScaleAdd( alpha, x, y );
 
 #if HAVE_MPI
             if( comm_ )
@@ -563,7 +550,7 @@ namespace Opm {
 
         protected:
           const matrix_type& A_ ;
-          WellModel& wellMod_;
+          const WellModel& wellMod_;
           std::unique_ptr< communication_type > comm_;
         };
 
diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp
index 8ba09e20e..5fa3b97c6 100644
--- a/opm/autodiff/StandardWellsDense.hpp
+++ b/opm/autodiff/StandardWellsDense.hpp
@@ -160,10 +160,10 @@ enum WellVariablePositions {
             void apply( BVector& r) const;
 
             // subtract B*inv(D)*C * x from A*x
-            void apply(const BVector& x, BVector& Ax);
+            void apply(const BVector& x, BVector& Ax) const;
 
             // apply well model with scaling of alpha
-            void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax);
+            void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const;
 
             // xw = inv(D)*(rw - C*x)
             void recoverVariable(const BVector& x, BVector& xw) const;
diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp
index 11781c98b..ef50b7999 100644
--- a/opm/autodiff/StandardWellsDense_impl.hpp
+++ b/opm/autodiff/StandardWellsDense_impl.hpp
@@ -397,7 +397,7 @@ namespace Opm {
     template<typename TypeTag>
     void
     StandardWellsDense<TypeTag>::
-    apply(const BVector& x, BVector& Ax)
+    apply(const BVector& x, BVector& Ax) const
     {
         if ( ! localWellsActive() ) {
             return;
@@ -420,7 +420,7 @@ namespace Opm {
     template<typename TypeTag>
     void
     StandardWellsDense<TypeTag>::
-    applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax)
+    applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const
     {
         if ( ! localWellsActive() ) {
             return;