diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index 0b6cb538e..66901ea4e 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -623,6 +623,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 exponents + * @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 * @@ -635,7 +665,7 @@ namespace Opm const double exp) { const typename AutoDiffBlock::V val = base.value().pow(exp); - const typename AutoDiffBlock::V derivative = exp * base.value().pow(exp-1.0); + const typename AutoDiffBlock::V derivative = exp * base.value().pow(exp - 1.0); const typename AutoDiffBlock::M derivative_diag(derivative.matrix().asDiagonal()); std::vector< typename AutoDiffBlock::M > jac (base.numBlocks()); diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index 9a2c5fee5..f7fda17b0 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -816,21 +816,21 @@ namespace Opm { Selector zero_selectorSsg(ssg_eff.value(), Selector::Zero); Selector zero_selectorSn(sn_eff.value(), Selector::Zero); - const ADB mu_s_pow = pow(mu_s,0.25); - const ADB mu_o_pow = pow(mu_o,0.25); - const ADB mu_g_pow = pow(mu_g,0.25); + const ADB mu_s_pow = pow(mu_s, 0.25); + const ADB mu_o_pow = pow(mu_o, 0.25); + const ADB mu_g_pow = pow(mu_g, 0.25); const ADB mu_mos = zero_selectorSos.select(mu_o + mu_s, mu_o * mu_s / pow( ( (so_eff / sos_eff) * mu_s_pow) + ( (ss_eff / sos_eff) * mu_o_pow) , 4.0)); const ADB mu_msg = zero_selectorSsg.select(mu_g + mu_s , mu_g * mu_s / pow( ( (sg_eff / ssg_eff) * mu_s_pow) + ( (ss_eff / ssg_eff) * mu_g_pow) , 4.0)); 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; @@ -852,9 +852,9 @@ namespace Opm { // Effective densities const ADB mu_sog_pow = mu_s_pow * ( (sgf * mu_o_pow) + (sof * mu_g_pow) ); - const ADB mu_o_eff_pow = pow(mu_o_eff,0.25); - const ADB mu_g_eff_pow = pow(mu_g_eff,0.25); - const ADB mu_s_eff_pow = pow(mu_s_eff,0.25); + const ADB mu_o_eff_pow = pow(mu_o_eff, 0.25); + const ADB mu_g_eff_pow = pow(mu_g_eff, 0.25); + const ADB mu_s_eff_pow = pow(mu_s_eff, 0.25); const ADB sfraction_oe = (mu_o_pow * (mu_o_eff_pow - mu_s_pow)) / (mu_o_eff_pow * (mu_o_pow - mu_s_pow)); const ADB sfraction_ge = (mu_g_pow * (mu_s_pow - mu_g_eff_pow)) / (mu_g_eff_pow * (mu_s_pow - mu_g_pow)); const ADB sfraction_se = (mu_sog_pow - ( mu_o_pow * mu_g_pow * mu_s_pow / mu_s_eff_pow) ) / ( mu_sog_pow - (mu_o_pow * mu_g_pow)); diff --git a/opm/autodiff/SolventPropsAdFromDeck.cpp b/opm/autodiff/SolventPropsAdFromDeck.cpp index 7bcb2f824..8f31d1524 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.cpp +++ b/opm/autodiff/SolventPropsAdFromDeck.cpp @@ -44,6 +44,7 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, // retrieve the cell specific PVT table index from the deck // and using the grid... extractPvtTableIndex(cellPvtRegionIdx_, eclState, number_of_cells, global_cell); + extractTableIndex("SATNUM", eclState, number_of_cells, global_cell, cellSatNumRegionIdx_); // surface densities if (deck->hasKeyword("SDENSITY")) { @@ -84,7 +85,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, inverseBmu[i] = 1.0 / (b[i] * visc[i]); } - b_[regionIdx] = NonuniformTableLinear(press, inverseB); viscosity_[regionIdx] = NonuniformTableLinear(press, visc); inverseBmu_[regionIdx] = NonuniformTableLinear(press, inverseBmu); @@ -99,9 +99,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, int numRegions = ssfnTables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table saturation function supported for SSFN"); - } // resize the attributes of the object krg_.resize(numRegions); krs_.resize(numRegions); @@ -123,15 +120,18 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, if (deck->hasKeyword("MISCIBLE") ) { + + + // retrieve the cell specific Misc table index from the deck + // and using the grid... + extractTableIndex("MISCNUM", eclState, number_of_cells, global_cell, cellMiscRegionIdx_); + // misicible hydrocabon relative permeability wrt water const TableContainer& sof2Tables = tables->getSof2Tables(); if (!sof2Tables.empty()) { int numRegions = sof2Tables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table saturation function supported for SOF2"); - } // resize the attributes of the object krn_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { @@ -154,9 +154,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, int numRegions = miscTables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table miscibility function supported for MISC"); - } // resize the attributes of the object misc_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { @@ -180,9 +177,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, int numRegions = msfnTables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table saturation function supported for MSFN"); - } // resize the attributes of the object mkrsg_.resize(numRegions); mkro_.resize(numRegions); @@ -206,9 +200,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, int numRegions = sorwmisTables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table miscibility function supported for SORWMIS"); - } // resize the attributes of the object sorwmis_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { @@ -227,9 +218,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, int numRegions = sgcwmisTables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table miscibility function supported for SGCWMIS"); - } // resize the attributes of the object sgcwmis_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { @@ -244,30 +232,27 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, } if (deck->hasKeyword("TLMIXPAR")) { - const auto tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(0); - std::vector mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData(); - const int numRegions = mix_params_viscosity.size(); - if (numRegions > 1) { - OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR."); - } - mix_param_viscosity_ = mix_params_viscosity[0]; + const int numRegions = deck->getKeyword("TLMIXPAR")->size(); - 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; } - - } } @@ -282,7 +267,7 @@ ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg, V dmudp(n); for (int i = 0; i < n; ++i) { const double& pg_i = pg.value()[i]; - int regionIdx = cellPvtRegionIdx_[i]; + int regionIdx = cellPvtRegionIdx_[cells[i]]; double tempInvB = b_[regionIdx](pg_i); double tempInvBmu = inverseBmu_[regionIdx](pg_i); mu[i] = tempInvB / tempInvBmu; @@ -302,35 +287,35 @@ ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg, ADB SolventPropsAdFromDeck::bSolvent(const ADB& pg, const Cells& cells) const { - return SolventPropsAdFromDeck::makeADBfromTables(pg, cells, b_); + return SolventPropsAdFromDeck::makeADBfromTables(pg, cells, cellPvtRegionIdx_, b_); } ADB SolventPropsAdFromDeck::gasRelPermMultiplier(const ADB& solventFraction, const Cells& cells) const { - return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krg_); + return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellSatNumRegionIdx_, krg_); } ADB SolventPropsAdFromDeck::solventRelPermMultiplier(const ADB& solventFraction, const Cells& cells) const { - return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krs_); + return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellSatNumRegionIdx_, krs_); } ADB SolventPropsAdFromDeck::misicibleHydrocarbonWaterRelPerm(const ADB& Sn, const Cells& cells) const { - return SolventPropsAdFromDeck::makeADBfromTables(Sn, cells, krn_); + return SolventPropsAdFromDeck::makeADBfromTables(Sn, cells, cellSatNumRegionIdx_, krn_); } ADB SolventPropsAdFromDeck::miscibleSolventGasRelPermMultiplier(const ADB& Ssg, const Cells& cells) const { if (mkrsg_.size() > 0) { - return SolventPropsAdFromDeck::makeADBfromTables(Ssg, cells, mkrsg_); + return SolventPropsAdFromDeck::makeADBfromTables(Ssg, cells, cellSatNumRegionIdx_, mkrsg_); } // trivial function if not specified return Ssg; @@ -340,7 +325,7 @@ ADB SolventPropsAdFromDeck::miscibleOilRelPermMultiplier(const ADB& So, const Cells& cells) const { if (mkro_.size() > 0) { - return SolventPropsAdFromDeck::makeADBfromTables(So, cells, mkro_); + return SolventPropsAdFromDeck::makeADBfromTables(So, cells, cellSatNumRegionIdx_, mkro_); } // trivial function if not specified return So; @@ -350,14 +335,14 @@ ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction, const Cells& cells) const { - return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, misc_); + return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellMiscRegionIdx_, misc_); } ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw, const Cells& cells) const { if (sgcwmis_.size()>0) { - return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sgcwmis_); + return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, cellMiscRegionIdx_, sgcwmis_); } // return zeros if not specified return ADB::constant(V::Zero(Sw.size())); @@ -367,7 +352,7 @@ ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw, const Cells& cells) const { if (sorwmis_.size()>0) { - return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sorwmis_); + return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, cellMiscRegionIdx_, sorwmis_); } // return zeros if not specified return ADB::constant(V::Zero(Sw.size())); @@ -375,6 +360,7 @@ ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD, const Cells& cells, + const std::vector& regionIdx, const std::vector>& tables) const { const int n = cells.size(); assert(X_AD.value().size() == n); @@ -382,9 +368,8 @@ ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD, V dx(n); for (int i = 0; i < n; ++i) { const double& X_i = X_AD.value()[i]; - int regionIdx = 0; // TODO add mapping from cells to sat function table - x[i] = tables[regionIdx](X_i); - dx[i] = tables[regionIdx].derivative(X_i); + x[i] = tables[regionIdx[cells[i]]](X_i); + dx[i] = tables[regionIdx[cells[i]]].derivative(X_i); } ADB::M dx_diag(dx.matrix().asDiagonal()); @@ -402,18 +387,56 @@ V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const { const int n = cells.size(); V density(n); for (int i = 0; i < n; ++i) { - int regionIdx = cellPvtRegionIdx_[i]; + int regionIdx = cellPvtRegionIdx_[cells[i]]; density[i] = solvent_surface_densities_[regionIdx]; } 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, + Opm::EclipseStateConstPtr eclState, + size_t numCompressed, + const int* compressedToCartesianCellIdx, + std::vector& tableIdx) const { + //Get the Region data + const std::vector& regionData = eclState->getIntGridProperty(keyword)->getData(); + // Convert this into an array of compressed cells + // Eclipse uses Fortran-style indices which start at 1 + // instead of 0, we subtract 1. + tableIdx.resize(numCompressed); + for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) { + size_t cartesianCellIdx = compressedToCartesianCellIdx ? compressedToCartesianCellIdx[cellIdx]:cellIdx; + assert(cartesianCellIdx < regionData.size()); + tableIdx[cellIdx] = regionData[cartesianCellIdx] - 1; + } } diff --git a/opm/autodiff/SolventPropsAdFromDeck.hpp b/opm/autodiff/SolventPropsAdFromDeck.hpp index 5cc4d0dab..4aa69b4c4 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: @@ -139,14 +143,30 @@ private: /// Makes ADB from table values /// \param[in] X Array of n table lookup values. /// \param[in] cells Array of n cell indices to be associated with the fraction values. - /// \param[in] tables Vector of tables, one for each PVT region. + /// \param[in] tables Vector of tables, one for each PVT region. /// \return Array of n solvent density values. ADB makeADBfromTables(const ADB& X, const Cells& cells, + const std::vector& regionIdx, const std::vector>& tables) const; + /// Helper function to create an array containing the + /// table index of for each compressed cell from an Eclipse deck. + /// \param[in] keyword eclKeyword specifying region (SATNUM etc. ) + /// \param[in] eclState eclState from opm-parser + /// \param[in] numCompressed number of compressed cells + /// \param[in] compressedToCartesianCellIdx cartesianCellIdx for each cell in the grid + /// \param[out] tableIdx table index for each compressed cell + void extractTableIndex(const std::string& keyword, + Opm::EclipseStateConstPtr eclState, + size_t numCompressed, + const int* compressedToCartesianCellIdx, + std::vector& tableIdx) const; + // The PVT region which is to be used for each cell std::vector cellPvtRegionIdx_; + std::vector cellMiscRegionIdx_; + std::vector cellSatNumRegionIdx_; std::vector > b_; std::vector > viscosity_; std::vector > inverseBmu_; @@ -159,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