add capillary pressure for incom solver.

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp

@@ -831,7 +831,7 @@ namespace {
 
         // convert the pressure offsets to the capillary pressures
         std::vector<ADB> pressure = fluid_.capPress(sw, so, sg, cells_);
-		pressure[0] = pressure[0] - pressure[0];
+		pressure[0] = pressure[0] - pressure[1];
 
         // add the total pressure to the capillary pressures
         for (int phaseIdx = 0; phaseIdx < 2; ++phaseIdx) {

opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp

@@ -548,6 +548,25 @@ namespace {
 
     }
    
+    std::vector<ADB>
+    FullyImplicitTwophasePolymerSolver::
+	computePressures(const SolutionState& state) const
+    {
+        const ADB sw = state.saturation[0];
+        const ADB so = state.saturation[1];
+
+        // convert the pressure offsets to the capillary pressures
+        std::vector<ADB> pressure = fluid_.capPress(sw, so, cells_);
+		pressure[0] = pressure[0] - pressure[1];
+
+        // add the total pressure to the capillary pressures
+        for (int phaseIdx = 0; phaseIdx < 2; ++phaseIdx) {
+            pressure[phaseIdx] += state.pressure;
+        }
+
+        return pressure;
+    }
+
     void
     FullyImplicitTwophasePolymerSolver::computeMassFlux(const V&                trans,
                                                         const ADB&              mc,
@@ -561,18 +580,18 @@ namespace {
         rq_[1].mob = kro / V::Constant(kro.size(), 1, mus[1]);
         rq_[2].mob =  mc * krw_eff * inv_wat_eff_vis;
 
-
         const int nc = grid_.number_of_cells; 
         V z(nc);
         // Compute z coordinates
         for (int c = 0; c < nc; ++c){
             z[c] = grid_.cell_centroids[c * 3 + 2];
         }
+		std::vector<ADB> press = computePressures(state);
         for (int phase = 0; phase < 2; ++phase) {
             const ADB rho = fluidDensity(phase, state.pressure);
 			ADB& head = rq_[phase].head;
             const ADB rhoavg = ops_.caver * rho;
-            const ADB dp = ops_.ngrad * state.pressure
+            const ADB dp = ops_.ngrad * press[phase]
                            - gravity_[2] * (rhoavg * (ops_.ngrad * z.matrix()));
             head = trans * dp;
             UpwindSelector<double> upwind(grid_, ops_, head.value());

opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp

@@ -102,6 +102,8 @@ namespace Opm {
         V
         transmissibility() const;
         
+		std::vector<ADB>
+		computePressures(const SolutionState& state) const;
 		void
         computeMassFlux(const V&                trans,
                         const ADB&              mc,

opm/polymer/fullyimplicit/IncompPropsAdFromDeck.cpp

@@ -132,5 +132,47 @@ namespace Opm
         }
         return relperms;
     }
+    std::vector<ADB> 
+	IncompPropsAdFromDeck::capPress(const ADB& sw,
+                                    const ADB& so,
+                                    const Cells& cells) const
+    {
+
+        const int n = cells.size();
+        const int np = numPhases();
+        const int num_blocks = so.numBlocks();
+        Block s_all(n, np);
+        assert(sw.size() == n && so.size() == n);
+        s_all.col(0) = sw.value();
+        s_all.col(1) = so.value();
+
+        Block pc(n, np);
+        Block dpc(n, np * np);
+
+        satprops_.capPress(n, s_all.data(), cells.data(), pc.data(), dpc.data());
+
+        std::vector<ADB> capPressures;
+        capPressures.reserve(2);
+        const ADB* s[2] = { &sw, &so};
+        for (int phase1 = 0; phase1 < 3; ++phase1) {
+        	const int phase1_pos = phase1;
+            std::vector<ADB::M> jacs(num_blocks);
+            for (int block = 0; block < num_blocks; ++block) {
+                jacs[block] = ADB::M(n, s[phase1]->derivative()[block].cols());
+            }
+            for (int phase2 = 0; phase2 < 3; ++phase2) {
+                const int phase2_pos = phase2;
+                // Assemble dpc1/ds2.
+                const int column = phase1_pos + phase2_pos; // Recall: Fortran ordering from props_.relperm()
+                ADB::M dpc1_ds2_diag = spdiag(dpc.col(column));
+                for (int block = 0; block < num_blocks; ++block) {
+                    jacs[block] += dpc1_ds2_diag * s[phase2]->derivative()[block];
+                }
+         	}
+            capPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs));
+		}
+        return capPressures;
+    }
+
 } //namespace Opm
 

opm/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp

@@ -46,6 +46,17 @@ namespace Opm
                                  const ADB& so,
                                  const Cells& cells) const;
 
+        /// Capillary pressure for all phases.
+        /// \param[in]  sw     Array of n water saturation values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the saturation values.
+        /// \return            An std::vector with 3 elements, each an array of n capillary pressure values,
+        ///                    containing the offsets for each p_g, p_o, p_w. The capillary pressure between
+        ///                    two arbitrary phases alpha and beta is then given as p_alpha - p_beta.
+        std::vector<ADB> capPress(const ADB& sw,
+                                  const ADB& so,
+                                  const Cells& cells) const;
+
     private:
         RockFromDeck rock_;
         PvtPropertiesIncompFromDeck pvt_;

opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp

@@ -82,6 +82,18 @@ namespace Opm
         std::vector<ADB> relperm(const ADB& sw,
                                  const ADB& so,
                                  const Cells& cells) const = 0;
+        /// Capillary pressure for all phases.
+        /// \param[in]  sw     Array of n water saturation values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the saturation values.
+        /// \return            An std::vector with 3 elements, each an array of n capillary pressure values,
+        ///                    containing the offsets for each p_g, p_o, p_w. The capillary pressure between
+        ///                    two arbitrary phases alpha and beta is then given as p_alpha - p_beta.
+        virtual
+        std::vector<ADB> capPress(const ADB& sw,
+                                  const ADB& so,
+                                  const Cells& cells) const = 0;
+
 
     };
 } // namespace Opm