Merge branch 'matrix-with-well-connections' into test-well-amg

CMakeLists_files.cmake

@@ -213,6 +213,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/SimulatorIncompTwophaseAd.hpp
   opm/autodiff/SimulatorSequentialBlackoil.hpp
   opm/autodiff/TransportSolverTwophaseAd.hpp
+  opm/autodiff/WellConnectionAuxiliaryModule.hpp
   opm/autodiff/WellDensitySegmented.hpp
   opm/autodiff/WellStateFullyImplicitBlackoil.hpp
   opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp

opm/autodiff/BlackoilModelEbos.hpp

@@ -29,6 +29,7 @@
 
 #include <opm/autodiff/BlackoilModelParameters.hpp>
 #include <opm/autodiff/StandardWellsDense.hpp>
+#include <opm/autodiff/WellConnectionAuxiliaryModule.hpp>
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
 #include <opm/autodiff/GridHelpers.hpp>
@@ -478,14 +479,17 @@ namespace Opm {
             if( isParallel() )
             {
                 typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, true > Operator;
-                Operator opA(ebosJac, const_cast< ThisType& > (*this), istlSolver().parallelInformation() );
+                Operator opA(ebosJac, const_cast< ThisType& > (*this),
+                             param_.matrix_add_well_contributions_,
+                             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) );
+                Operator opA(ebosJac, const_cast< ThisType& > (*this),
+                             param_.matrix_add_well_contributions_ );
                 istlSolver().solve( opA, x, ebosResid );
             }
 
@@ -531,8 +535,11 @@ 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, WellModel& wellMod,
-              : A_( A ), wellMod_( wellMod ), comm_()
+                                  bool matrix_add_well_contributions,
+                                  const boost::any& parallelInformation = boost::any() )
+              : A_( A ), wellMod_( wellMod ), comm_(),
+                matrix_add_well_contributions_(matrix_add_well_contributions)
           {
 #if HAVE_MPI
             if( parallelInformation.type() == typeid(ParallelISTLInformation) )
@@ -547,22 +554,29 @@ namespace Opm {
           virtual void apply( const X& x, Y& y ) const
           {
             A_.mv( x, y );
-            // add well model modification to y
+
-            wellMod_.applyWellModelAdd(x, y );
+            if ( ! matrix_add_well_contributions_ )
+            {
+              // add well model modification to y
+              wellMod_.applyWellModelAdd(x, y );
 
 #if HAVE_MPI
-            if( comm_ )
+              if( comm_ )
-              comm_->project( y );
+                comm_->project( y );
 #endif
+            }
           }
 
           // y += \alpha * A * x
           virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
           {
             A_.usmv(alpha,x,y);
-            // add scaled well model modification to y
-            wellMod_.applyWellModelScaleAdd( alpha, x, y );
 
+            if ( ! matrix_add_well_contributions_ )
+            {
+                // add scaled well model modification to y
+                wellMod_.applyWellModelScaleAdd( alpha, x, y );
+            }
 #if HAVE_MPI
             if( comm_ )
               comm_->project( y );
@@ -580,6 +594,7 @@ namespace Opm {
           const matrix_type& A_ ;
           WellModel& wellMod_;
           std::unique_ptr< communication_type > comm_;
+          bool matrix_add_well_contributions_;
         };
 
         /// Apply an update to the primary variables, chopped if appropriate.

opm/autodiff/BlackoilModelParameters.cpp

@@ -56,6 +56,7 @@ namespace Opm
         update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_);
         use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_);
         deck_file_name_ = param.template get<std::string>("deck_filename");
+        matrix_add_well_contributions_ = param.getDefault("matrix_add_well_contributions", matrix_add_well_contributions_);
     }
 
 
@@ -78,6 +79,7 @@ namespace Opm
         solve_welleq_initially_ = true;
         update_equations_scaling_ = false;
         use_update_stabilization_ = true;
+        matrix_add_well_contributions_=false;
     }
 
 

