diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index d8c33322a..3f212b63c 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/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 + AutoDiffBlock pow(const AutoDiffBlock& base, + const typename AutoDiffBlock::V& exp) + { + const int num_elem = base.value().size(); + typename AutoDiffBlock::V val (num_elem); + typename AutoDiffBlock::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::M derivative_diag(derivative.matrix().asDiagonal()); + + std::vector< typename AutoDiffBlock::M > jac (base.numBlocks()); + for (int block = 0; block < base.numBlocks(); block++) { + fastSparseProduct(derivative_diag, base.derivative()[block], jac[block]); + } + + return AutoDiffBlock::function( std::move(val), std::move(jac) ); + } + + /** * @brief Computes the value of base raised to the power of exp * diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index 9a2c5fee5..a6eed3b13 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/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[ Gas ]] = pow(mu_g,1.0 - 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[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,ones - mix_param_mu) * pow(mu_msg,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_g_eff = pow(mu_g,1.0 - 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_o_eff = pow(mu_o,ones - mix_param_rho) * pow(mu_mos,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,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; diff --git a/opm/autodiff/SolventPropsAdFromDeck.cpp b/opm/autodiff/SolventPropsAdFromDeck.cpp index f4bd097e5..f4b1d31d6 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.cpp +++ b/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 mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData(); - mix_param_viscosity_ = mix_params_viscosity[0]; - - std::vector mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData(); - const int numDensityItems = mix_params_density.size(); - if (numDensityItems == 0) { - mix_param_density_ = mix_param_viscosity_; - } else if (numDensityItems == 1) { - mix_param_density_ = mix_params_density[0]; - } else { - OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR."); + // resize the attributes of the object + mix_param_viscosity_.resize(numRegions); + mix_param_density_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + const auto& tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(regionIdx); + const auto& mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData(); + mix_param_viscosity_[regionIdx] = mix_params_viscosity[0]; + const auto& mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData(); + const int numDensityItems = mix_params_density.size(); + 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 { - return mix_param_viscosity_; +V SolventPropsAdFromDeck::mixingParameterViscosity(const Cells& cells) const { + 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 { - return mix_param_density_; +V SolventPropsAdFromDeck::mixingParameterDensity(const Cells& cells) const { + 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, diff --git a/opm/autodiff/SolventPropsAdFromDeck.hpp b/opm/autodiff/SolventPropsAdFromDeck.hpp index c2e439900..ad03a80b3 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.hpp +++ b/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 > misc_; std::vector > sorwmis_; std::vector > sgcwmis_; - double mix_param_viscosity_; - double mix_param_density_; + std::vector mix_param_viscosity_; + std::vector mix_param_density_; }; } // namespace OPM