Add support for MISC regions

- element wise power operator impemented in AutoDiffBlock
- TL parameters are given pr cell

opm/autodiff/AutoDiffBlock.hpp

@@ -622,6 +622,36 @@ namespace Opm
         return rhs * lhs; // Commutative operation.
     }
 
+    /**
+     * @brief Computes the value of base raised to the power of exp elementwise
+     *
+     * @param base The AD forward block
+     * @param exp  array of exponent
+     * @return The value of base raised to the power of exp elementwise
+     */
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
+                              const typename AutoDiffBlock<Scalar>::V& exp)
+    {
+        const int num_elem = base.value().size();
+        typename AutoDiffBlock<Scalar>::V val (num_elem);
+        typename AutoDiffBlock<Scalar>::V derivative = exp;
+        assert(exp.size() == num_elem);
+        for (int i = 0; i < num_elem; ++i) {
+            val[i] = std::pow(base.value()[i], exp[i]);
+            derivative[i] *= std::pow(base.value()[i],exp[i]-1.0);
+        }
+        const typename AutoDiffBlock<Scalar>::M derivative_diag(derivative.matrix().asDiagonal());
+
+        std::vector< typename AutoDiffBlock<Scalar>::M > jac (base.numBlocks());
+        for (int block = 0; block < base.numBlocks(); block++) {
+             fastSparseProduct(derivative_diag, base.derivative()[block], jac[block]);
+        }
+
+        return AutoDiffBlock<Scalar>::function( std::move(val), std::move(jac) );
+    }
+
+
     /**
      * @brief Computes the value of base raised to the power of exp
      *

opm/autodiff/BlackoilSolventModel_impl.hpp

@@ -825,12 +825,12 @@ 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_.mixingParameterViscosity();
+        const V mix_param_mu = solvent_props_.mixingParameterViscosity(cells_);
 
         // Update viscosities
-        viscosity[pu.phase_pos[ Oil ]] = pow(mu_o,1.0 - mix_param_mu) * pow(mu_mos,mix_param_mu);
+        viscosity[pu.phase_pos[ Oil ]] = pow(mu_o,ones - mix_param_mu) * pow(mu_mos,mix_param_mu);
-        viscosity[pu.phase_pos[ Gas ]] = pow(mu_g,1.0 - mix_param_mu) * pow(mu_msg,mix_param_mu);
+        viscosity[pu.phase_pos[ Gas ]] = pow(mu_g,ones - mix_param_mu) * pow(mu_msg,mix_param_mu);
-        viscosity[solvent_pos_] = pow(mu_s,1.0 - mix_param_mu) * pow(mu_m,mix_param_mu);
+        viscosity[solvent_pos_] = pow(mu_s,ones - mix_param_mu) * pow(mu_m,mix_param_mu);
 
         // Density
         ADB& rho_o = density[pu.phase_pos[ Oil ]];
@@ -838,13 +838,13 @@ namespace Opm {
         ADB& rho_s = density[solvent_pos_];
 
         // mixing parameter for density
-        const double mix_param_rho = solvent_props_.mixingParameterDensity();
+        const V mix_param_rho = solvent_props_.mixingParameterDensity(cells_);
 
         // compute effective viscosities for density calculations. These have to
         // be recomputed as a different mixing parameter may be used.
-        const ADB mu_o_eff = pow(mu_o,1.0 - mix_param_rho) * pow(mu_mos,mix_param_rho);
+        const ADB mu_o_eff = pow(mu_o,ones - mix_param_rho) * pow(mu_mos,mix_param_rho);
-        const ADB mu_g_eff = pow(mu_g,1.0 - mix_param_rho) * pow(mu_msg,mix_param_rho);
+        const ADB mu_g_eff = pow(mu_g,ones - mix_param_rho) * pow(mu_msg,mix_param_rho);
-        const ADB mu_s_eff = pow(mu_s,1.0 - mix_param_rho) * pow(mu_m,mix_param_rho);
+        const ADB mu_s_eff = pow(mu_s,ones - mix_param_rho) * pow(mu_m,mix_param_rho);
 
         const ADB sog_eff = so_eff + sg_eff;
         const ADB sof = so_eff / sog_eff;

opm/autodiff/SolventPropsAdFromDeck.cpp

@@ -233,30 +233,26 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
 
             if (deck->hasKeyword("TLMIXPAR")) {
                 const int numRegions = deck->getKeyword("TLMIXPAR")->size();
-                if (numRegions > 1) {
-                    OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR.");
-                }
-                const auto tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(0);
-                std::vector<double> mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData();
 
-                mix_param_viscosity_ = mix_params_viscosity[0];
+                // resize the attributes of the object
-
+                mix_param_viscosity_.resize(numRegions);
-                std::vector<double> mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData();
+                mix_param_density_.resize(numRegions);
-                const int numDensityItems = mix_params_density.size();
+                for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
-                if (numDensityItems == 0) {
+                    const auto& tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(regionIdx);
-                    mix_param_density_ = mix_param_viscosity_;
+                    const auto& mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData();
-                } else if (numDensityItems == 1) {
+                    mix_param_viscosity_[regionIdx] = mix_params_viscosity[0];
-                    mix_param_density_ = mix_params_density[0];
+                    const auto& mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData();
-                } else {
+                    const int numDensityItems = mix_params_density.size();
-                    OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR.");
+                    if (numDensityItems == 0) {
+                        mix_param_density_[regionIdx] = mix_param_viscosity_[regionIdx];
+                    } else if (numDensityItems == 1) {
+                        mix_param_density_[regionIdx] = mix_params_density[0];
+                    } else {
+                        OPM_THROW(std::runtime_error, "Only one value can be entered for the TL parameter pr MISC region.");
+                    }
                 }
-            } else {
-                mix_param_viscosity_ = 0.0;
-                mix_param_density_ = 0.0;
             }
 
-
-
         }
     }
 
@@ -397,12 +393,32 @@ V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const {
     return density;
 }
 
-double SolventPropsAdFromDeck::mixingParameterViscosity() const {
+V SolventPropsAdFromDeck::mixingParameterViscosity(const Cells& cells) const {
-    return mix_param_viscosity_;
+    const int n = cells.size();
+    if (mix_param_viscosity_.size() > 0) {
+        V mix_param(n);
+        for (int i = 0; i < n; ++i) {
+            int regionIdx = cellMiscRegionIdx_[cells[i]];
+            mix_param[i] = mix_param_viscosity_[regionIdx];
+        }
+        return mix_param;
+    }
+    // return zeros if not specified
+    return V::Zero(n);
 }
 
-double SolventPropsAdFromDeck::mixingParameterDensity() const {
+V SolventPropsAdFromDeck::mixingParameterDensity(const Cells& cells) const {
-    return mix_param_density_;
+    const int n = cells.size();
+    if (mix_param_viscosity_.size() > 0) {
+        V mix_param(n);
+        for (int i = 0; i < n; ++i) {
+            int regionIdx = cellMiscRegionIdx_[cells[i]];
+            mix_param[i] = mix_param_density_[regionIdx];
+        }
+        return mix_param;
+    }
+    // return zeros if not specified
+    return V::Zero(n);
 }
 
 void SolventPropsAdFromDeck::extractTableIndex(const std::string& keyword,

opm/autodiff/SolventPropsAdFromDeck.hpp

@@ -128,10 +128,14 @@ public:
     V solventSurfaceDensity(const Cells& cells) const;
 
     /// Todd-Longstaff mixing parameter for viscosity calculation
-    double mixingParameterViscosity() const;
+    /// \param[in]  cells           Array of n cell indices to be associated with the fraction values.
+    /// return                      Array of n mixing paramters for viscosity calculation
+    V mixingParameterViscosity(const Cells& cells) const;
 
     /// Todd-Longstaff mixing parameter for density calculation
-    double mixingParameterDensity() const;
+    /// \param[in]  cells           Array of n cell indices to be associated with the fraction values.
+    /// return                      Array of n mixing paramters for density calculation
+    V mixingParameterDensity(const Cells& cells) const;
 
 
 private:
@@ -175,8 +179,8 @@ private:
     std::vector<NonuniformTableLinear<double> > misc_;
     std::vector<NonuniformTableLinear<double> > sorwmis_;
     std::vector<NonuniformTableLinear<double> > sgcwmis_;
-    double mix_param_viscosity_;
+    std::vector<double> mix_param_viscosity_;
-    double mix_param_density_;
+    std::vector<double> mix_param_density_;
 };
 
 } // namespace OPM