Checkpoint

opm/autodiff/AutoDiffBlock.hpp

@@ -292,6 +292,14 @@ namespace Opm
         /// Elementwise operator -
         AutoDiffBlock operator-(const AutoDiffBlock& rhs) const
         {
+            assert(val_.rows() > 0 || rhs.val_.rows() > 0);
+            if (val_.rows() == 0) {
+                return rhs*(-1.0);
+            }
+            if (rhs.val_.rows() == 0) {
+                return *this;
+            }
+
             if (jac_.empty() && rhs.jac_.empty()) {
                 return constant(val_ - rhs.val_);
             }

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -658,12 +658,16 @@ namespace detail {
                 sg = isSg_*xvar + isRv_*so;
                 so -= sg;
 
-                if (active_[ Oil ]) {
+                //Compute the phase pressures before computing RS/RV
-                    // RS and RV is only defined if both oil and gas phase are active.
+                {
                     const ADB& sw = (active_[ Water ]
                                              ? state.saturation[ pu.phase_pos[ Water ] ]
-                                             : ADB::constant(V::Zero(nc, 1)));
+                                             : ADB::constant(V::Zero(nc, 1), sg.blockPattern()));
                     state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
+                }
+
+                if (active_[ Oil ]) {
+                    // RS and RV is only defined if both oil and gas phase are active.
                     const ADB rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_);
                     if (has_disgas_) {
                         state.rs = (1-isRs_)*rsSat + isRs_*xvar;
@@ -678,6 +682,14 @@ namespace detail {
                     }
                 }
             }
+            else {
+                // Compute phase pressures also if gas phase is not active
+                const ADB& sw = (active_[ Water ]
+                                         ? state.saturation[ pu.phase_pos[ Water ] ]
+                                         : ADB::constant(V::Zero(nc, 1), so.blockPattern()));
+                const ADB& sg = ADB::constant(V::Zero(nc, 1), so.blockPattern());
+                state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
+            }
 
             if (active_[ Oil ]) {
                 // Note that so is never a primary variable.

opm/autodiff/BlackoilPropsAdFromDeck.cpp

@@ -373,7 +373,9 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
                 muLad = oilPvt_->saturatedViscosity(pvtRegionIdx, TLad, pLad);
             }
             else {
-                RsLad.value = rs.value()[i];
+                if (phase_usage_.phase_used[Gas]) {
+                    RsLad.value = rs.value()[i];
+                }
                 muLad = oilPvt_->viscosity(pvtRegionIdx, TLad, pLad, RsLad);
             }
 
@@ -388,9 +390,11 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
         std::vector<ADB::M> jacs(num_blocks);
         for (int block = 0; block < num_blocks; ++block) {
             fastSparseProduct(dmudp_diag, po.derivative()[block], jacs[block]);
-            ADB::M temp;
+            if (phase_usage_.phase_used[Gas]) {
-            fastSparseProduct(dmudr_diag, rs.derivative()[block], temp);
+                ADB::M temp;
-            jacs[block] += temp;
+                fastSparseProduct(dmudr_diag, rs.derivative()[block], temp);
+                jacs[block] += temp;
+            }
         }
         return ADB::function(std::move(mu), std::move(jacs));
     }
@@ -532,6 +536,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
 
         typedef Opm::LocalAd::Evaluation<double, /*size=*/2> LadEval;
 
+        // FIXME: What is Lad?
         LadEval pLad = 0.0;
         LadEval TLad = 0.0;
         LadEval RsLad = 0.0;
@@ -545,12 +550,18 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
             pLad.value = po.value()[i];
             TLad.value = T.value()[i];
 
-            if (cond[i].hasFreeGas()) {
+            if (phase_usage_.phase_used[Gas]) {
-                bLad = oilPvt_->saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
+                //RS/RV only makes sense when gas phase is active
+                if (cond[i].hasFreeGas()) {
+                    bLad = oilPvt_->saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
+                }
+                else {
+                    RsLad.value = rs.value()[i];
+                    bLad = oilPvt_->inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad);
+                }
             }
             else {
-                RsLad.value = rs.value()[i];
+                bLad = 0.0;
-                bLad = oilPvt_->inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad);
             }
 
             b[i] = bLad.value;
@@ -562,11 +573,14 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
         ADB::M dbdr_diag(dbdr.matrix().asDiagonal());
         const int num_blocks = po.numBlocks();
         std::vector<ADB::M> jacs(num_blocks);
+
         for (int block = 0; block < num_blocks; ++block) {
             fastSparseProduct(dbdp_diag, po.derivative()[block], jacs[block]);
-            ADB::M temp;
+            if (phase_usage_.phase_used[Gas]) {
-            fastSparseProduct(dbdr_diag, rs.derivative()[block], temp);
+                ADB::M temp;
-            jacs[block] += temp;
+                fastSparseProduct(dbdr_diag, rs.derivative()[block], temp);
+                jacs[block] += temp;
+            }
         }
         return ADB::function(std::move(b), std::move(jacs));
     }
