Merge pull request #670 from GitPaean/well_refactoring_multi_segment

Well refactoring multi segment

CMakeLists_files.cmake

@@ -47,6 +47,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/autodiff/VFPProdProperties.cpp
   opm/autodiff/VFPInjProperties.cpp
   opm/autodiff/WellMultiSegment.cpp
+	opm/autodiff/MultisegmentWells.cpp
   opm/autodiff/BlackoilSolventState.cpp
   opm/autodiff/ThreadHandle.hpp
   opm/polymer/PolymerState.cpp
@@ -83,6 +84,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_vfpproperties.cpp
   tests/test_singlecellsolves.cpp
   tests/test_solventprops_ad.cpp
+	tests/test_multisegmentwells.cpp
 	# tests/test_thresholdpressure.cpp
   )
 
@@ -90,6 +92,7 @@ list (APPEND TEST_DATA_FILES
   tests/fluid.data
   tests/VFPPROD1
   tests/VFPPROD2
+	tests/msw.data
   )
 
 
@@ -193,9 +196,13 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/VFPInjProperties.hpp
   opm/autodiff/WellStateMultiSegment.hpp
   opm/autodiff/WellMultiSegment.hpp
+  opm/autodiff/MultisegmentWells.hpp
+  opm/autodiff/MultisegmentWells_impl.hpp
   opm/autodiff/WellHelpers.hpp
   opm/autodiff/StandardWells.hpp
+  opm/autodiff/StandardWells_impl.hpp
   opm/autodiff/StandardWellsSolvent.hpp
+  opm/autodiff/StandardWellsSolvent_impl.hpp
   opm/polymer/CompressibleTpfaPolymer.hpp
   opm/polymer/GravityColumnSolverPolymer.hpp
   opm/polymer/GravityColumnSolverPolymer_impl.hpp

opm/autodiff/BlackoilMultiSegmentModel.hpp

@@ -27,6 +27,8 @@
 #include <opm/autodiff/WellMultiSegment.hpp>
 #include <opm/autodiff/StandardWells.hpp>
 