opm/autodiff/BlackoilModelParameters.hpp

@@ -68,6 +68,9 @@ namespace Opm
         // The file name of the deck
         std::string deck_file_name_;
 
+        // Whether to add influences of wells between cells to the matrix
+        bool matrix_add_well_contributions_;
+
         /// Construct from user parameters or defaults.
         explicit BlackoilModelParameters( const ParameterGroup& param );
 

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -217,6 +217,13 @@ public:
         std::vector<std::vector<double>> originalFluidInPlace;
         std::vector<double> originalFluidInPlaceTotals;
 
+        if ( model_param_.matrix_add_well_contributions_ )
+        {
+            ebosSimulator_.model().clearAuxiliaryModules();
+            auto auxMod = std::make_shared<WellConnectionAuxiliaryModule<TypeTag> >(eclState(), grid());
+            ebosSimulator_.model().addAuxiliaryModule(auxMod);
+        }
+
         // Main simulation loop.
         while (!timer.done()) {
             // Report timestep.

opm/autodiff/StandardWellsDense.hpp

@@ -58,6 +58,27 @@
 #include <opm/simulators/WellSwitchingLogger.hpp>
 
 namespace Opm {
+namespace Detail
+{
+    //! calculates ret = A^T * B
+    template< class K, int m, int n, int p >
+    static inline void multMatrixTransposed ( const Dune::FieldMatrix< K, m, n > &A,
+                                              const Dune::FieldMatrix< K, n, p > &B,
+                                              Dune::FieldMatrix< K, m, p > &ret )
+    {
+        typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type;
+
+        for( size_type i = 0; i < m; ++i )
+        {
+            for( size_type j = 0; j < p; ++j )
+            {
+                ret[ i ][ j ] = K( 0 );
+                for( size_type k = 0; k < n; ++k )
+                    ret[ i ][ j ] += A[ k ][ i ] * B[ k ][ j ];
+            }
+        }
+    }
+}
 
 enum WellVariablePositions {
     XvarWell = 0,
@@ -338,6 +359,12 @@ enum WellVariablePositions {
                                      const WellState& well_state,
                                      const int well_number) const;
 
+            /// \! brief Modifies matrix to include influences of the well perforations.
+            ///
+            /// \param mat The linear system with the assembled mass balance
+            ///            equations
+            void addWellContributions(Mat& mat) const;
+
             using WellMapType = typename WellState::WellMapType;
             using WellMapEntryType = typename WellState::mapentry_t;
 

opm/autodiff/StandardWellsDense_impl.hpp

@@ -116,6 +116,42 @@ namespace Opm {
     }
 
 
+    template<class TypeTag>
+    void
+    StandardWellsDense<TypeTag>::
+    addWellContributions(Mat& mat) const
+    {
+        // We need to change matrx A as follows
+        // A -= B^T D^-1 C
+        // D is diagonal
+        // B and C have nw rows, nc colums and nonzero
+        // at (i,j) only if well i has perforation at cell j.
+        const int nw = wells().number_of_wells;
+        for(int w = 0; w < nw; ++w)
+        {
+            for(auto colB = duneB_[w].begin(), endB =  duneB_[w].end();
+                colB != endB; ++colB) {
+                auto pi = colB.index();
+                auto& row = mat[pi];
+                assert(colB.index() ==  pi);
+                auto col = row.begin();
+
+                for(auto colC = duneC_[w].begin(), endC=duneC_[w].end();
+                    colC != endC; ++colC) {
+                    auto pj = colC.index();
+                    // Move to col index pj
+                    while(col.index()<pj) ++col;
+                    assert(col.index() == pj);
+
+                    typename Mat::block_type tmp, tmp1;
+                    Dune::FMatrixHelp::multMatrix(*invDuneD_[w].begin(),  (*colC),
+                                                  tmp);
+                    Detail::multMatrixTransposed((*colB), tmp, tmp1);
+                    (*col) -= tmp1;
+                }
+            }
+        }
+    }
 
     template<typename TypeTag>
     SimulatorReport
@@ -243,6 +279,11 @@ namespace Opm {
 
         // do the local inversion of D.
         localInvert( invDuneD_ );
+
+        if ( param_.matrix_add_well_contributions_ )
+        {
+            addWellContributions( ebosJac );
+        }
     }
 
 

opm/autodiff/WellConnectionAuxiliaryModule.hpp

@@ -0,0 +1,128 @@
+/*
+  Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
+  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_WELLCONNECTIONAUXILIARYMODULE_HEADER_INCLUDED
+#define OPM_WELLCONNECTIONAUXILIARYMODULE_HEADER_INCLUDED
+
+#include <ewoms/aux/baseauxiliarymodule.hh>
+
+#include <dune/grid/CpGrid.hpp>
+
+#include <opm/core/wells.h>
+
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
+namespace Opm
+{
+template<class TypeTag>
+class WellConnectionAuxiliaryModule
+    : public Ewoms::BaseAuxiliaryModule<TypeTag>
+{
+    typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+
+public:
+
+    using NeighborSet = typename
+        Ewoms::BaseAuxiliaryModule<TypeTag>::NeighborSet;
+
+    WellConnectionAuxiliaryModule(const EclipseState& eclipseState,
+                                  const Dune::CpGrid& grid)
+    {
+        // Create cartesian to compressed mapping
+        const auto& globalCell = grid.globalCell();
+        const auto& cartesianSize = grid.logicalCartesianSize();
+
+        auto size = cartesianSize[0]*cartesianSize[1]*cartesianSize[2];
+
+        std::vector<int> cartesianToCompressed(size, -1);
+        auto begin = globalCell.begin();
+
+        for(auto cell = begin, end= globalCell.end(); cell != end; ++cell)
+        {
+            cartesianToCompressed[ *cell ] = cell - begin;
+        }
+
+        std::vector<const Opm::Well*>  wells  = eclipseState.getSchedule().getWells();
+        int last_time_step = eclipseState.getSchedule().getTimeMap().size()-1;
+        wells_.reserve(wells.size());
+
+        // initialize the additional cell connections introduced by wells.
+        for (const auto well : wells) {
+            std::vector<int> compressed_well_perforations;
+            // All possible completions of the well
+            const auto& completionSet = well->getCompletions(last_time_step);
+            compressed_well_perforations.reserve(completionSet.size());
+
+            for (size_t c=0; c<completionSet.size(); c++) {
+                const auto& completion = completionSet.get(c);
+                int i = completion.getI();
+                int j = completion.getJ();
+                int k = completion.getK();
+                int cart_grid_idx = i + cartesianSize[0]*(j + cartesianSize[1]*k);
+                int compressed_idx = cartesianToCompressed[cart_grid_idx];
+
+                if ( compressed_idx >= 0 ) // Ignore completions in inactive/remote cells.
+                {
+                    compressed_well_perforations.push_back(compressed_idx);
+                }
+            }
+
+            if( ! compressed_well_perforations.empty() )
+            {
+                std::sort(compressed_well_perforations.begin(),
+                          compressed_well_perforations.end());
+
+                wells_.push_back(compressed_well_perforations);
+            }
+        }
+    }
+
+    unsigned numDofs() const
+    {
+        // No extra dofs are inserted for wells.
+        return 0;
+    }
+
+    void addNeighbors(std::vector<NeighborSet>& neighbors) const
+    {
+        for(const auto well_perforations : wells_)
+        {
+            for(const auto& perforation : well_perforations)
+                neighbors[perforation].insert(well_perforations.begin(),
+                                              well_perforations.end());
+        }
+    }
+
+    void applyInitial()
+    {}
+
+    void linearize(JacobianMatrix& matrix, GlobalEqVector& residual)
+    {
+        // Linearization is done in StandardDenseWells
+    }
+
+    private:
+
+    std::vector<std::vector<int> > wells_;
+};
+
+} // end namespace OPM
+#endif