Cleaning up various issues.

ebos/eclfluxmodule.hh

@@ -272,7 +272,8 @@ protected:
 
             const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
             Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx));
-            if (enableExtbo) //z-dependency ...
+            if (enableExtbo) // added stability; particulary useful for solvent migrating in pure water
+                             // where the solvent fraction displays a 0/1 behaviour ...
                 pressureExterior += Toolbox::value(rhoAvg)*(distZ*g);
             else
                 pressureExterior += rhoAvg*(distZ*g);

ebos/ecloutputblackoilmodule.hh

@@ -637,7 +637,12 @@ public:
                 if (viscosity_[phaseIdx].size() == 0)
                     continue;
 
-                viscosity_[phaseIdx][globalDofIdx] = Opm::getValue(fs.viscosity(phaseIdx));
+                if (extboX_.size() > 0 && phaseIdx==oilPhaseIdx)
+                    viscosity_[phaseIdx][globalDofIdx] = Opm::getValue(intQuants.oilViscosity());
+                else if (extboX_.size() > 0 && phaseIdx==gasPhaseIdx)
+                    viscosity_[phaseIdx][globalDofIdx] = Opm::getValue(intQuants.gasViscosity());
+                else
+                    viscosity_[phaseIdx][globalDofIdx] = Opm::getValue(fs.viscosity(phaseIdx));
                 Opm::Valgrind::CheckDefined(viscosity_[phaseIdx][globalDofIdx]);
             }
 
