Merge remote-tracking branch 'atgeirr/fully-implicit' into fully-implicit

opm/autodiff/AutoDiffHelpers.hpp

@@ -22,6 +22,7 @@
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/core/grid.h>
+#include <opm/core/utility/ErrorMacros.hpp>
 
 
 // -------------------- class HelperOps --------------------
@@ -424,4 +425,94 @@ vertcat(const AutoDiff::ForwardBlock<double>& x,
 
 
 
+class Span
+{
+public:
+    explicit Span(const int num)
+    : num_(num),
+      stride_(1),
+      start_(0)
+    {
+    }
+    Span(const int num, const int stride, const int start)
+        : num_(num),
+          stride_(stride),
+          start_(start)
+    {
+    }
+    int operator[](const int i) const
+    {
+        ASSERT(i >= 0 && i < num_);
+        return start_ + i*stride_;
+    }
+    int size() const
+    {
+        return num_;
+    }
+
+
+    class SpanIterator
+    {
+    public:
+        SpanIterator(const Span* span, const int index)
+            : span_(span),
+              index_(index)
+        {
+        }
+        SpanIterator operator++()
+        {
+            ++index_;
+            return *this;
+        }
+        SpanIterator operator++(int)
+        {
+            SpanIterator before_increment(*this);
+            ++index_;
+            return before_increment;
+        }
+        bool operator<(const SpanIterator& rhs) const
+        {
+            ASSERT(span_ == rhs.span_);
+            return index_ < rhs.index_;
+        }
+        bool operator==(const SpanIterator& rhs) const
+        {
+            ASSERT(span_ == rhs.span_);
+            return index_ == rhs.index_;
+        }
+        bool operator!=(const SpanIterator& rhs) const
+        {
+            ASSERT(span_ == rhs.span_);
+            return index_ != rhs.index_;
+        }
+    private:
+        const Span* span_;
+        int index_;
+    };
+
+    typedef SpanIterator iterator;
+    typedef SpanIterator const_iterator;
+
+    SpanIterator begin() const
+    {
+        return SpanIterator(this, 0);
+    }
+
+    SpanIterator end() const
+    {
+        return SpanIterator(this, num_);
+    }
+
+    bool operator==(const Span& rhs)
+    {
+        return num_ == rhs.num_ && start_ == rhs.start_ && stride_ == rhs.stride_;
+    }
+
+private:
+    const int num_;
+    const int stride_;
+    const int start_;
+};
+
+
 #endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED

opm/autodiff/ImpesTPFAAD.cpp

@@ -363,11 +363,7 @@ namespace Opm {
             divcontrib_sum = divcontrib_sum - divcontrib/cell_b;
             cell_residual_ = cell_residual_ - (component_contrib/cell_b);
             const ADB well_rates = perf_to_well * (perf_flux*perf_b);
-            std::vector<int> well_flow_res_phase_idx(nw);
-            for (int w = 0; w < nw; ++w) {
-                well_flow_res_phase_idx[w] = w + phase*nw;
-            }
-            qs_ = qs_ +  superset(well_rates, well_flow_res_phase_idx, nw*np);
+            qs_ = qs_ +  superset(well_rates, Span(nw, 1, phase*nw), nw*np);
         }
         cell_residual_ = cell_residual_ + divcontrib_sum;
         // Handling BHP and SURFACE_RATE wells.
@@ -426,15 +422,11 @@ namespace Opm {
             THROW("ImpesTPFAAD::solve(): Linear solver convergence failure.");
         }
         const V p0 = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
-        const V dp = subset(dx, buildAllCells(nc));
+        const V dp = subset(dx, Span(nc));
         const V p = p0 - dp;
         std::copy(&p[0], &p[0] + nc, state.pressure().begin());
         const V bhp0 = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
-        std::vector<int> bhp_dofs(nw);
-        for (int w = 0; w < nw; ++w) {
-            bhp_dofs[w] = nc + w;
-        }
-        ASSERT(bhp_dofs.back() + 1 == total_residual_.size());
+        Span bhp_dofs(nw, 1, nc);
         const V dbhp = subset(dx, bhp_dofs);
         const V bhp = bhp0 - dbhp;
         std::copy(&bhp[0], &bhp[0] + nw, well_state.bhp().begin());