Add non neighboring connections induced by well perforations to matrix sparsity pattern

This is only done upon request and uses the auxiliary module approach
provided by ewoms.

CMakeLists_files.cmake

@@ -315,6 +315,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/SimulatorFullyImplicitBlackoilOutput.hpp

opm/autodiff/BlackoilModelEbos.hpp

@@ -29,6 +29,7 @@
 
 #include <opm/autodiff/BlackoilModelParameters.hpp>
 #include <opm/autodiff/BlackoilWellModel.hpp>
+#include <opm/autodiff/WellConnectionAuxiliaryModule.hpp>
 #include <opm/autodiff/BlackoilDetails.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 
@@ -353,6 +354,16 @@ namespace Opm {
         SimulatorReport assemble(const SimulatorTimerInterface& timer,
                                  const int iterationIdx)
         {
+            if ( param_.matrix_add_well_contributions_ )
+            {
+                // This might be dangerous?!
+                ebosSimulator_.model().clearAuxiliaryModules();
+                const auto* wells = wellModel().wellsPointer();
+                auto auxMod = std::make_shared<WellConnectionAuxiliaryModule<TypeTag> >(wells);
+                ebosSimulator_.model().addAuxiliaryModule(auxMod);
+            }
+
+
             // -------- Mass balance equations --------
             ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx);
             ebosSimulator_.problem().beginIteration();
@@ -1067,6 +1078,7 @@ namespace Opm {
 
     private:
 
+
         double dpMaxRel() const { return param_.dp_max_rel_; }
         double dsMax() const { return param_.ds_max_; }
         double drMaxRel() const { return param_.dr_max_rel_; }

opm/autodiff/BlackoilModelParameters.cpp

@@ -65,6 +65,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_);
     }
 
 
@@ -93,6 +94,7 @@ namespace Opm
         update_equations_scaling_ = false;
         use_update_stabilization_ = true;
         use_multisegment_well_ = false;
+        matrix_add_well_contributions_=false;
     }
 
 

opm/autodiff/BlackoilModelParameters.hpp

@@ -91,6 +91,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/WellConnectionAuxiliaryModule.hpp

@@ -0,0 +1,91 @@
+/*
+  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 <opm/core/wells.h>
+
+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 Wells* wells)
+        : wells_(wells)
+    {
+    }
+
+    unsigned numDofs() const
+    {
+        // No extra dofs are inserted for wells.
+        return 0;
+    }
+
+    void addNeighbors(std::vector<NeighborSet>& neighbors) const
+    {
+        const int nw = wells().number_of_wells;
+
+        for (int w = 0; w < nw; ++w)
+        {
+            const int nperfs = wells().well_connpos[w+1];
+            for (int perf = wells().well_connpos[w] ; perf < nperfs; ++perf) {
+                const auto cell1_idx = wells().well_cells[perf];
+                for(int perf1 = perf; perf1 < nperfs; ++perf1)
+                {
+                    const auto cell2_idx = wells().well_cells[perf1];
+                    neighbors[cell1_idx].insert(cell2_idx);
+                    neighbors[cell2_idx].insert(cell1_idx);
+                }
+            }
+        }
+    }
+
+    void applyInitial()
+    {}
+
+    void linearize(JacobianMatrix& matrix, GlobalEqVector& residual)
+    {
+        // Linearization is done in StandardDenseWells
+    }
+
+    private:
+
+    const Wells& wells() const
+    {
+        return *wells_;
+    }
+
+    const Wells* wells_;
+};
+
+} // end namespace OPM
+#endif