@@ -853,6 +867,9 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
         std::vector<ADB> adbCapPressures;
         adbCapPressures.reserve(3);
         const ADB* s[3] = { &sw, &so, &sg };
+        //const std::vector<int>& block_pattern = sw.blockPattern();
+        assert(sw.blockPattern() == so.blockPattern());
+        assert(sw.blockPattern() == sg.blockPattern());
         for (int phase1 = 0; phase1 < 3; ++phase1) {
             if (phase_usage_.phase_used[phase1]) {
                 const int phase1_pos = phase_usage_.phase_pos[phase1];

opm/autodiff/StandardWells_impl.hpp

@@ -219,26 +219,22 @@ namespace Opm
         }
         assert(active[Oil]);
         const Vector perf_so =  subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
-        if (pu.phase_used[BlackoilPhases::Liquid]) {
+        if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
             const ADB perf_rs = subset(state.rs, well_cells);
             const Vector bo = fluid.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
             b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
             // const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
             const Vector rssat = fluid.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
             rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
-        } else {
+
-            rsmax_perf.assign(nperf, 0.0);
-        }
-        if (pu.phase_used[BlackoilPhases::Vapour]) {
             const ADB perf_rv = subset(state.rv, well_cells);
             const Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
             b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
             // const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
             const Vector rvsat = fluid.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
             rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
-        } else {
-            rvmax_perf.assign(nperf, 0.0);
         }
+
         // b is row major, so can just copy data.
         b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
 
@@ -347,8 +343,6 @@ namespace Opm
         const Vector& cdp = wellPerforationPressureDiffs();
         // Extract needed quantities for the perforation cells
         const ADB& p_perfcells = subset(state.pressure, well_cells);
-        const ADB& rv_perfcells = subset(state.rv, well_cells);
-        const ADB& rs_perfcells = subset(state.rs, well_cells);
 
         // Perforation pressure
         const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp;
@@ -401,6 +395,8 @@ namespace Opm
             const int gaspos = pu.phase_pos[Gas];
             const ADB cq_psOil = cq_ps[oilpos];
             const ADB cq_psGas = cq_ps[gaspos];
+            const ADB& rv_perfcells = subset(state.rv, well_cells);
+            const ADB& rs_perfcells = subset(state.rs, well_cells);
             cq_ps[gaspos] += rs_perfcells * cq_psOil;
             cq_ps[oilpos] += rv_perfcells * cq_psGas;
         }
@@ -440,21 +436,46 @@ namespace Opm
 
         // compute volume ratio 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]) {
+        if (active[Water]) {
-                const int gaspos = pu.phase_pos[Gas];
+            const int watpos = pu.phase_pos[Water];
-                tmp -= rv_perfcells * cmix_s[gaspos] / d;
+            volumeRatio += cmix_s[watpos] / b_perfcells[watpos];
-            }
-            if (phase == Gas && active[Oil]) {
-                const int oilpos = pu.phase_pos[Oil];
-                tmp -= rs_perfcells * cmix_s[oilpos] / d;
-            }
-            volumeRatio += tmp / b_perfcells[phase];
         }
 
+        if (active[Oil] && active[Gas]) {
+            // Incorporate RS/RV factors if both oil and gas active
+            const ADB& rv_perfcells = subset(state.rv, well_cells);
+            const ADB& rs_perfcells = subset(state.rs, well_cells);
+            const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
+
+            const int oilpos = pu.phase_pos[Oil];
+            const int gaspos = pu.phase_pos[Gas];
+
+
+            // First compute for Oil
+            {
+                const ADB tmp = cmix_s[oilpos] - (rv_perfcells * cmix_s[gaspos] / d);
+                volumeRatio += tmp / b_perfcells[oilpos];
+            }
+
+            // Then compute for Gas
+            {
+                const ADB tmp = cmix_s[gaspos] - (rs_perfcells * cmix_s[oilpos] / d);
+                volumeRatio += tmp / b_perfcells[gaspos];
+            }
+        }
+        else {
+            if (active[Oil]) {
+                const int oilpos = pu.phase_pos[Oil];
+                volumeRatio += cmix_s[oilpos] / b_perfcells[oilpos];
+            }
+            if (active[Gas]) {
+                const int gaspos = pu.phase_pos[Gas];
+                volumeRatio += cmix_s[gaspos] / b_perfcells[gaspos];
+            }
+        }
+
+
         // injecting connections total volumerates at standard conditions
         ADB cqt_is = cqt_i/volumeRatio;