Cleaning and adding comments

opm/autodiff/BlackoilSolventModel.hpp

@@ -147,9 +147,6 @@ namespace Opm {
         using Base::computePressures;
         using Base::computeGasPressure;
         using Base::applyThresholdPressures;
-        //using Base::fluidViscosity;
-        //using Base::fluidReciprocFVF;
-        //using Base::fluidDensity;
         using Base::fluidRsSat;
         using Base::fluidRvSat;
         using Base::poroMult;
@@ -168,9 +165,6 @@ namespace Opm {
         std::vector<ADB>
         computeRelPerm(const SolutionState& state) const;
 
-        void calculateEffectiveProperties(const SolutionState&  state);
-
-
         ADB
         fluidViscosity(const int               phase,
                        const ADB&              p    ,
@@ -237,10 +231,12 @@ namespace Opm {
         const std::vector<PhasePresence>
         phaseCondition() const {return this->phaseCondition_;}
 
+        // compute effective viscosities (mu_eff_) and effective b factors (b_eff_)  using the ToddLongstaff model
+        void calculateEffectiveProperties(const SolutionState&  state);
+
+        // compute density and viscosity using the ToddLongstaff mixing model
         void ToddLongstaffModel(std::vector<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const Opm::PhaseUsage pu);
 
-
-
     };
 
 

opm/autodiff/BlackoilSolventModel_impl.hpp

@@ -387,12 +387,16 @@ namespace Opm {
         }
 
         if (pu.phase_used[BlackoilPhases::Vapour]) {
+            // Unclear wether the effective or the pure values should be used for the wells
+            // the current usage of unmodified properties values gives best match.
             //V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value();
             V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
             V rhog = fluid_.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
             if (has_solvent_) {
+
                 const V bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value();
                 //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value();
+
                 // A weighted sum of the b-factors of gas and solvent are used.
                 const int nc = Opm::AutoDiffGrid::numCells(grid_);
 
@@ -529,49 +533,38 @@ namespace Opm {
                                                 const SolutionState&    state)
     {
 
+        const int canonicalPhaseIdx = canph_[ actph ];
         // make a copy to make it possible to modify it
         ADB kr_mod = kr;
+        if (canonicalPhaseIdx == Gas) {
             if (has_solvent_) {
                 const int  nc   = Opm::UgGridHelpers::numCells(grid_);
                 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
                 const ADB zero = ADB::constant(V::Zero(nc));
                 const V ones = V::Constant(nc, 1.0);
-            const int canonicalPhaseIdx = canph_[ actph ];
 
                 const ADB& ss = state.solvent_saturation;
                 const ADB& sg = (active_[ Gas ]
                                  ? state.saturation[ pu.phase_pos[ Gas ] ]
                                  : zero);
-
                 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
                 const ADB F_solvent = zero_selector.select(zero, ss / (ss + sg));
 
+                // Compute solvent properties
                 const std::vector<PhasePresence>& cond = phaseCondition();
                 ADB mu_s = fluidViscosity(Solvent, phasePressure,state.temperature, state.rs, state.rv, cond);
                 ADB rho_s = fluidDensity(Solvent,rq_[solvent_pos_].b, state.rs, state.rv);
 
-            if (canonicalPhaseIdx == Gas) {
+                // Compute solvent relperm and mass flux
-
-                // compute solvent mobility and flux
-                //const ADB tr_mult = transMult(state.pressure);
-
                 ADB krs = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * kr_mod;
                 Base::computeMassFlux(solvent_pos_, transi, krs, mu_s, rho_s, phasePressure, state);
 
-                //rq_[solvent_pos_].mob = krs * tr_mult / mu_s;
+                // Modify gas relperm
-
-                //const ADB rhoavg_solvent = ops_.caver * rho_s;
-                //rq_[ solvent_pos_ ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg_solvent * (ops_.ngrad * geo_.z().matrix()));
-
-                //UpwindSelector<double> upwind_solvent(grid_, ops_, rq_[solvent_pos_].dh.value());
-                //rq_[solvent_pos_].mflux = upwind_solvent.select(rq_[solvent_pos_].b * rq_[solvent_pos_].mob) * (transi * rq_[solvent_pos_].dh);
-
-                // modify gas relperm
                 kr_mod = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * kr_mod;
 
             }
         }
-        // compute mobility and flux
+        // Compute mobility and flux
         Base::computeMassFlux(actph, transi, kr_mod, mu, rho, phasePressure, state);
 
     }
@@ -832,7 +825,7 @@ namespace Opm {
         const ADB mu_m = zero_selectorSn.select(mu_s + mu_o + mu_g, mu_o * mu_s * mu_g / pow( ( (so_eff / sn_eff) * mu_s_pow *  mu_g_pow)
                                                        + ( (ss_eff / sn_eff) * mu_o_pow *  mu_g_pow) + ( (sg_eff / sn_eff) * mu_s_pow * mu_o_pow), 4.0));
         // Mixing parameter for viscosity
-        const double mix_param_mu = solvent_props_.mixingParamterViscosity();
+        const double mix_param_mu = solvent_props_.mixingParameterViscosity();
 
         // Update viscosities
         viscosity[pu.phase_pos[ Oil ]] = pow(mu_o,1.0 - mix_param_mu) * pow(mu_mos,mix_param_mu);
@@ -844,8 +837,8 @@ namespace Opm {
         ADB& rho_g = density[pu.phase_pos[ Gas ]];
         ADB& rho_s = density[solvent_pos_];
 
-        // mixing paramter for density
+        // mixing parameter for density
-        const double mix_param_rho = solvent_props_.mixingParamterDensity();
+        const double mix_param_rho = solvent_props_.mixingParameterDensity();
 
         // compute effective viscosities for density calculations. These have to
         // be recomputed as a different mixing parameter may be used.

opm/autodiff/SolventPropsAdFromDeck.cpp

@@ -142,11 +142,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
                     const auto& sn = sof2Table.getSoColumn();
                     const auto& krn = sof2Table.getKroColumn();
 
-                    for (size_t i = 0; i < sn.size(); ++i) {
-                        std::cout << sn[i] << " " << krn[i] <<std::endl;
-                    }
-
-
                     krn_[regionIdx] = NonuniformTableLinear<double>(sn, krn);
                 }
 
@@ -307,35 +302,35 @@ ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg,
 ADB SolventPropsAdFromDeck::bSolvent(const ADB& pg,
                                 const Cells& cells) const
 {
-    return SolventPropsAdFromDeck::makeAD(pg, cells, b_);
+    return SolventPropsAdFromDeck::makeADBfromTables(pg, cells, b_);
 
 }
 
 ADB SolventPropsAdFromDeck::gasRelPermMultiplier(const ADB& solventFraction,
                                  const Cells& cells) const
 {
-    return SolventPropsAdFromDeck::makeAD(solventFraction, cells, krg_);
+    return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krg_);
 
 }
 
 ADB SolventPropsAdFromDeck::solventRelPermMultiplier(const ADB& solventFraction,
                                  const Cells& cells) const
 {
-    return SolventPropsAdFromDeck::makeAD(solventFraction, cells, krs_);
+    return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krs_);
 }
 
 
 ADB SolventPropsAdFromDeck::misicibleHydrocarbonWaterRelPerm(const ADB& Sn,
                                  const Cells& cells) const
 {
-    return SolventPropsAdFromDeck::makeAD(Sn, cells, krn_);
+    return SolventPropsAdFromDeck::makeADBfromTables(Sn, cells, krn_);
 }
 
 ADB SolventPropsAdFromDeck::miscibleSolventGasRelPermMultiplier(const ADB& Ssg,
                                  const Cells& cells) const
 {
     if (mkrsg_.size() > 0) {
-        return SolventPropsAdFromDeck::makeAD(Ssg, cells, mkrsg_);
+        return SolventPropsAdFromDeck::makeADBfromTables(Ssg, cells, mkrsg_);
     }
     // trivial function if not specified
     return Ssg;
@@ -345,7 +340,7 @@ ADB SolventPropsAdFromDeck::miscibleOilRelPermMultiplier(const ADB& So,
                                  const Cells& cells) const
 {
     if (mkro_.size() > 0) {
-        return SolventPropsAdFromDeck::makeAD(So, cells, mkro_);
+        return SolventPropsAdFromDeck::makeADBfromTables(So, cells, mkro_);
     }
     // trivial function if not specified
     return So;
@@ -355,14 +350,14 @@ ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction,
                                  const Cells& cells) const
 {
 
-    return SolventPropsAdFromDeck::makeAD(solventFraction, cells, misc_);
+    return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, misc_);
 }
 
 
 ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw,
                                            const Cells& cells) const {
     if (sgcwmis_.size()>0) {
-        return SolventPropsAdFromDeck::makeAD(Sw, cells, sgcwmis_);
+        return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sgcwmis_);
     }
     // return zeros if not specified
     return ADB::constant(V::Zero(Sw.size()));
@@ -372,15 +367,15 @@ ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw
 ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw,
                                            const Cells& cells) const {
     if (sorwmis_.size()>0) {
-        return SolventPropsAdFromDeck::makeAD(Sw, cells, sorwmis_);
+        return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sorwmis_);
     }
     // return zeros if not specified
     return ADB::constant(V::Zero(Sw.size()));
 }
 
-ADB SolventPropsAdFromDeck::makeAD(const ADB& X_AD, const Cells& cells, std::vector<NonuniformTableLinear<double>> table) const {
+ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD, const Cells& cells, std::vector<NonuniformTableLinear<double>> table) const {
     const int n = cells.size();
-    assert(Sn.value().size() == n);
+    assert(X_AD.value().size() == n);
     V x(n);
     V dx(n);
     for (int i = 0; i < n; ++i) {
@@ -411,11 +406,11 @@ V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const {
     return density;
 }
 
-double SolventPropsAdFromDeck::mixingParamterViscosity() const {
+double SolventPropsAdFromDeck::mixingParameterViscosity() const {
     return mix_param_viscosity_;
 }
 
-double SolventPropsAdFromDeck::mixingParamterDensity() const {
+double SolventPropsAdFromDeck::mixingParameterDensity() const {
     return mix_param_density_;
 }
 

opm/autodiff/SolventPropsAdFromDeck.hpp

@@ -127,16 +127,17 @@ public:
     /// \return                     Array of n solvent density values.
     V solventSurfaceDensity(const Cells& cells) const;
 
-    /// Todd-Longstaff mixing paramter for viscosity calculation
+    /// Todd-Longstaff mixing parameter for viscosity calculation
-    double mixingParamterViscosity() const;
+    double mixingParameterViscosity() const;
 
-    /// Todd-Longstaff mixing paramter for density calculation
+    /// Todd-Longstaff mixing parameter for density calculation
-    double mixingParamterDensity() const;
+    double mixingParameterDensity() const;
 
 
 private:
 
-    ADB makeAD(const ADB& X, const Cells& cells,     std::vector<NonuniformTableLinear<double> > table) const;
+    /// Makes ADB from table values
+    ADB makeADBfromTables(const ADB& X, const Cells& cells,     std::vector<NonuniformTableLinear<double> > table) const;
 
     // The PVT region which is to be used for each cell
     std::vector<int> cellPvtRegionIdx_;

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -427,7 +427,6 @@ namespace Opm {
         const int canonicalPhaseIdx = canph_[ actph ];
         if (canonicalPhaseIdx == Water) {
             if (has_polymer_) {
-                const std::vector<PhasePresence>& cond = phaseCondition();
                 const ADB tr_mult = transMult(state.pressure);
                 const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
                 const ADB mc = computeMc(state);