@@ -666,11 +671,11 @@ public:
             }
 
             if (extboX_.size() > 0) {
-                extboX_[globalDofIdx] = intQuants.xvalue().value();
+                extboX_[globalDofIdx] = intQuants.xVolume().value();
             }
 
             if (extboY_.size() > 0) {
-                extboY_[globalDofIdx] = intQuants.yvalue().value();
+                extboY_[globalDofIdx] = intQuants.yVolume().value();
             }
 
             if (extboZ_.size() > 0) {
@@ -680,21 +685,19 @@ public:
             if (mFracCo2_.size() > 0) {
                 const Scalar stdVolOil = Opm::getValue(fs.saturation(oilPhaseIdx))*Opm::getValue(fs.invB(oilPhaseIdx))
                                        + Opm::getValue(fs.saturation(gasPhaseIdx))*Opm::getValue(fs.invB(gasPhaseIdx))*Opm::getValue(fs.Rv());
-                const Scalar stdVolGas = Opm::getValue(fs.saturation(gasPhaseIdx))*Opm::getValue(fs.invB(gasPhaseIdx))*(1.0-intQuants.yvalue().value())
+                const Scalar stdVolGas = Opm::getValue(fs.saturation(gasPhaseIdx))*Opm::getValue(fs.invB(gasPhaseIdx))*(1.0-intQuants.yVolume().value())
-                                       + Opm::getValue(fs.saturation(oilPhaseIdx))*Opm::getValue(fs.invB(oilPhaseIdx))*Opm::getValue(fs.Rs())*(1.0-intQuants.xvalue().value());
+                                       + Opm::getValue(fs.saturation(oilPhaseIdx))*Opm::getValue(fs.invB(oilPhaseIdx))*Opm::getValue(fs.Rs())*(1.0-intQuants.xVolume().value());
-                const Scalar stdVolCo2 = Opm::getValue(fs.saturation(gasPhaseIdx))*Opm::getValue(fs.invB(gasPhaseIdx))*intQuants.yvalue().value()
+                const Scalar stdVolCo2 = Opm::getValue(fs.saturation(gasPhaseIdx))*Opm::getValue(fs.invB(gasPhaseIdx))*intQuants.yVolume().value()
-                                       + Opm::getValue(fs.saturation(oilPhaseIdx))*Opm::getValue(fs.invB(oilPhaseIdx))*Opm::getValue(fs.Rs())*intQuants.xvalue().value();
+                                       + Opm::getValue(fs.saturation(oilPhaseIdx))*Opm::getValue(fs.invB(oilPhaseIdx))*Opm::getValue(fs.Rs())*intQuants.xVolume().value();
-                const Scalar rhoO= 836.19;
+                const Scalar rhoO= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
-                const Scalar rhoG= 1.2984;
+                const Scalar rhoG= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
-                const Scalar rhoCO2=1.8727048125;
+                const Scalar rhoCO2= intQuants.zRefDensity();
                 const Scalar stdMassTotal= 1.0e-10 + stdVolOil*rhoO + stdVolGas*rhoG + stdVolCo2*rhoCO2;
                 mFracOil_[globalDofIdx] = stdVolOil*rhoO/stdMassTotal;
                 mFracGas_[globalDofIdx] = stdVolGas*rhoG/stdMassTotal;
                 mFracCo2_[globalDofIdx] = stdVolCo2*rhoCO2/stdMassTotal;
             }
 
-
-
             if (bubblePointPressure_.size() > 0) {
                 try {
                     bubblePointPressure_[globalDofIdx] = Opm::getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()));
@@ -1084,22 +1087,22 @@ public:
             sol.insert ("SSOLVENT", Opm::UnitSystem::measure::identity, std::move(sSol_), Opm::data::TargetType::RESTART_SOLUTION);
 
         if (extboX_.size() > 0)
-            sol.insert ("SS_X", Opm::UnitSystem::measure::identity, std::move(extboX_), Opm::data::TargetType::RESTART_AUXILIARY);
+            sol.insert ("SS_X", Opm::UnitSystem::measure::identity, std::move(extboX_), Opm::data::TargetType::RESTART_SOLUTION);
 
         if (extboY_.size() > 0)
-            sol.insert ("SS_Y", Opm::UnitSystem::measure::identity, std::move(extboY_), Opm::data::TargetType::RESTART_AUXILIARY);
+            sol.insert ("SS_Y", Opm::UnitSystem::measure::identity, std::move(extboY_), Opm::data::TargetType::RESTART_SOLUTION);
 
         if (extboZ_.size() > 0)
-            sol.insert ("SS_Z", Opm::UnitSystem::measure::identity, std::move(extboZ_), Opm::data::TargetType::RESTART_AUXILIARY);
+            sol.insert ("SS_Z", Opm::UnitSystem::measure::identity, std::move(extboZ_), Opm::data::TargetType::RESTART_SOLUTION);
 
         if (mFracOil_.size() > 0)
-            sol.insert ("STD_OIL", Opm::UnitSystem::measure::identity, std::move(mFracOil_), Opm::data::TargetType::RESTART_AUXILIARY);
+            sol.insert ("STD_OIL", Opm::UnitSystem::measure::identity, std::move(mFracOil_), Opm::data::TargetType::RESTART_SOLUTION);
 
         if (mFracGas_.size() > 0)
-            sol.insert ("STD_GAS", Opm::UnitSystem::measure::identity, std::move(mFracGas_), Opm::data::TargetType::RESTART_AUXILIARY);
+            sol.insert ("STD_GAS", Opm::UnitSystem::measure::identity, std::move(mFracGas_), Opm::data::TargetType::RESTART_SOLUTION);
 
         if (mFracCo2_.size() > 0)
-            sol.insert ("STD_CO2", Opm::UnitSystem::measure::identity, std::move(mFracCo2_), Opm::data::TargetType::RESTART_AUXILIARY);
+            sol.insert ("STD_CO2", Opm::UnitSystem::measure::identity, std::move(mFracCo2_), Opm::data::TargetType::RESTART_SOLUTION);
 
         if (cPolymer_.size() > 0)
             sol.insert ("POLYMER", Opm::UnitSystem::measure::identity, std::move(cPolymer_), Opm::data::TargetType::RESTART_SOLUTION);

opm/simulators/wells/StandardWell_impl.hpp

@@ -826,9 +826,9 @@ namespace Opm
                 cq_s_zfrac *= wsolvent();
             } else if (cq_s_zfrac.value() != 0.0) {
                 const double dis_gas_frac = perf_dis_gas_rate / cq_s_zfrac.value();
-                cq_s_zfrac *= extendEval(dis_gas_frac*intQuants.xvalue() + (1.0-dis_gas_frac)*intQuants.yvalue());
+                cq_s_zfrac *= extendEval(dis_gas_frac*intQuants.xVolume() + (1.0-dis_gas_frac)*intQuants.yVolume());
             }
-            well_state.perfRateZFrac()[first_perf_ + perf] = cq_s_zfrac.value();
+            well_state.perfRateSolvent()[first_perf_ + perf] = cq_s_zfrac.value();
 
             cq_s_zfrac *= well_efficiency_factor_;
             connectionRates_[perf][contiZfracEqIdx] = Base::restrictEval(cq_s_zfrac);