+#include <opm/autodiff/MultisegmentWells.hpp>
+
 namespace Opm {
 
     struct BlackoilMultiSegmentSolutionState : public DefaultBlackoilSolutionState
@@ -136,76 +138,7 @@ namespace Opm {
         using Base::param_;
         using Base::linsolver_;
 
-        // Pressure correction due to the different depth of the perforation
+        MultisegmentWells ms_wells_;
-        // and the cell center of the grid block
-        // For the non-segmented wells, since the perforation are forced to be
-        // at the center of the grid cell, it should be ZERO.
-        // It only applies to the mutli-segmented wells.
-        V well_perforation_cell_pressure_diffs_;
-
-        // Pressure correction due to the depth differennce between segment depth and perforation depth.
-        ADB well_segment_perforation_pressure_diffs_;
-
-        // The depth difference between segment nodes and perforations
-        V well_segment_perforation_depth_diffs_;
-
-        // the average of the fluid densities in the grid block
-        // which is used to calculate the hydrostatic head correction due to the depth difference of the perforation
-        // and the cell center of the grid block
-        V well_perforation_cell_densities_;
-
-        // the density of the fluid mixture in the segments
-        // which is calculated in an implicit way
-        ADB well_segment_densities_;
-
-        // the hydrostatic pressure drop between segment nodes
-        // calculated with the above density of fluid mixtures
-        // for the top segment, they should always be zero for the moment.
-        ADB well_segment_pressures_delta_;
-
-        // the surface volume of components in the segments
-        // the initial value at the beginning of the time step
-        std::vector<V>   segment_comp_surf_volume_initial_;
-
-        // the value within the current iteration.
-        std::vector<ADB> segment_comp_surf_volume_current_;
-
-        // the mass flow rate in the segments
-        ADB segment_mass_flow_rates_;
-
-        // the viscosity of the fluid mixture in the segments
-        // TODO: it is only used to calculate the Reynolds number as we know
-        //       maybe it is not better just to store the Reynolds number?
-        ADB segment_viscosities_;
-
-        const std::vector<WellMultiSegmentConstPtr> wells_multisegment_;
-
-        std::vector<int> top_well_segments_;
-
-        // segment volume by dt (time step)
-        // to handle the volume effects of the segment
-        V segvdt_;
-
-        // Well operations and data needed.
-        struct MultiSegmentWellOps {
-            explicit MultiSegmentWellOps(const std::vector<WellMultiSegmentConstPtr>& wells_ms);
-            Eigen::SparseMatrix<double> w2p;              // well -> perf (scatter)
-            Eigen::SparseMatrix<double> p2w;              // perf -> well (gather)
-            Eigen::SparseMatrix<double> w2s;              // well -> segment (scatter)
-            Eigen::SparseMatrix<double> s2w;              // segment -> well (gather)
-            Eigen::SparseMatrix<double> s2p;              // segment -> perf (scatter)
-            Eigen::SparseMatrix<double> p2s;              // perf -> segment (gather)
-            Eigen::SparseMatrix<double> s2s_inlets;       // segment -> its inlet segments
-            Eigen::SparseMatrix<double> s2s_outlet;       // segment -> its outlet segment
-            Eigen::SparseMatrix<double> topseg2w;         // top segment -> well
-            AutoDiffMatrix eliminate_topseg;              // change the top segment related to be zero
-            std::vector<int> well_cells;                  // the set of perforated cells
-            V conn_trans_factors;                         // connection transmissibility factors
-            bool has_multisegment_wells;                  // flag indicating whether there is any muli-segment well
-        };
-
-        MultiSegmentWellOps wops_ms_;
-
 
         using Base::stdWells;
         using Base::wells;
@@ -227,11 +160,15 @@ namespace Opm {
         using Base::asImpl;
         using Base::variableReservoirStateInitials;
 
-        const std::vector<WellMultiSegmentConstPtr>& wellsMultiSegment() const { return wells_multisegment_; }
+        // TODO: fixing the confusing naming
+        const MultisegmentWells& msWells() const { return ms_wells_; }
+        MultisegmentWells& msWells() { return ms_wells_; }
 
-        void updateWellControls(WellState& xw) const;
+        const std::vector<WellMultiSegmentConstPtr>& wellsMultiSegment() const { return msWells().wells(); }
 
+        const MultisegmentWells::MultisegmentWellOps& msWellOps() const { return msWells().wellOps(); }
 
+        // TODO: kept for now. to be removed soon.
         void updateWellState(const V& dwells,
                              WellState& well_state);
 
@@ -259,27 +196,11 @@ namespace Opm {
                     SolutionState& state,
                     WellState& well_state);
 
-        void
-        computeWellFlux(const SolutionState& state,
-                        const std::vector<ADB>& mob_perfcells,
-                        const std::vector<ADB>& b_perfcells,
-                        V& aliveWells,
-                        std::vector<ADB>& cq_s) const;
-
         void
         updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
                                          const SolutionState& state,
                                          WellState& xw) const;
 
-        void
-        addWellFluxEq(const std::vector<ADB>& cq_s,
-                      const SolutionState& state);
-
-        void
-        addWellControlEq(const SolutionState& state,
-                         const WellState& xw,
-                         const V& aliveWells);
-
         int numWellVars() const;
 
         void
@@ -290,14 +211,6 @@ namespace Opm {
                                       std::vector<ADB>& vars,
                                       SolutionState& state) const;
 
-        // Calculate the density of the mixture in the segments
-        // And the surface volume of the components in the segments by dt
-        void
-        computeSegmentFluidProperties(const SolutionState& state);
-
-        void
-        computeSegmentPressuresDelta(const SolutionState& state);
-
 
     };
 

opm/autodiff/MultisegmentWells.cpp

@@ -0,0 +1,234 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 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/>.
+*/
+
+#include <config.h>
+
+#include <opm/autodiff/MultisegmentWells.hpp>
+
+namespace Opm {
+
+
+    MultisegmentWells::
+    MultisegmentWellOps::MultisegmentWellOps(const std::vector<WellMultiSegmentConstPtr>& wells_ms)
+    {
+        // no multi-segment wells are involved by default.
+        has_multisegment_wells = false;
+
+        if (wells_ms.empty()) {
+            return;
+        }
+
+        // Count the total number of perforations and segments.
+        const int nw = wells_ms.size();
+        int total_nperf = 0;
+        int total_nseg = 0;
+        for (int w = 0; w < nw; ++w) {
+            total_nperf += wells_ms[w]->numberOfPerforations();
+            total_nseg += wells_ms[w]->numberOfSegments();
+            if (wells_ms[w]->isMultiSegmented()) {
+                has_multisegment_wells = true;
+            }
+        }
+
+        // Create well_cells and conn_trans_factors.
+        well_cells.reserve(total_nperf);
+        conn_trans_factors.resize(total_nperf);
+        int well_perf_start = 0;
+        for (int w = 0; w < nw; ++w) {
+            WellMultiSegmentConstPtr well = wells_ms[w];
+            well_cells.insert(well_cells.end(), well->wellCells().begin(), well->wellCells().end());
+            const std::vector<double>& perf_trans = well->wellIndex();
+            std::copy(perf_trans.begin(), perf_trans.end(), conn_trans_factors.data() + well_perf_start);
+            well_perf_start += well->numberOfPerforations();
+        }
+        assert(well_perf_start == total_nperf);
+        assert(int(well_cells.size()) == total_nperf);
+
+        // Create all the operator matrices,
+        // using the setFromTriplets() method.
+        s2s_inlets = Eigen::SparseMatrix<double>(total_nseg, total_nseg);
+        s2s_outlet = Eigen::SparseMatrix<double>(total_nseg, total_nseg);
+        s2w        = Eigen::SparseMatrix<double>(nw, total_nseg);
+        w2s        = Eigen::SparseMatrix<double>(total_nseg, nw);
+        topseg2w   = Eigen::SparseMatrix<double>(nw, total_nseg);
+        s2p = Eigen::SparseMatrix<double>(total_nperf, total_nseg);
+        p2s = Eigen::SparseMatrix<double>(total_nseg, total_nperf);
+        typedef Eigen::Triplet<double> Tri;
+        std::vector<Tri> s2s_inlets_vector;
+        std::vector<Tri> s2s_outlet_vector;
+        std::vector<Tri> s2w_vector;
+        std::vector<Tri> w2s_vector;
+        std::vector<Tri> topseg2w_vector;
+        std::vector<Tri> s2p_vector;
+        std::vector<Tri> p2s_vector;
+        Vector topseg_zero = Vector::Ones(total_nseg);
+        s2s_inlets_vector.reserve(total_nseg);
+        s2s_outlet_vector.reserve(total_nseg);
+        s2w_vector.reserve(total_nseg);
+        w2s_vector.reserve(total_nseg);
+        topseg2w_vector.reserve(nw);
+        s2p_vector.reserve(total_nperf);
+        p2s_vector.reserve(total_nperf);
+        int seg_start = 0;
+        int perf_start = 0;
+        for (int w = 0; w < nw; ++w) {
+            const int ns = wells_ms[w]->numberOfSegments();
+            const int np = wells_ms[w]->numberOfPerforations();
+            for (int seg = 0; seg < ns; ++seg) {
+                const int seg_ind = seg_start + seg;
+                w2s_vector.push_back(Tri(seg_ind, w, 1.0));
+                s2w_vector.push_back(Tri(w, seg_ind, 1.0));
+                if (seg == 0) {
+                    topseg2w_vector.push_back(Tri(w, seg_ind, 1.0));
+                    topseg_zero(seg_ind) = 0.0;
+                }
+                int seg_outlet = wells_ms[w]->outletSegment()[seg];
+                if (seg_outlet >= 0) {
+                    const int outlet_ind = seg_start + seg_outlet;
+                    s2s_inlets_vector.push_back(Tri(outlet_ind, seg_ind, 1.0));
+                    s2s_outlet_vector.push_back(Tri(seg_ind, outlet_ind, 1.0));
+                }
+
+                const auto& seg_perf = wells_ms[w]->segmentPerforations()[seg];
+                // the number of perforations related to this segment
+                const int npseg = seg_perf.size();
+                for (int perf = 0; perf < npseg; ++perf) {
+                    const int perf_ind = perf_start + seg_perf[perf];
+                    s2p_vector.push_back(Tri(perf_ind, seg_ind, 1.0));
+                    p2s_vector.push_back(Tri(seg_ind, perf_ind, 1.0));
+                }
+            }
+            seg_start += ns;
+            perf_start += np;
+        }
+
+        s2s_inlets.setFromTriplets(s2s_inlets_vector.begin(), s2s_inlets_vector.end());
+        s2s_outlet.setFromTriplets(s2s_outlet_vector.begin(), s2s_outlet_vector.end());
+        w2s.setFromTriplets(w2s_vector.begin(), w2s_vector.end());
+        s2w.setFromTriplets(s2w_vector.begin(), s2w_vector.end());
+        topseg2w.setFromTriplets(topseg2w_vector.begin(), topseg2w_vector.end());
+        s2p.setFromTriplets(s2p_vector.begin(), s2p_vector.end());
+        p2s.setFromTriplets(p2s_vector.begin(), p2s_vector.end());
+
+        w2p = Eigen::SparseMatrix<double>(total_nperf, nw);
+        p2w = Eigen::SparseMatrix<double>(nw, total_nperf);
+        w2p = s2p * w2s;
+        p2w = s2w * p2s;
+
+        eliminate_topseg = AutoDiffMatrix(topseg_zero.matrix().asDiagonal());
+    }
+
+
+
+
+
+    MultisegmentWells::
+    MultisegmentWells(const std::vector<WellMultiSegmentConstPtr>& wells_ms, const int np)
+      : wells_multisegment_(wells_ms)
+      , wops_ms_(wells_ms)
+      , num_phases_(np)
+      , well_segment_perforation_pressure_diffs_(ADB::null())
+      , well_segment_densities_(ADB::null())
+      , well_segment_pressures_delta_(ADB::null())
+      , segment_comp_surf_volume_initial_(num_phases_)
+      , segment_comp_surf_volume_current_(num_phases_, ADB::null())
+      , segment_mass_flow_rates_(ADB::null())
+      , segment_viscosities_(ADB::null())
+    {
+        // TODO: repeated with the wellOps's initialization, delete one soon.
+        // Count the total number of perforations and segments.
+        const int nw = wells_ms.size();
+        top_well_segments_.resize(nw);
+        int nperf_total = 0;
+        int nseg_total = 0;
+        for (int w = 0; w < nw; ++w) {
+            top_well_segments_[w] = nseg_total;
+            nperf_total += wells_ms[w]->numberOfPerforations();
+            nseg_total += wells_ms[w]->numberOfSegments();
+        }
+
+        nperf_total_ = nperf_total;
+        nseg_total_ = nseg_total;
+    }
+
+
+
+
+
+    const std::vector<WellMultiSegmentConstPtr>&
+    MultisegmentWells::wells() const
+    {
+        return wells_multisegment_;
+    }
+
+
+
+
+
+    const MultisegmentWells::MultisegmentWellOps&
+    MultisegmentWells::wellOps() const
+    {
+        return wops_ms_;
+    }
+
+
+
+
+
+    void
+    MultisegmentWells::
+    computeSegmentPressuresDelta(const double grav)
+    {
+        const int nw = wells().size();
+        const int nseg_total = nseg_total_;
+
+        if ( !wellOps().has_multisegment_wells ) {
+            well_segment_pressures_delta_ = ADB::constant(Vector::Zero(nseg_total));
+            well_segment_perforation_pressure_diffs_ = wellOps().s2p * well_segment_pressures_delta_;
+            return;
+        }
+
+        // calculate the depth difference of the segments
+        // TODO: we need to store the following values somewhere to avoid recomputation.
+        Vector segment_depth_delta = Vector::Zero(nseg_total);
+        int start_segment = 0;
+        for (int w = 0; w < nw; ++w) {
+            WellMultiSegmentConstPtr well = wells()[w];
+            const int nseg = well->numberOfSegments();
+            for (int s = 1; s < nseg; ++s) {
+                const int s_outlet = well->outletSegment()[s];
+                assert(s_outlet >= 0 && s_outlet < nseg);
+                segment_depth_delta[s + start_segment] = well->segmentDepth()[s_outlet] - well->segmentDepth()[s];
+            }
+            start_segment += nseg;
+        }
+        assert(start_segment == nseg_total);
+
+        const ADB grav_adb = ADB::constant(Vector::Constant(nseg_total, grav));
+        well_segment_pressures_delta_ = segment_depth_delta * grav_adb * well_segment_densities_;
+
+        ADB well_segment_perforation_densities = wellOps().s2p * well_segment_densities_;
+        well_segment_perforation_pressure_diffs_ = grav * well_segment_perforation_depth_diffs_ * well_segment_perforation_densities;
+    }
+
+} // end of namespace Opm
+
+
+

opm/autodiff/MultisegmentWells.hpp

@@ -0,0 +1,227 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 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_MULTISEGMENTWELLS_HEADER_INCLUDED
+#define OPM_MULTISEGMENTWELLS_HEADER_INCLUDED
+
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <Eigen/Eigen>
+#include <Eigen/Sparse>
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+
+#include <cassert>
+
+#include <opm/core/props/BlackoilPhases.hpp>
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/BlackoilModelEnums.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
+#include <opm/autodiff/WellHelpers.hpp>
+
+#include <opm/autodiff/WellMultiSegment.hpp>
+
+
+
+namespace Opm {
+
+
+    /// Class for handling the multi-segment well model
+    class MultisegmentWells {
+    public:
+            // ---------      Types      ---------
+            using ADB = AutoDiffBlock<double>;
+            using Vector = ADB::V;
+
+            // Well operations and data needed.
+            struct MultisegmentWellOps {
+                explicit MultisegmentWellOps(const std::vector<WellMultiSegmentConstPtr>& wells_ms);
+                Eigen::SparseMatrix<double> w2p;              // well -> perf (scatter)
+                Eigen::SparseMatrix<double> p2w;              // perf -> well (gather)
+                Eigen::SparseMatrix<double> w2s;              // well -> segment (scatter)
+                Eigen::SparseMatrix<double> s2w;              // segment -> well (gather)
+                Eigen::SparseMatrix<double> s2p;              // segment -> perf (scatter)
+                Eigen::SparseMatrix<double> p2s;              // perf -> segment (gather)
+                Eigen::SparseMatrix<double> s2s_inlets;       // segment -> its inlet segments
+                Eigen::SparseMatrix<double> s2s_outlet;       // segment -> its outlet segment
+                Eigen::SparseMatrix<double> topseg2w;         // top segment -> well
+                AutoDiffMatrix eliminate_topseg;              // change the top segment related to be zero
+                std::vector<int> well_cells;                  // the set of perforated cells
+                Vector conn_trans_factors;                         // connection transmissibility factors
+                bool has_multisegment_wells;                  // flag indicating whether there is any muli-segment well
+            };
+
+            // copied from BlackoilModelBase
+            // should put to somewhere better
+            using DataBlock =  Eigen::Array<double,
+                                            Eigen::Dynamic,
+                                            Eigen::Dynamic,
+                                            Eigen::RowMajor>;
+
+            // ---------  Public methods  ---------
+            // TODO: using a vector of WellMultiSegmentConstPtr for now
+            // TODO: it should use const Wells or something else later.
+            MultisegmentWells(const std::vector<WellMultiSegmentConstPtr>& wells_multisegment, const int np);
+
+            const std::vector<WellMultiSegmentConstPtr>& wells() const;
+            const MultisegmentWellOps& wellOps() const;
+
+            int numPhases() const { return num_phases_; };
+
+            int numSegment() const { return nseg_total_; };
+
+            int numPerf() const { return nperf_total_; };
+
+            Vector& wellPerforationCellPressureDiffs() { return well_perforation_cell_pressure_diffs_; };
+
+            Vector& wellSegmentPerforationDepthDiffs() { return well_segment_perforation_depth_diffs_; };
+
+            const Vector& wellPerforationCellDensities() const { return well_perforation_cell_densities_; };
+            Vector& wellPerforationCellDensities() { return well_perforation_cell_densities_; };
+
+            const std::vector<Vector>& segmentCompSurfVolumeInitial() const { return segment_comp_surf_volume_initial_; };
+            std::vector<Vector>& segmentCompSurfVolumeInitial() { return segment_comp_surf_volume_initial_; };
+
+            const std::vector<ADB>& segmentCompSurfVolumeCurrent() const { return segment_comp_surf_volume_current_; };
+
+            const std::vector<int>& topWellSegments() const { return top_well_segments_; };
+            std::vector<int>& topWellSegments() { return top_well_segments_; };
+
+            Vector& segVDt() { return segvdt_; };
+
+
+
+            template <class WellState>
+            void
+            updateWellState(const Vector& dwells,
+                            const double dpmaxrel,
+                            WellState& well_state) const;
+
+            // TODO: some arguments can be removed later
+            // TODO: compi will be required in the multisegment wells
+            template <class SolutionState>
+            void
+            computeWellFlux(const SolutionState& state,
+                            const Opm::PhaseUsage& pu,
+                            const std::vector<bool>& active,
+                            const Vector& well_perforation_pressure_diffs,
+                            const DataBlock& compi,
+                            const std::vector<ADB>& mob_perfcells,
+                            const std::vector<ADB>& b_perfcells,
+                            Vector& aliveWells,
+                            std::vector<ADB>& cq_s) const;
+
+
+            // Calculate the density of the mixture in the segments
+            // And the surface volume of the components in the segments by dt
+            template <class SolutionState>
+            void
+            computeSegmentFluidProperties(const SolutionState& state,
+                                          const std::vector<PhasePresence>& pc,
+                                          const std::vector<bool>& active,
+                                          const BlackoilPropsAdInterface& fluid);
+
+            void
+            computeSegmentPressuresDelta(const double grav);
+
+            template <class SolutionState>
+            void
+            addWellFluxEq(const std::vector<ADB>& cq_s,
+                          const SolutionState& state,
+                          LinearisedBlackoilResidual& residual);
+
+            template <class SolutionState, class WellState>
+            void
+            addWellControlEq(const SolutionState& state,
+                             const WellState& xw,
+                             const Vector& aliveWells,
+                             const std::vector<bool>& active,
+                             LinearisedBlackoilResidual& residual);
+
+            template <class WellState>
+            void
+            updateWellControls(const bool terminal_output,
+                               WellState& xw) const;
+
+    protected:
+        // TODO: probably a wells_active_ will be required here.
+        const std::vector<WellMultiSegmentConstPtr> wells_multisegment_;
+        const MultisegmentWellOps wops_ms_;
+        const int num_phases_;
+        int nseg_total_;
+        int nperf_total_;
+
+        // Pressure correction due to the different depth of the perforation
+        // and the cell center of the grid block
+        // For the non-segmented wells, since the perforation are forced to be
+        // at the center of the grid cell, it should be ZERO.
+        // It only applies to the mutli-segmented wells.
+        Vector well_perforation_cell_pressure_diffs_;
+
+        // Pressure correction due to the depth differennce between segment depth and perforation depth.
+        ADB well_segment_perforation_pressure_diffs_;
+
+        // The depth difference between segment nodes and perforations
+        Vector well_segment_perforation_depth_diffs_;
+
+        // the average of the fluid densities in the grid block
+        // which is used to calculate the hydrostatic head correction due to the depth difference of the perforation
+        // and the cell center of the grid block
+        Vector well_perforation_cell_densities_;
+
+        // the density of the fluid mixture in the segments
+        // which is calculated in an implicit way
+        ADB well_segment_densities_;
+
+        // the hydrostatic pressure drop between segment nodes
+        // calculated with the above density of fluid mixtures
+        // for the top segment, they should always be zero for the moment.
+        ADB well_segment_pressures_delta_;
+
+        // the surface volume of components in the segments
+        // the initial value at the beginning of the time step
+        std::vector<Vector>   segment_comp_surf_volume_initial_;
+
+        // the value within the current iteration.
+        std::vector<ADB> segment_comp_surf_volume_current_;
+
+        // the mass flow rate in the segments
+        ADB segment_mass_flow_rates_;
+
+        // the viscosity of the fluid mixture in the segments
+        // TODO: it is only used to calculate the Reynolds number as we know
+        //       maybe it is not better just to store the Reynolds number?
+        ADB segment_viscosities_;
+
+        std::vector<int> top_well_segments_;
+
+        // segment volume by dt (time step)
+        // to handle the volume effects of the segment
+        Vector segvdt_;
+
+    };
+
+} // namespace Opm
+
+#include "MultisegmentWells_impl.hpp"
+
+#endif // OPM_MULTISEGMENTWELLS_HEADER_INCLUDED

opm/autodiff/MultisegmentWells_impl.hpp

@@ -0,0 +1,838 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 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_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED
+#define OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED
+
+
+namespace Opm
+{
+
+
+
+    namespace wellhelpers {
+
+        using ADB = MultisegmentWells::ADB;
+        using Vector = MultisegmentWells::Vector;
+
+        inline
+        ADB onlyWellDerivs(const ADB& x)
+        {
+            Vector val = x.value();
+            const int nb = x.numBlocks();
+            if (nb < 2) {
+                OPM_THROW(std::logic_error, "Called onlyWellDerivs() with argument that has " << nb << " blocks.");
+            }
+            std::vector<ADB::M> derivs = { x.derivative()[nb - 2], x.derivative()[nb - 1] };
+            return ADB::function(std::move(val), std::move(derivs));
+        }
+    }
+
+
+
+    template <class WellState>
+    void
+    MultisegmentWells::
+    updateWellState(const Vector& dwells,
+                    const double dpmaxrel,
+                    WellState& well_state) const
+    {
+        if (!wells().empty())
+        {
+            const int nw = wells().size();
+            const int nseg_total = nseg_total_;
+            const int np = numPhases();
+
+            // Extract parts of dwells corresponding to each part.
+            int varstart = 0;
+            const Vector dsegqs = subset(dwells, Span(np * nseg_total, 1, varstart));
+            varstart += dsegqs.size();
+            const Vector dsegp = subset(dwells, Span(nseg_total, 1, varstart));
+            varstart += dsegp.size();
+            assert(varstart == dwells.size());
+
+
+            // segment phase rates update
+            // in dwells, the phase rates are ordered by phase.
+            // while in WellStateMultiSegment, the phase rates are ordered by segments
+            const DataBlock wsr = Eigen::Map<const DataBlock>(dsegqs.data(), np, nseg_total).transpose();
+            const Vector dwsr = Eigen::Map<const Vector>(wsr.data(), nseg_total * np);
+            const Vector wsr_old = Eigen::Map<const Vector>(&well_state.segPhaseRates()[0], nseg_total * np);
+            const Vector sr = wsr_old - dwsr;
+            std::copy(&sr[0], &sr[0] + sr.size(), well_state.segPhaseRates().begin());
+
+
+            // segment pressure updates
+            const Vector segp_old = Eigen::Map<const Vector>(&well_state.segPress()[0], nseg_total, 1);
+            // TODO: applying the pressure change limiter to all the segments, not sure if it is the correct thing to do
+            const Vector dsegp_limited = sign(dsegp) * dsegp.abs().min(segp_old.abs() * dpmaxrel);
+            const Vector segp = segp_old - dsegp_limited;
+            std::copy(&segp[0], &segp[0] + segp.size(), well_state.segPress().begin());
+
+            // update the well rates and bhps, which are not anymore primary vabriables.
+            // they are updated directly from the updated segment phase rates and segment pressures.
+
+            // Bhp update.
+            Vector bhp = Vector::Zero(nw);
+            Vector wr = Vector::Zero(nw * np);
+            // it is better to use subset
+
+            int start_segment = 0;
+            for (int w = 0; w < nw; ++w) {
+                bhp[w] = well_state.segPress()[start_segment];
+                // insert can be faster
+                for (int p = 0; p < np; ++p) {
+                    wr[p + np * w] = well_state.segPhaseRates()[p + np * start_segment];
+                }
+
+                const int nseg = wells()[w]->numberOfSegments();
+                start_segment += nseg;
+            }
+
+            assert(start_segment == nseg_total);
+            std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
+            std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
+
+            // TODO: handling the THP control related.
+        }
+    }
+
+
+
+
+
+    template <class SolutionState>
+    void
+    MultisegmentWells::
+    computeWellFlux(const SolutionState& state,
+                    const Opm::PhaseUsage& pu,
+                    const std::vector<bool>& active,
+                    const Vector& well_perforation_pressure_diffs,
+                    const DataBlock& compi,
+                    const std::vector<ADB>& mob_perfcells,
+                    const std::vector<ADB>& b_perfcells,
+                    Vector& aliveWells,
+                    std::vector<ADB>& cq_s) const
+    {
+        if (wells().size() == 0) return;
+
+        const int np = numPhases();
+        const int nw = wells().size();
+
+        aliveWells = Vector::Constant(nw, 1.0);
+
+        const int nseg = nseg_total_;
+        const int nperf = nperf_total_;
+
+        cq_s.resize(np, ADB::null());
+
+        {
+            const Vector& Tw = wellOps().conn_trans_factors;
+            const std::vector<int>& well_cells = wellOps().well_cells;
+
+            // determining in-flow (towards well-bore) or out-flow (towards reservoir)
+            // for mutli-segmented wells and non-segmented wells, the calculation of the drawdown are different.
+            const ADB& p_perfcells = subset(state.pressure, well_cells);
+            const ADB& rs_perfcells = subset(state.rs, well_cells);
+            const ADB& rv_perfcells = subset(state.rv, well_cells);
+
+            const ADB& seg_pressures = state.segp;
+
+            const ADB seg_pressures_perf = wellOps().s2p * seg_pressures;
+
+            // Create selector for perforations of multi-segment vs. regular wells.
+            Vector is_multisegment_well(nw);
+            for (int w = 0; w < nw; ++w) {
+                is_multisegment_well[w] = double(wells()[w]->isMultiSegmented());
+            }
+            // Take one flag per well and expand to one flag per perforation.
+            Vector is_multisegment_perf = wellOps().w2p * is_multisegment_well.matrix();
+            Selector<double> msperf_selector(is_multisegment_perf, Selector<double>::NotEqualZero);
+
+            // Compute drawdown.
+            ADB h_nc = msperf_selector.select(well_segment_perforation_pressure_diffs_,
+                                              ADB::constant(well_perforation_pressure_diffs));
+            const Vector h_cj = msperf_selector.select(well_perforation_cell_pressure_diffs_, Vector::Zero(nperf));
+
+            // Special handling for when we are called from solveWellEq().
+            // TODO: restructure to eliminate need for special treatmemt.
+            if ((h_nc.numBlocks() != 0) && (h_nc.numBlocks() != seg_pressures_perf.numBlocks())) {
+                assert(seg_pressures_perf.numBlocks() == 2);
+                assert(h_nc.numBlocks() > 2);
+                h_nc = wellhelpers::onlyWellDerivs(h_nc);
+                assert(h_nc.numBlocks() == 2);
+            }
+
+            ADB drawdown = (p_perfcells + h_cj - seg_pressures_perf - h_nc);
+
+            // selects injection perforations
+            Vector selectInjectingPerforations = Vector::Zero(nperf);
+            // selects producing perforations
+            Vector selectProducingPerforations = Vector::Zero(nperf);
+            for (int c = 0; c < nperf; ++c){
+                if (drawdown.value()[c] < 0)
+                    selectInjectingPerforations[c] = 1;
+                else
+                    selectProducingPerforations[c] = 1;
+            }
+
+            // handling flow into wellbore
+            // maybe there are something to do there make the procedure easier.
+            std::vector<ADB> cq_ps(np, ADB::null());
+            for (int phase = 0; phase < np; ++phase) {
+                const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
+                cq_ps[phase] = b_perfcells[phase] * cq_p;
+            }
+
+            if (active[Oil] && active[Gas]) {
+                const int oilpos = pu.phase_pos[Oil];
+                const int gaspos = pu.phase_pos[Gas];
+                const ADB cq_psOil = cq_ps[oilpos];
+                const ADB cq_psGas = cq_ps[gaspos];
+                cq_ps[gaspos] += rs_perfcells * cq_psOil;
+                cq_ps[oilpos] += rv_perfcells * cq_psGas;
+            }
+
+            // hadling flow out from wellbore
+            ADB total_mob = mob_perfcells[0];
+            for (int phase = 1; phase < np; ++phase) {
+                total_mob += mob_perfcells[phase];
+            }
+
+            // injection perforations total volume rates
+            const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
+
+            // compute wellbore mixture for injecting perforations
+            // The wellbore mixture depends on the inflow from the reservoir
+            // and the well injection rates.
+            // TODO: should this based on the segments?
+            // TODO: for the usual wells, the well rates are the sum of the perforations.
+            // TODO: for multi-segmented wells, the segment rates are not the sum of the perforations.
+
+            // TODO: two options here
+            // TODO: 1. for each segment, only the inflow from the perforations related to this segment are considered.
+            // TODO: 2. for each segment, the inflow from the perforrations related to this segment and also all the inflow
+            // TODO: from the upstreaming sgments and their perforations need to be considered.
+            // TODO: This way can be the more consistent way, while let us begin with the first option. The second option
+            // TODO: involves one operations that are not valid now. (i.e. how to transverse from the leaves to the root,
+            // TODO: although we can begin from the brutal force way)
+
+            // TODO: stop using wells() here.
+            std::vector<ADB> wbq(np, ADB::null());
+            ADB wbqt = ADB::constant(Vector::Zero(nseg));
+
+            for (int phase = 0; phase < np; ++phase) {
+                const ADB& q_ps = wellOps().p2s * cq_ps[phase];
+                const ADB& q_s = subset(state.segqs, Span(nseg, 1, phase * nseg));
+                Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
+
+                const int pos = pu.phase_pos[phase];
+
+                // this is per segment
+                wbq[phase] = (wellOps().w2s * ADB::constant(compi.col(pos)) * injectingPhase_selector.select(q_s, ADB::constant(Vector::Zero(nseg)))) - q_ps;
+
+                // TODO: it should be a single value for this certain well.
+                // TODO: it need to be changed later to handle things more consistently
+                // or there should be an earsier way to decide if the well is dead.
+                wbqt += wbq[phase];
+            }
+
+            // Set aliveWells.
+            // the first value of the wbqt is the one to decide if the well is dead
+            // or there should be some dead segments?
+            {
+                int topseg = 0;
+                for (int w = 0; w < nw; ++w) {
+                    if (wbqt.value()[topseg] == 0.0) { // yes we really mean == here, no fuzzyness
+                        aliveWells[w] = 0.0;
+                    }
+                    topseg += wells()[w]->numberOfSegments();
+                }
+            }
+
+            // compute wellbore mixture at standard conditions.
+            // before, the determination of alive wells is based on wells.
+            // now, will there be any dead segment? I think no.
+            // TODO: it is not clear if the cmix_s should be based on segment or the well
+            std::vector<ADB> cmix_s(np, ADB::null());
+            Selector<double> aliveWells_selector(aliveWells, Selector<double>::NotEqualZero);
+            for (int phase = 0; phase < np; ++phase) {
+                const int pos = pu.phase_pos[phase];
+                const ADB phase_fraction = wellOps().topseg2w * (wbq[phase] / wbqt);
+                cmix_s[phase] = wellOps().w2p * aliveWells_selector.select(phase_fraction, ADB::constant(compi.col(pos)));
+            }
+
+            // compute volume ration between connection at standard conditions
+            ADB volumeRatio = ADB::constant(Vector::Zero(nperf));
+            const ADB d = Vector::Constant(nperf,1.0) -  rv_perfcells * rs_perfcells;
+
+            for (int phase = 0; phase < np; ++phase) {
+                ADB tmp = cmix_s[phase];
+                if (phase == Oil && active[Gas]) {
+                    const int gaspos = pu.phase_pos[Gas];
+                    tmp = tmp - rv_perfcells * cmix_s[gaspos] / d;
+                }
+                if (phase == Gas && active[Oil]) {
+                    const int oilpos = pu.phase_pos[Oil];
+                    tmp = tmp - rs_perfcells * cmix_s[oilpos] / d;
+                }
+                volumeRatio += tmp / b_perfcells[phase];
+            }
+
+            // injecting connections total volumerates at standard conditions
+            ADB cqt_is = cqt_i/volumeRatio;
+
+            // connection phase volumerates at standard conditions
+            for (int phase = 0; phase < np; ++phase) {
+                cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
+            }
+        }
+    }
+
+
+
+
+
+    template <class SolutionState>
+    void
+    MultisegmentWells::
+    computeSegmentFluidProperties(const SolutionState& state,
+                                  const std::vector<PhasePresence>& pc,
+                                  const std::vector<bool>& active,
+                                  const BlackoilPropsAdInterface& fluid)
+    {
+        const int np = numPhases();
+        const int nw = wells().size();
+        const int nseg_total = nseg_total_;
+
+        if ( !wellOps().has_multisegment_wells ){
+            // not sure if this is needed actually
+            // TODO: to check later if this is really necessary.
+            well_segment_densities_ = ADB::constant(Vector::Zero(nseg_total));
+            segment_mass_flow_rates_ = ADB::constant(Vector::Zero(nseg_total));
+            segment_viscosities_ = ADB::constant(Vector::Zero(nseg_total));
+            for (int phase = 0; phase < np; ++phase) {
+                segment_comp_surf_volume_current_[phase] = ADB::constant(Vector::Zero(nseg_total));
+                segmentCompSurfVolumeInitial()[phase] = Vector::Zero(nseg_total);
+            }
+            return;
+        }
+
+        // although we will calculate segment density for non-segmented wells at the same time,
+        // while under most of the cases, they will not be used,
+        // since for most of the cases, the density calculation for non-segment wells are
+        // set to be 'SEG' way, which is not a option for multi-segment wells.
+        // When the density calcuation for non-segmented wells are set to 'AVG', then
+        // the density calculation of the mixtures can be the same, while it remains to be verified.
+
+        // The grid cells associated with segments.
+        // TODO: shoud be computed once and stored in WellState or global Wells structure or class.
+        std::vector<int> segment_cells;
+        segment_cells.reserve(nseg_total);
+        for (int w = 0; w < nw; ++w) {
+            const std::vector<int>& segment_cells_well = wells()[w]->segmentCells();
+            segment_cells.insert(segment_cells.end(), segment_cells_well.begin(), segment_cells_well.end());
+        }
+        assert(int(segment_cells.size()) == nseg_total);
+
+        const ADB segment_temp = subset(state.temperature, segment_cells);
+        // using the segment pressure or the average pressure
+        // using the segment pressure first
+        const ADB& segment_press = state.segp;
+
+        // Compute PVT properties for segments.
+        std::vector<PhasePresence> segment_cond(nseg_total);
+        for (int s = 0; s < nseg_total; ++s) {
+            segment_cond[s] = pc[segment_cells[s]];
+        }
+        std::vector<ADB> b_seg(np, ADB::null());
+        // Viscosities for different phases
+        std::vector<ADB> mu_seg(np, ADB::null());
+        ADB rsmax_seg = ADB::null();
+        ADB rvmax_seg = ADB::null();
+        const PhaseUsage& pu = fluid.phaseUsage();
+        if (pu.phase_used[Water]) {
+            b_seg[pu.phase_pos[Water]] = fluid.bWat(segment_press, segment_temp, segment_cells);
+            mu_seg[pu.phase_pos[Water]] = fluid.muWat(segment_press, segment_temp, segment_cells);
+        }
+        assert(active[Oil]);
+        const ADB segment_so = subset(state.saturation[pu.phase_pos[Oil]], segment_cells);
+        if (pu.phase_used[Oil]) {
+            const ADB segment_rs = subset(state.rs, segment_cells);
+            b_seg[pu.phase_pos[Oil]] = fluid.bOil(segment_press, segment_temp, segment_rs,
+                                                   segment_cond, segment_cells);
+            // rsmax_seg = fluidRsSat(segment_press, segment_so, segment_cells);
+            rsmax_seg = fluid.rsSat(segment_press, segment_so, segment_cells);
+            mu_seg[pu.phase_pos[Oil]] = fluid.muOil(segment_press, segment_temp, segment_rs,
+                                                     segment_cond, segment_cells);
+        }
+        assert(active[Gas]);
+        if (pu.phase_used[Gas]) {
+            const ADB segment_rv = subset(state.rv, segment_cells);
+            b_seg[pu.phase_pos[Gas]] = fluid.bGas(segment_press, segment_temp, segment_rv,
+                                                   segment_cond, segment_cells);
+            // rvmax_seg = fluidRvSat(segment_press, segment_so, segment_cells);
+            rvmax_seg = fluid.rvSat(segment_press, segment_so, segment_cells);
+            mu_seg[pu.phase_pos[Gas]] = fluid.muGas(segment_press, segment_temp, segment_rv,
+                                                   segment_cond, segment_cells);
+        }
+
+        // Extract segment flow by phase (segqs) and compute total surface rate.
+        ADB tot_surface_rate = ADB::constant(Vector::Zero(nseg_total));
+        std::vector<ADB> segqs(np, ADB::null());
+        for (int phase = 0; phase < np; ++phase) {
+            segqs[phase] = subset(state.segqs, Span(nseg_total, 1, phase * nseg_total));
+            tot_surface_rate += segqs[phase];
+        }
+
+        // TODO: later this will be implmented as a global mapping
+        std::vector<std::vector<double>> comp_frac(np, std::vector<double>(nseg_total, 0.0));
+        int start_segment = 0;
+        for (int w = 0; w < nw; ++w) {
+            WellMultiSegmentConstPtr well = wells()[w];
+            const int nseg = well->numberOfSegments();
+            const std::vector<double>& comp_frac_well = well->compFrac();
+            for (int phase = 0; phase < np; ++phase) {
+                for (int s = 0; s < nseg; ++s) {
+                    comp_frac[phase][s + start_segment] = comp_frac_well[phase];
+                }
+            }
+            start_segment += nseg;
+        }
+        assert(start_segment == nseg_total);
+
+        // Compute mix.
+        // 'mix' contains the component fractions under surface conditions.
+        std::vector<ADB> mix(np, ADB::null());
+        for (int phase = 0; phase < np; ++phase) {
+            // initialize to be the compFrac for each well,
+            // then update only the one with non-zero total volume rate
+            mix[phase] = ADB::constant(Eigen::Map<Vector>(comp_frac[phase].data(), nseg_total));
+        }
+        // There should be a better way to do this.
+        Selector<double> non_zero_tot_rate(tot_surface_rate.value(), Selector<double>::NotEqualZero);
+        for (int phase = 0; phase < np; ++phase) {
+            mix[phase] = non_zero_tot_rate.select(segqs[phase] / tot_surface_rate, mix[phase]);
+        }
+
+        // Calculate rs and rv.
+        ADB rs = ADB::constant(Vector::Zero(nseg_total));
+        ADB rv = rs;
+        const int gaspos = pu.phase_pos[Gas];
+        const int oilpos = pu.phase_pos[Oil];
+        Selector<double> non_zero_mix_oilpos(mix[oilpos].value(), Selector<double>::GreaterZero);
+        Selector<double> non_zero_mix_gaspos(mix[gaspos].value(), Selector<double>::GreaterZero);
+        // What is the better way to do this?
+        // big values should not be necessary
+        ADB big_values = ADB::constant(Vector::Constant(nseg_total, 1.e100));
+        ADB mix_gas_oil = non_zero_mix_oilpos.select(mix[gaspos] / mix[oilpos], big_values);
+        ADB mix_oil_gas = non_zero_mix_gaspos.select(mix[oilpos] / mix[gaspos], big_values);
+        if (active[Oil]) {
+            Vector selectorUnderRsmax = Vector::Zero(nseg_total);
+            Vector selectorAboveRsmax = Vector::Zero(nseg_total);
+            for (int s = 0; s < nseg_total; ++s) {
+                if (mix_gas_oil.value()[s] > rsmax_seg.value()[s]) {
+                    selectorAboveRsmax[s] = 1.0;
+                } else {
+                    selectorUnderRsmax[s] = 1.0;
+                }
+            }
+            rs = non_zero_mix_oilpos.select(selectorAboveRsmax * rsmax_seg + selectorUnderRsmax * mix_gas_oil, rs);
+        }
+        if (active[Gas]) {
+            Vector selectorUnderRvmax = Vector::Zero(nseg_total);
+            Vector selectorAboveRvmax = Vector::Zero(nseg_total);
+            for (int s = 0; s < nseg_total; ++s) {
+                if (mix_oil_gas.value()[s] > rvmax_seg.value()[s]) {
+                    selectorAboveRvmax[s] = 1.0;
+                } else {
+                    selectorUnderRvmax[s] = 1.0;
+                }
+            }
+            rv = non_zero_mix_gaspos.select(selectorAboveRvmax * rvmax_seg + selectorUnderRvmax * mix_oil_gas, rv);
+        }
+
+        // Calculate the phase fraction under reservoir conditions.
+        std::vector<ADB> x(np, ADB::null());
+        for (int phase = 0; phase < np; ++phase) {
+            x[phase] = mix[phase];
+        }
+        if (active[Gas] && active[Oil]) {
+            x[gaspos] = (mix[gaspos] - mix[oilpos] * rs) / (Vector::Ones(nseg_total) - rs * rv);
+            x[oilpos] = (mix[oilpos] - mix[gaspos] * rv) / (Vector::Ones(nseg_total) - rs * rv);
+        }
+
+        // Compute total reservoir volume to surface volume ratio.
+        ADB volrat = ADB::constant(Vector::Zero(nseg_total));
+        for (int phase = 0; phase < np; ++phase) {
+            volrat += x[phase] / b_seg[phase];
+        }
+
+        // Compute segment densities.
+        ADB dens = ADB::constant(Vector::Zero(nseg_total));
+        for (int phase = 0; phase < np; ++phase) {
+            const Vector surface_density = fluid.surfaceDensity(phase, segment_cells);
+            dens += surface_density * mix[phase];
+        }
+        well_segment_densities_ = dens / volrat;
+
+        // Calculating the surface volume of each component in the segment
+        assert(np == int(segment_comp_surf_volume_current_.size()));
+        const ADB segment_surface_volume = segvdt_ / volrat;
+        for (int phase = 0; phase < np; ++phase) {
+            segment_comp_surf_volume_current_[phase] = segment_surface_volume * mix[phase];
+        }
+
+        // Mass flow rate of the segments
+        segment_mass_flow_rates_ = ADB::constant(Vector::Zero(nseg_total));
+        for (int phase = 0; phase < np; ++phase) {
+            // TODO: how to remove one repeated surfaceDensity()
+            const Vector surface_density = fluid.surfaceDensity(phase, segment_cells);
+            segment_mass_flow_rates_ += surface_density * segqs[phase];
+        }
+
+        // Viscosity of the fluid mixture in the segments
+        segment_viscosities_ = ADB::constant(Vector::Zero(nseg_total));
+        for (int phase = 0; phase < np; ++phase) {
+            segment_viscosities_ += x[phase] * mu_seg[phase];
+        }
+    }
+
+
+
+
+
+    template <class SolutionState>
+    void
+    MultisegmentWells::
+    addWellFluxEq(const std::vector<ADB>& cq_s,
+                  const SolutionState& state,
+                  LinearisedBlackoilResidual& residual)
+    {
+        // the well flux equations are for each segment and each phase.
+        //    /delta m_p_n / dt  - /sigma Q_pi - /sigma q_pj + Q_pn = 0
+        // 1. It is the gain of the amount of the component p in the segment n during the
+        //    current time step under stock-tank conditions.
+        //    It is used to handle the volume storage effects of the wellbore.
+        //    We need the information from the previous step and the crrent time step.
+        // 2. for the second term, it is flow into the segment from the inlet segments,
+        //    which are unknown and treated implictly.
+        // 3. for the third term, it is the inflow through the perforations.
+        // 4. for the last term, it is the outlet rates and also the segment rates,
+        //    which are the primary variable.
+        const int np = numPhases();
+        const int nseg_total = nseg_total_;
+
+        ADB segqs = state.segqs;
+
+        std::vector<ADB> segment_volume_change_dt(np, ADB::null());
+        for (int phase = 0; phase < np; ++phase) {
+            if ( wellOps().has_multisegment_wells ) {
+                // Gain of the surface volume of each component in the segment by dt
+                segment_volume_change_dt[phase] = segment_comp_surf_volume_current_[phase] -
+                                                  segmentCompSurfVolumeInitial()[phase];
+
+                // Special handling for when we are called from solveWellEq().
+                // TODO: restructure to eliminate need for special treatmemt.
+                if (segment_volume_change_dt[phase].numBlocks() != segqs.numBlocks()) {
+                    assert(segment_volume_change_dt[phase].numBlocks() > 2);
+                    assert(segqs.numBlocks() == 2);
+                    segment_volume_change_dt[phase] = wellhelpers::onlyWellDerivs(segment_volume_change_dt[phase]);
+                    assert(segment_volume_change_dt[phase].numBlocks() == 2);
+                }
+
+                const ADB cq_s_seg = wellOps().p2s * cq_s[phase];
+                const ADB segqs_phase = subset(segqs, Span(nseg_total, 1, phase * nseg_total));
+                segqs -= superset(cq_s_seg + wellOps().s2s_inlets * segqs_phase + segment_volume_change_dt[phase],
+                                  Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
+            } else {
+                segqs -= superset(wellOps().p2s * cq_s[phase], Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
+            }
+        }
+
+        residual.well_flux_eq = segqs;
+    }
+
+
+
+
+
+    template <class SolutionState, class WellState>
+    void
+    MultisegmentWells::
+    addWellControlEq(const SolutionState& state,
+                     const WellState& xw,
+                     const Vector& aliveWells,
+                     const std::vector<bool>& active,
+                     LinearisedBlackoilResidual& residual)
+    {
+        // the name of the function is a a little misleading.
+        // Basically it is the function for the pressure equation.
+        // And also, it work as the control equation when it is the segment
+        if( wells().empty() ) return;
+
+        const int np = numPhases();
+        const int nw = wells().size();
+        const int nseg_total = nseg_total_;
+
+        ADB aqua   = ADB::constant(Vector::Zero(nseg_total));
+        ADB liquid = ADB::constant(Vector::Zero(nseg_total));
+        ADB vapour = ADB::constant(Vector::Zero(nseg_total));
+
+        if (active[Water]) {
+            aqua += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Aqua * nseg_total));
+        }
+        if (active[Oil]) {
+            liquid += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Liquid * nseg_total));
+        }
+        if (active[Gas]) {
+            vapour += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Vapour * nseg_total));
+        }
+
+        // THP control is not implemented for the moment.
+
+        // Hydrostatic correction variables
+        Vector rho_v = Vector::Zero(nw);
+        Vector vfp_ref_depth_v = Vector::Zero(nw);
+
+        // Target vars
+        Vector bhp_targets  = Vector::Zero(nw);
+        Vector rate_targets = Vector::Zero(nw);
+        Eigen::SparseMatrix<double>  rate_distr(nw, np*nw);
+
+        // Selection variables
+        // well selectors
+        std::vector<int> bhp_well_elems;
+        std::vector<int> rate_well_elems;
+        // segment selectors
+        std::vector<int> bhp_top_elems;
+        std::vector<int> rate_top_elems;
+        std::vector<int> rate_top_phase_elems;
+        std::vector<int> others_elems;
+
+        //Run through all wells to calculate BHP/RATE targets
+        //and gather info about current control
+        int start_segment = 0;
+        for (int w = 0; w < nw; ++w) {
+            const struct WellControls* wc = wells()[w]->wellControls();
+
+            // The current control in the well state overrides
+            // the current control set in the Wells struct, which
+            // is instead treated as a default.
+            const int current = xw.currentControls()[w];
+
+            const int nseg = wells()[w]->numberOfSegments();
+
+            switch (well_controls_iget_type(wc, current)) {
+            case BHP:
+            {
+                bhp_well_elems.push_back(w);
+                bhp_top_elems.push_back(start_segment);
+                bhp_targets(w)  = well_controls_iget_target(wc, current);
+                rate_targets(w) = -1e100;
+                for (int p = 0; p < np; ++p) {
+                    rate_top_phase_elems.push_back(np * start_segment + p);
+                }
+            }
+            break;
+
+            case THP:
+            {
+                OPM_THROW(std::runtime_error, "THP control is not implemented for multi-sgement wells yet!!");
+            }
+            break;
+
+            case RESERVOIR_RATE: // Intentional fall-through
+            case SURFACE_RATE:
+            {
+                rate_well_elems.push_back(w);
+                rate_top_elems.push_back(start_segment);
+                for (int p = 0; p < np; ++p) {
+                    rate_top_phase_elems.push_back(np * start_segment + p);
+                }
+                // RESERVOIR and SURFACE rates look the same, from a
+                // high-level point of view, in the system of
+                // simultaneous linear equations.
+
+                const double* const distr =
+                    well_controls_iget_distr(wc, current);
+
+                for (int p = 0; p < np; ++p) {
+                    rate_distr.insert(w, p*nw + w) = distr[p];
+                }
+
+                bhp_targets(w)  = -1.0e100;
+                rate_targets(w) = well_controls_iget_target(wc, current);
+            }
+            break;
+
+            }
+
+            for (int i = 1; i < nseg; ++i) {
+                others_elems.push_back(i + start_segment);
+            }
+            start_segment += nseg;
+        }
+
+        // for each segment: 1, if the segment is the top segment, then control equation
+        //                   2, if the segment is not the top segment, then the pressure equation
+        const ADB bhp_residual = subset(state.segp, bhp_top_elems) - subset(bhp_targets, bhp_well_elems);
+        const ADB rate_residual = subset(rate_distr * subset(state.segqs, rate_top_phase_elems) - rate_targets, rate_well_elems);
+
+        ADB others_residual = ADB::constant(Vector::Zero(nseg_total));
+
+        if ( wellOps().has_multisegment_wells ) {
+            // Special handling for when we are called from solveWellEq().
+            // TODO: restructure to eliminate need for special treatmemt.
+            ADB wspd = (state.segp.numBlocks() == 2)
+                ? wellhelpers::onlyWellDerivs(well_segment_pressures_delta_)
+                : well_segment_pressures_delta_;
+
+            others_residual = wellOps().eliminate_topseg * (state.segp - wellOps().s2s_outlet * state.segp + wspd);
+        } else {
+            others_residual = wellOps().eliminate_topseg * (state.segp - wellOps().s2s_outlet * state.segp);
+        }
+
+        //       all the control equations
+        // TODO: can be optimized better
+        ADB well_eq_topsegment = subset(superset(bhp_residual, bhp_top_elems, nseg_total) +
+                                        superset(rate_residual, rate_top_elems, nseg_total), top_well_segments_);
+
+        // For wells that are dead (not flowing), and therefore not communicating
+        // with the reservoir, we set the equation to be equal to the well's total
+        // flow. This will be a solution only if the target rate is also zero.
+        Eigen::SparseMatrix<double> rate_summer(nw, np*nw);
+        for (int w = 0; w < nw; ++w) {
+            for (int phase = 0; phase < np; ++phase) {
+                rate_summer.insert(w, phase*nw + w) = 1.0;
+            }
+        }
+        Selector<double> alive_selector(aliveWells, Selector<double>::NotEqualZero);
+        // TODO: Here only handles the wells, or the top segments
+        // should we also handle some non-alive non-top segments?
+        // should we introduce the cocept of non-alive segments?
+        // At the moment, we only handle the control equations
+        well_eq_topsegment = alive_selector.select(well_eq_topsegment, rate_summer * subset(state.segqs, rate_top_phase_elems));
+
+        /* residual_.well_eq = superset(bhp_residual, bhp_top_elems, nseg_total) +
+                            superset(rate_residual, rate_top_elems, nseg_total) +
+                            superset(others_residual, others_elems, nseg_total); */
+        residual.well_eq = superset(well_eq_topsegment, top_well_segments_, nseg_total) +
+                            others_residual;
+    }
+
+
+
+
+
+    template <class WellState>
+    void
+    MultisegmentWells::
+    updateWellControls(const bool terminal_output,
+                       WellState& xw) const
+    {
+        if( wells().empty() ) return ;
+
+        std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" };
+        // Find, for each well, if any constraints are broken. If so,
+        // switch control to first broken constraint.
+        const int np = numPhases();
+        const int nw = wells().size();
+        for (int w = 0; w < nw; ++w) {
+            const WellControls* wc = wells()[w]->wellControls();
+            // The current control in the well state overrides
+            // the current control set in the Wells struct, which
+            // is instead treated as a default.
+            int current = xw.currentControls()[w];
+            // Loop over all controls except the current one, and also
+            // skip any RESERVOIR_RATE controls, since we cannot
+            // handle those.
+            const int nwc = well_controls_get_num(wc);
+            int ctrl_index = 0;
+            for (; ctrl_index < nwc; ++ctrl_index) {
+                if (ctrl_index == current) {
+                    // This is the currently used control, so it is
+                    // used as an equation. So this is not used as an
+                    // inequality constraint, and therefore skipped.
+                    continue;
+                }
+                if (wellhelpers::constraintBroken(
+                        xw.bhp(), xw.thp(), xw.wellRates(),
+                        w, np, wells()[w]->wellType(), wc, ctrl_index)) {
+                    // ctrl_index will be the index of the broken constraint after the loop.
+                    break;
+                }
+            }
+
+            if (ctrl_index != nwc) {
+                // Constraint number ctrl_index was broken, switch to it.
+                if (terminal_output)
+                {
+                    std::cout << "Switching control mode for well " << wells()[w]->name()
+                              << " from " << modestring[well_controls_iget_type(wc, current)]
+                              << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
+                }
+                xw.currentControls()[w] = ctrl_index;
+                current = xw.currentControls()[w];
+            }
+
+            // Get gravity for THP hydrostatic corrrection
+            // const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
+
+            // Updating well state and primary variables.
+            // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
+            const double target = well_controls_iget_target(wc, current);
+            const double* distr = well_controls_iget_distr(wc, current);
+            switch (well_controls_iget_type(wc, current)) {
+            case BHP:
+                xw.bhp()[w] = target;
+                xw.segPress()[top_well_segments_[w]] = target;
+                break;
+
+            case THP: {
+                OPM_THROW(std::runtime_error, "THP control is not implemented for multi-sgement wells yet!!");
+            }
+
+            case RESERVOIR_RATE:
+                // No direct change to any observable quantity at
+                // surface condition.  In this case, use existing
+                // flow rates as initial conditions as reservoir
+                // rate acts only in aggregate.
+                break;
+
+            case SURFACE_RATE:
+                for (int phase = 0; phase < np; ++phase) {
+                    if (distr[phase] > 0.0) {
+                        xw.wellRates()[np * w + phase] = target * distr[phase];
+                        // TODO: consider changing all (not just top) segment rates
+                        // to make them consistent, it could possibly improve convergence.
+                        xw.segPhaseRates()[np * xw.topSegmentLoc()[w] + phase] = target * distr[phase];
+                    }
+                }
+                break;
+            }
+
+        }
+
+    }
+
+
+}
+#endif // OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED

