move invDX to MSWellHelpers.hpp.

CMakeLists_files.cmake

@@ -252,6 +252,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/StandardWell_impl.hpp
   opm/autodiff/MultisegmentWell.hpp
   opm/autodiff/MultisegmentWell_impl.hpp
+  opm/autodiff/MSWellHelpers.hpp
   opm/autodiff/BlackoilWellModel.hpp
   opm/autodiff/BlackoilWellModel_impl.hpp
   opm/autodiff/StandardWellsSolvent.hpp

opm/autodiff/MSWellHelpers.hpp

@@ -0,0 +1,72 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef OPM_MSWELLHELPERS_HEADER_INCLUDED
+#define OPM_MSWELLHELPERS_HEADER_INCLUDED
+
+#include <dune/istl/solvers.hh>
+
+namespace Opm {
+
+namespace mswellhelpers
+{
+
+    // obtain y = D^-1 * x
+    template<typename MatrixType, typename VectorType>
+    VectorType
+    invDX(const MatrixType& D, VectorType x)
+    {
+        // the function will change the value of x, so we should not use reference of x here.
+
+        // TODO: store some of the following information to avoid to call it again and again for
+        // efficiency improvement.
+        // Bassically, only the solve / apply step is different.
+
+        VectorType y(x.size());
+        y = 0.;
+
+        Dune::MatrixAdapter<MatrixType, VectorType, VectorType> linearOperator(D);
+
+        // Sequential incomplete LU decomposition as the preconditioner
+        Dune::SeqILU0<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
+
+        // Preconditioned BICGSTAB solver
+        Dune::BiCGSTABSolver<VectorType> linsolver(linearOperator,
+                                                   preconditioner,
+                                                   1.e-6, // desired residual reduction factor
+                                                   50, // maximum number of iterations
+                                                   0); // verbosity of the solver
+
+        // Object storing some statistics about the solving process
+        Dune::InverseOperatorResult statistics ;
+
+        // Solve
+        linsolver.apply(y, x, statistics );
+
+        return y;
+    }
+
+
+} // namespace mswellhelpers
+
+}
+
+#endif

opm/autodiff/MultisegmentWell.hpp

@@ -25,8 +25,6 @@
 
 #include <opm/autodiff/WellInterface.hpp>
 
-#include <dune/istl/solvers.hh>
-
 namespace Opm
 {
 
@@ -334,42 +332,8 @@ namespace Opm
         void updateWellStateFromPrimaryVariables(WellState& well_state) const;
 
         double scalingFactor(const int comp_idx) const;
-    };
 
-    // obtain y = D^-1 * x
-    template<typename MatrixType, typename VectorType>
-    VectorType
-    invDX(const MatrixType& D, VectorType x)
-    {
-        // the function will change the value of x, so we should not use reference of x here.
 
-        // TODO: store some of the following information to avoid to call it again and again for
-        // efficiency improvement.
-        // Bassically, only the solve / apply step is different.
-
-        VectorType y(x.size());
-        y = 0.;
-
-        Dune::MatrixAdapter<MatrixType, VectorType, VectorType> linearOperator(D);
-
-        // Sequential incomplete LU decomposition as the preconditioner
-        Dune::SeqILU0<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
-
-        // Preconditioned BICGSTAB solver
-        Dune::BiCGSTABSolver<VectorType> linsolver(linearOperator,
-                                                   preconditioner,
-                                                   1.e-6, // desired residual reduction factor
-                                                   50, // maximum number of iterations
-                                                   0); // verbosity of the solver
-
-        // Object storing some statistics about the solving process
-        Dune::InverseOperatorResult statistics ;
-
-        // Solve
-        linsolver.apply(y, x, statistics );
-
-        return y;
-    }
 
 }
 

opm/autodiff/MultisegmentWell_impl.hpp

@@ -19,7 +19,7 @@
 */
 
 
-
+#include <opm/autodiff/MSWellHelpers.hpp>
 
 namespace Opm
 {
@@ -574,7 +574,7 @@ namespace Opm
         duneB_.mv(x, Bx);
 
         // invDBx = duneD^-1 * Bx_
-        BVectorWell invDBx = invDX(duneD_, Bx);
+        BVectorWell invDBx = mswellhelpers::invDX(duneD_, Bx);
 
         // Ax = Ax - duneC_^T * invDBx
         duneC_.mmtv(invDBx,Ax);
@@ -590,7 +590,7 @@ namespace Opm
     apply(BVector& r) const
     {
         // invDrw_ = duneD^-1 * resWell_
-        BVectorWell invDrw = invDX(duneD_, resWell_);
+        BVectorWell invDrw = mswellhelpers::invDX(duneD_, resWell_);
         // r = r - duneC_^T * invDrw
         duneC_.mmtv(invDrw, r);
     }
@@ -713,7 +713,7 @@ namespace Opm
         // resWell = resWell - B * x
         duneB_.mmv(x, resWell);
         // xw = D^-1 * resWell
-        xw = invDX(duneD_, resWell);
+        xw = mswellhelpers::invDX(duneD_, resWell);
     }
 
 
@@ -728,7 +728,7 @@ namespace Opm
     {
         // We assemble the well equations, then we check the convergence,
         // which is why we do not put the assembleWellEq here.
-        const BVectorWell dx_well = invDX(duneD_, resWell_);
+        const BVectorWell dx_well = mswellhelpers::invDX(duneD_, resWell_);
 
         updateWellState(dx_well, param, well_state);
     }