diff --git a/opm/polymer/fullyimplicit/.AutoDiffBlock.hpp.swp b/opm/polymer/fullyimplicit/.AutoDiffBlock.hpp.swp new file mode 100644 index 000000000..30f7c63bd Binary files /dev/null and b/opm/polymer/fullyimplicit/.AutoDiffBlock.hpp.swp differ diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp index 5e16ce992..a101c5512 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -246,7 +246,6 @@ typedef Eigen::Array jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dads_diag * c.derivative()[block]; + } + + return ADB::function(ads, jacs); + } + + + V + PolymerPropsAd::effectiveRelPerm(const V& c, + const V& cmax_cells, + const V& krw) const + { + const int nc = c.size(); + + V one = V::Ones(nc); + V ads = adsorption(c); + double max_ads = polymer_props_.cMaxAds(); + double res_factor = polymer_props_.resFactor(); + double factor = (res_factor -1.) / max_ads; + V rk = one + factor * ads; + V krw_eff = krw / rk; + + return eff_relperm; + } + + + ADB + PolymerPropsAd::effectiveRelPerm(const ADB& c, + const ADB& cmax_cells, + const ADB& krw, + const ADB& sw) const + { + const int nc = c.value().size(); + + V one = V::Ones(nc); + + ADB ads = adsorption(c); + V krw_eff = effectiveRelPerm(c.value(), cmax_cells.value(), krw.value()); + + double max_ads = polymer_props_.cMaxAds(); + double res_factor = polymer_props_.resFactor(); + double factor = (res_factor - 1.) / max_ads; + ADB rk = one + ads * factor; + ADB::M dkrw_ds = krw.derivative() / rk.derivative(); + ADB::M dkrw_dc = -krw.value() / (rk.value * rk.value()) * ads.derivative() * factor; + + const int num_blocks = c.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jac[block] = dkrw_ds * sw.derivative()[block] + dkrw_dc * c.derivative()[block]; + } + + return ADB::function(krw_eff, jacs); + } + }// namespace Opm diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp index 69936843d..b964a241d 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp @@ -81,6 +81,17 @@ namespace Opm { ADB polymerWaterVelocityRatio(const ADB& c) const; + V + adsorption(const V& c) const; + + ADB + adsorption(const ADB& c) const; + + V + effectiveRelPerm(const V& c, const V& cmax_cells, const V& relperm) const; + + ADB + effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw, const ADB& sw) const; private: const PolymerProperties& polymer_props_; };