@@ -2004,7 +2004,7 @@ namespace Opm
                     const double oilrate = std::abs(well_state.wellRates()[oilpos_well]); //in order to handle negative rates in producers
                     rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
                     if (oilrate > 0) {
-                        const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - well_state.solventWellRate(w);
+                        const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
                         double rv = 0.0;
                         if (gasrate > 0) {
                             rv = oilrate / gasrate;
@@ -2029,7 +2029,7 @@ namespace Opm
                 if (gasPresent) {
                     rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
                     const int gaspos_well = pu.phase_pos[Gas] + w * pu.num_phases;
-                    const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - well_state.solventWellRate(w);
+                    const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
                     if (gasrate > 0) {
                         const double oilrate = std::abs(well_state.wellRates()[oilpos_well]);
                         double rs = 0.0;
@@ -2926,7 +2926,8 @@ namespace Opm
                 primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate;
             }
             if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / total_well_rate ;
+                primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]]
+                                                 - (has_solvent ? well_state.solventWellRate(well_index) : 0.0) ) / total_well_rate ;
             }
             if (has_solvent) {
                 primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ;

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -163,10 +163,6 @@ namespace Opm
 
             perfRateBrine_.clear();
             perfRateBrine_.resize(nperf, 0.0);
-
-            perfRateZFrac_.clear();
-            perfRateZFrac_.resize(nperf, 0.0);
-
             // intialize wells that have been there before
             // order may change so the mapping is based on the well name
             if (prevState && !prevState->wellMap().empty()) {
@@ -579,7 +575,7 @@ namespace Opm
                     well.rates.set( rt::well_potential_gas, this->well_potentials_[well_rate_index + pu.phase_pos[Gas]] );
                 }
 
-                if ( pu.has_solvent ) {
+                if ( pu.has_solvent || pu.has_zFraction) {
                     well.rates.set( rt::solvent, solventWellRate(w) );
                 }
 
@@ -598,10 +594,6 @@ namespace Opm
                     well.rates.set( rt::alq, 0.0 );
                 }
 
-                if ( pu.has_zFraction ) { // shared with the standard solvent model ...
-                    well.rates.set( rt::solvent, zFracWellRate(w) );
-                }
-
                 well.rates.set( rt::dissolved_gas, this->well_dissolved_gas_rates_[w] );
                 well.rates.set( rt::vaporized_oil, this->well_vaporized_oil_rates_[w] );
 
@@ -628,11 +620,8 @@ namespace Opm
                     if ( pu.has_brine ) {
                         comp.rates.set( rt::brine, this->perfRateBrine()[wt.second[1] + local_comp_index]);
                     }
-                    if ( pu.has_zFraction ) { // shared with the standard solvent model ...
+                    if ( pu.has_solvent || pu.has_zFraction) {
-                        comp.rates.set( rt::solvent, this->perfRateZFrac()[local_comp_index]);
+                        comp.rates.set( rt::solvent, this->perfRateSolvent()[local_comp_index]);
-                    }
-                    if ( pu.has_solvent ) {
-                        comp.rates.set( rt::solvent, this->perfRateSolvent()[wt.second[1] + local_comp_index]);
                     }
 
                     ++local_comp_index;
@@ -886,15 +875,6 @@ namespace Opm
             return std::accumulate(&perfRateBrine_[0] + first_perf_index_[w], &perfRateBrine_[0] + first_perf_index_[w+1], 0.0);
         }
 
-        /// One rate pr well connection.
-        std::vector<double>& perfRateZFrac() { return perfRateZFrac_; }
-        const std::vector<double>& perfRateZFrac() const { return perfRateZFrac_; }
-
-        /// One rate pr well
-        double zFracWellRate(const int w) const {
-            return std::accumulate(&perfRateZFrac_[0] + first_perf_index_[w], &perfRateZFrac_[0] + first_perf_index_[w+1], 0.0);
-        }
-
         std::vector<double>& wellReservoirRates()
         {
             return well_reservoir_rates_;
@@ -1226,7 +1206,6 @@ namespace Opm
         // only for output
         std::vector<double> perfRatePolymer_;
         std::vector<double> perfRateBrine_;
-        std::vector<double> perfRateZFrac_;
 
         // it is the throughput of water flow through the perforations
         // it is used as a measure of formation damage around well-bore due to particle deposition