modification to to well assembly form well side to avoid asking for wells on all cells

ebos/eclproblem.hh

@@ -1785,6 +1785,7 @@ public:
         OPM_TIMEBLOCK_LOCAL(eclProblemSource);
         rate = 0.0;
 
+        // Add well contribution to source here.
         wellModel_.computeTotalRatesForDof(rate, globalDofIdx);
 
         // convert the source term from the total mass rate of the
@@ -1796,6 +1797,14 @@ public:
             assert(isfinite(rate[eqIdx]));
         }
 
+        // Add non-well sources.
+        addToSourceDense(rate, globalDofIdx, timeIdx);
+    }
+
+    void addToSourceDense(RateVector& rate,
+                          unsigned globalDofIdx,
+                          unsigned timeIdx) const
+    {
         if (enableAquifers_)
             aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
 

opm/simulators/wells/BlackoilWellModel.hpp

@@ -271,6 +271,10 @@ namespace Opm {
 
             void addWellContributions(SparseMatrixAdapter& jacobian) const;
 
+            // add source from wells to the reservoir matrix
+            void addReservoirSourceTerms(GlobalEqVector& residual,
+                                         SparseMatrixAdapter& jacobian) const;
+
             // called at the beginning of a report step
             void beginReportStep(const int time_step);
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1248,6 +1248,35 @@ namespace Opm {
         }
     }
 
+    template <typename TypeTag>
+    void BlackoilWellModel<TypeTag>::
+    addReservoirSourceTerms(GlobalEqVector& residual, SparseMatrixAdapter& jacobian) const
+    {
+        // NB this loop may write multiple times to the same element
+        // if a cell is perforated by more than one well, so it should
+        // not be OpenMP-parallelized.
+        for (const auto& well : well_container_) {
+            if (!well->isOperableAndSolvable() && !well->wellIsStopped()) {
+                continue;
+            }
+            const auto& cells = well->cells();
+            const auto& rates = well->connectionRates();
+            for (unsigned perfIdx = 0; perfIdx < rates.size(); ++perfIdx) {
+                unsigned cellIdx = cells[perfIdx];
+                auto rate = rates[perfIdx];
+                // Scalar volume =  ebosSimulator_.problem().volume(cellIdx,0);
+                rate *= -1.0;
+                VectorBlockType res(0.0);
+		using  MatrixBlockType = Opm::MatrixBlock<Scalar, numEq, numEq >;
+                MatrixBlockType bMat(0.0);
+                ebosSimulator_.model().linearizer().setResAndJacobi(res, bMat, rate);
+                residual[cellIdx] += res;
+                jacobian.addToBlock(cellIdx, cellIdx, bMat);
+            }
+        }
+    }
+
+
     template<typename TypeTag>
     int
     BlackoilWellModel<TypeTag>::

opm/simulators/wells/WellInterface.hpp

@@ -293,6 +293,11 @@ public:
                            const GroupState& group_state,
                            DeferredLogger& deferred_logger);
 
+    const std::vector<RateVector>& connectionRates() const
+    {
+        return connectionRates_;
+    }
+
 protected:
 
     // simulation parameters