tests/test_multisegmentwells.cpp

@@ -0,0 +1,181 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 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/>.
+*/
+
+
+#include <config.h>
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define BOOST_TEST_MODULE MultisegmentWellsTest
+
+#include <vector>
+#include <memory>
+#include <array>
+
+
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <boost/test/unit_test.hpp>
+#include <boost/filesystem.hpp>
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/wells.h>
+
+#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
+#include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/createGlobalCellArray.hpp>
+#include <opm/autodiff/GridInit.hpp>
+
+#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
+#include <opm/autodiff/MultisegmentWells.hpp>
+
+
+
+struct SetupMSW {
+
+    using Grid = UnstructuredGrid;
+    using GridInit = Opm::GridInit<Grid>;
+    using FluidProps = Opm::BlackoilPropsAdFromDeck;
+    using MaterialLawManager = FluidProps::MaterialLawManager;
+
+    SetupMSW()
+    {
+        Opm::ParseContext parse_context;
+        Opm::ParserPtr parser(new Opm::Parser());
+        Opm::DeckConstPtr deck = parser->parseFile("msw.data", parse_context);
+        Opm::EclipseStateConstPtr ecl_state(new Opm::EclipseState(deck , parse_context));
+
+        // Create grid.
+        const std::vector<double>& porv =
+                            ecl_state->get3DProperties().getDoubleGridProperty("PORV").getData();
+
+        std::unique_ptr<GridInit> grid_init(new GridInit(ecl_state, porv));
+        const Grid& grid = grid_init->grid();
+
+        // Create material law manager.
+        std::vector<int> compressed_to_cartesianIdx;
+        Opm::createGlobalCellArray(grid, compressed_to_cartesianIdx);
+
+        std::shared_ptr<MaterialLawManager> material_law_manager(new MaterialLawManager());
+        material_law_manager->initFromDeck(deck, ecl_state, compressed_to_cartesianIdx);
+
+        std::unique_ptr<FluidProps> fluidprops(new FluidProps(deck, ecl_state, material_law_manager, grid));
+
+        const size_t current_timestep = 0;
+
+        // Create wells.
+        Opm::WellsManager wells_manager(ecl_state,
+                                        current_timestep,
+                                        Opm::UgGridHelpers::numCells(grid),
+                                        Opm::UgGridHelpers::globalCell(grid),
+                                        Opm::UgGridHelpers::cartDims(grid),
+                                        Opm::UgGridHelpers::dimensions(grid),
+                                        Opm::UgGridHelpers::cell2Faces(grid),
+                                        Opm::UgGridHelpers::beginFaceCentroids(grid),
+                                        fluidprops->permeability(),
+                                        false);
+
+        const Wells* wells = wells_manager.c_wells();
+        const std::vector<Opm::WellConstPtr>& wells_ecl = ecl_state->getSchedule()->getWells(current_timestep);
+
+        std::vector<Opm::WellMultiSegmentConstPtr> wells_multisegment;
+        wells_multisegment.reserve(wells_ecl.size());
+        for (size_t i = 0; i < wells_ecl.size(); ++i) {
+            // not processing SHUT wells.
+            if (wells_ecl[i]->getStatus(current_timestep) == Opm::WellCommon::SHUT) {
+                continue;
+            }
+            // checking if the well can be found in the wells
+            const std::string& well_name = wells_ecl[i]->name();
+            // number of wells in wells
+            const int nw_wells = wells->number_of_wells;
+            int index_well;
+            for (index_well = 0; index_well < nw_wells; ++index_well) {
+                if (well_name == std::string(wells->name[index_well])) {
+                    break;
+                }
+            }
+
+            // have to be able to be found
+            assert(index_well != nw_wells);
+
+            wells_multisegment.push_back(std::make_shared<Opm::WellMultiSegment>(wells_ecl[i], current_timestep, wells));
+        }
+
+        ms_wells.reset(new Opm::MultisegmentWells(wells_multisegment, wells->number_of_phases));
+    };
+
+    std::shared_ptr<const Opm::MultisegmentWells> ms_wells;
+};
+
+// number of wells for this case
+const int nw = 2;
+// number of segments for this case
+const int nseg = 7;
+// number of perforations for this case
+const int nperf = 8;
+
+BOOST_AUTO_TEST_CASE(testOperators)
+{
+    SetupMSW msw_setup;
+    const std::shared_ptr<const Opm::MultisegmentWells>& ms_wells = msw_setup.ms_wells;
+
+    const Opm::MultisegmentWells::MultisegmentWellOps& wops_ms = ms_wells->wellOps();
+    BOOST_CHECK_EQUAL(true, wops_ms.has_multisegment_wells);
+    BOOST_CHECK_EQUAL(nperf, wops_ms.well_cells.size());
+    BOOST_CHECK_EQUAL(nperf, wops_ms.conn_trans_factors.size());
+    BOOST_CHECK_EQUAL(nseg, wops_ms.s2s_outlet.rows());
+    BOOST_CHECK_EQUAL(nseg, wops_ms.s2s_inlets.cols());
+    BOOST_CHECK_EQUAL(nw, wops_ms.s2w.rows());
+    BOOST_CHECK_EQUAL(nseg, wops_ms.s2w.cols());
+    BOOST_CHECK_EQUAL(nseg, wops_ms.topseg2w.cols());
+    BOOST_CHECK_EQUAL(nperf, wops_ms.s2p.rows());
+    BOOST_CHECK_EQUAL(nseg, wops_ms.p2s.rows());
+    BOOST_CHECK_EQUAL(nperf, wops_ms.w2p.rows());
+    BOOST_CHECK_EQUAL(nw, wops_ms.w2p.cols());
+    BOOST_CHECK_EQUAL(nperf, wops_ms.p2w.cols());
+}
+
+
+BOOST_AUTO_TEST_CASE(testStructure)
+{
+    SetupMSW msw_setup;
+    const std::shared_ptr<const Opm::MultisegmentWells>& ms_wells = msw_setup.ms_wells;
+
+    BOOST_CHECK_EQUAL(3, ms_wells->numPhases());
+    BOOST_CHECK_EQUAL(nseg, ms_wells->numSegment());
+    BOOST_CHECK_EQUAL(nperf, ms_wells->numPerf());
+
+    BOOST_CHECK_EQUAL(nw, ms_wells->topWellSegments().size());
+    BOOST_CHECK_EQUAL(nw, ms_wells->wells().size());
+    BOOST_CHECK_EQUAL(0, ms_wells->topWellSegments()[0]);
+    BOOST_CHECK_EQUAL(1, ms_wells->topWellSegments()[1]);
+}