Can now solve with bhp-controlled wells.

Simple initial code. Assumes that well_state.bhp() contains
well bhp targets, does not check control structures.

opm/autodiff/ImpesTPFAAD.hpp

@@ -33,6 +33,7 @@
 #include <algorithm>
 #include <cassert>
 #include <vector>
+#include <iterator>
 
 #include <boost/shared_ptr.hpp>
 
@@ -129,9 +130,24 @@ namespace Opm {
         {
             const int nc = grid_.number_of_cells;
             const int np = state.numPhases();
+
+            // Compute relperms once and for all (since saturations are explicit).
             DataBlock s = Eigen::Map<const DataBlock>(state.saturation().data(), nc, np);
             ASSERT(np == 2);
             kr_ = fluid_.relperm(s.col(0), s.col(1), V::Zero(nc,1), buildAllCells(nc));
+            // Compute relperms for wells. This must be revisited for crossflow.
+            const int nw = wells_.number_of_wells;
+            const int nperf = wells_.well_connpos[nw];
+            DataBlock well_s(nperf, np);
+            for (int w = 0; w < nw; ++w) {
+                const double* comp_frac = &wells_.comp_frac[np*w];
+                for (int j = wells_.well_connpos[w]; j < wells_.well_connpos[w+1]; ++j) {
+                    well_s.row(j) = Eigen::Map<const DataBlock>(comp_frac, 1, np);
+                }
+            }
+            const std::vector<int> well_cells(wells_.well_cells,
+                                              wells_.well_cells + nperf);
+            well_kr_ = fluid_.relperm(well_s.col(0), well_s.col(1), V::Zero(nperf,1), well_cells);
 
             const double atol  = 1.0e-15;
             const double rtol  = 5.0e-10;
@@ -177,22 +193,26 @@ namespace Opm {
                              Eigen::Dynamic,
                              Eigen::RowMajor> DataBlock;
 
-        const UnstructuredGrid& grid_;
+        const UnstructuredGrid&      grid_;
-        const BOFluid&          fluid_;
+        const BOFluid&               fluid_;
-        const GeoProps&         geo_ ;
+        const GeoProps&              geo_ ;
-        const Wells&            wells_;
+        const Wells&                 wells_;
         const LinearSolverInterface& linsolver_;
-        // PDepFData               pdepfdata_;
+        HelperOps                    ops_;
-        HelperOps               ops_;
+        const M                      grav_;
-        const M                 grav_;
+        ADB                          cell_residual_;
-        ADB                     cell_residual_;
+        ADB                          well_residual_;
-        ADB                     well_residual_;
+        std::vector<V>               kr_;
-        std::vector<V> kr_;
+        std::vector<V>               well_kr_;
 
         enum { Water = BOFluid::Water,
                Oil = BOFluid::Oil,
                Gas = BOFluid::Gas };
 
+
+
+
+
         void
         assemble(const double         dt,
                  const BlackoilState& state,
@@ -203,8 +223,7 @@ namespace Opm {
             const int nc = grid_.number_of_cells;
             const int np = state.numPhases();
             const int nw = wells_.number_of_wells;
-
+            const int nperf = wells_.well_connpos[nw];
-            // pdepfdata_.computePressQuant(state);
 
             const std::vector<int> cells = buildAllCells(nc);
 
@@ -213,10 +232,9 @@ namespace Opm {
             const V delta_t = dt * V::Ones(nc, 1);
             const V transi = subset(geo_.transmissibility(),
                                     ops_.internal_faces);
-            const int num_perf = wells_.well_connpos[nw];
             const std::vector<int> well_cells(wells_.well_cells,
-                                              wells_.well_cells + num_perf);
+                                              wells_.well_cells + nperf);
-            const V transw = Eigen::Map<const V>(wells_.WI, num_perf, 1);
+            const V transw = Eigen::Map<const V>(wells_.WI, nperf, 1);
 
             // Initialize AD variables: p (cell pressures) and bhp (well bhp).
             const V p0 = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
@@ -246,10 +264,10 @@ namespace Opm {
             well_to_perf.setFromTriplets(w2p.begin(), w2p.end());
             // Construct pressure difference vector for wells.
             const V well_perf_dp = V::Zero(well_cells.size()); // No gravity yet!
-            // Finally construct well perforation pressures.
+            // Finally construct well perforation pressures and well flows.
             const ADB p_perfwell = well_to_perf*bhp + well_perf_dp;
-
             const ADB nkgradp_well = transw * (p_perfcell - p_perfwell);
+            const Selector<double> cell_to_well_selector(nkgradp_well.value());
 
             cell_residual_ = ADB::constant(pv, bpat);
 #define COMPENSATE_FLOAT_PRECISION 0
@@ -259,22 +277,35 @@ namespace Opm {
             for (int phase = 0; phase < np; ++phase) {
                 const ADB cell_b = fluidFvf(phase, p, cells);
                 const ADB cell_rho = fluidRho(phase, p, cells);
+                const ADB well_b = fluidFvf(phase, p_perfwell, well_cells);
                 const V kr = fluidKr(phase);
+                // Explicitly not asking for derivatives of viscosity,
+                // since they are not available yet.
                 const V mu = fluidMu(phase, p.value(), cells);
-                const V mf = upwind.select(kr / mu);
+                const V cell_mob = kr / mu;
-                const ADB flux = mf * (nkgradp + (grav_ * cell_rho));
+                const V face_mob = upwind.select(cell_mob);
+                const V well_kr = fluidKrWell(phase);
+                const V well_mu = fluidMu(phase, p_perfwell.value(), well_cells);
+                const V well_mob = well_kr / well_mu;
+                const V perf_mob = cell_to_well_selector.select(subset(cell_mob, well_cells), well_mob);
+                const ADB flux = face_mob * (nkgradp + (grav_ * cell_rho));
+                const ADB perf_flux = perf_mob * (nkgradp_well); // No gravity term for perforations.
                 const ADB face_b = upwind.select(cell_b);
+                const ADB perf_b = cell_to_well_selector.select(subset(cell_b, well_cells), well_b);
                 const V z0 = z0all.block(0, phase, nc, 1);
                 const V q  = qall .block(0, phase, nc, 1);
+                const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc);
 #if COMPENSATE_FLOAT_PRECISION
-                const ADB divcontrib = delta_t * (ops_.div * (flux * face_b));
+                const ADB divcontrib = delta_t * (ops_.div * (flux * face_b)
+                                                  + well_contrib);
                 const V qcontrib = delta_t * q;
                 const ADB pvcontrib = ADB::constant(pv*z0, bpat);
                 const ADB component_contrib = pvcontrib + qcontrib;
                 divcontrib_sum = divcontrib_sum - divcontrib/cell_b;
                 cell_residual_ = cell_residual_ - (component_contrib/cell_b);
 #else
-                const ADB component_contrib = pv*z0 + delta_t*(q - (ops_.div * (flux * face_b)));
+                const ADB component_contrib = pv*z0
+                    + delta_t*(q - (ops_.div * (flux * face_b) + well_contrib));
                 cell_residual_ = cell_residual_ - (component_contrib / cell_b);
 #endif
             }
@@ -283,6 +314,10 @@ namespace Opm {
 #endif
         }
 
+
+
+
+
         void
         solveJacobianSystem(BlackoilState& state) const
         {
@@ -306,12 +341,20 @@ namespace Opm {
             std::copy(&p[0], &p[0] + nc, state.pressure().begin());
         }
 
+
+
+
+
         double
         residualNorm() const
         {
             return cell_residual_.value().matrix().norm();
         }
 
+
+
+
+
         void
         computeFluxes(BlackoilState& state) const
         {
@@ -344,6 +387,9 @@ namespace Opm {
         }
 
 
+
+
+
         V fluidMu(const int phase, const V& p, const std::vector<int>& cells) const
         {
             switch (phase) {
@@ -359,6 +405,11 @@ namespace Opm {
                 THROW("Unknown phase index " << phase);
             }
         }
+
+
+
+
+
         ADB fluidMu(const int phase, const ADB& p, const std::vector<int>& cells) const
         {
             switch (phase) {
@@ -375,6 +426,10 @@ namespace Opm {
             }
         }
 
+
+
+
+
         ADB fluidFvf(const int phase, const ADB& p, const std::vector<int>& cells) const
         {
             switch (phase) {
@@ -391,6 +446,10 @@ namespace Opm {
             }
         }
 
+
+
+
+
         ADB fluidRho(const int phase, const ADB& p, const std::vector<int>& cells) const
         {
             const double* rhos = fluid_.surfaceDensity();
@@ -400,10 +459,23 @@ namespace Opm {
         }
 
 
+
+
+
         V fluidKr(const int phase) const
         {
             return kr_[phase];
         }
+
+
+
+
+
+        V fluidKrWell(const int phase) const
+        {
+            return well_kr_[phase];
+        }
+
     };
 } // namespace Opm