From a0794117c9c5e8b8a5818159cb8507c9548b1aac Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Thu, 19 Apr 2012 17:21:08 +0200 Subject: [PATCH] Cleaned code for computation of residual, many changes, not tested! --- Makefile.am | 3 +- examples/polymer_reorder.cpp | 15 +- opm/polymer/PolymerProperties.cpp | 308 ++++++++++++ opm/polymer/PolymerProperties.hpp | 199 +++----- .../SinglePointUpwindTwoPhasePolymer.hpp | 30 +- opm/polymer/TransportModelPolymer.cpp | 444 ++++++++---------- opm/polymer/TransportModelPolymer.hpp | 11 +- opm/polymer/polymerUtilities.cpp | 4 +- 8 files changed, 605 insertions(+), 409 deletions(-) create mode 100644 opm/polymer/PolymerProperties.cpp diff --git a/Makefile.am b/Makefile.am index 241ae66b6..8b168d890 100644 --- a/Makefile.am +++ b/Makefile.am @@ -8,7 +8,8 @@ lib_LTLIBRARIES = libopmpolymer.la libopmpolymer_la_SOURCES = \ opm/polymer/TransportModelPolymer.cpp \ -opm/polymer/polymerUtilities.cpp +opm/polymer/PolymerProperties.cpp \ +opm/polymer/polymerUtilities.cpp nobase_include_HEADERS = \ opm/polymer/GravityColumnSolverPolymer.hpp \ diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index b00c3daa7..9556f2f93 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -169,7 +169,7 @@ public: void adsorption(const PolyC& c, const PolyC& cmax, CAds& cads, DCAdsDc& dcadsdc) { - cads = polyprops_.adsorptionWithDer(c, cmax, &dcadsdc); + polyprops_.adsorptionWithDer(c, cmax, cads, dcadsdc); } const double* porosity() const @@ -187,14 +187,14 @@ public: class Mob, class DMobDs, class DMobWatDc> - void mobility(int cell, const Sat& s, const PolyC& c, + void mobility(int cell, const Sat& s, const PolyC& c, const PolyC& cmax, Mob& mob, DMobDs& dmobds, DMobWatDc& dmobwatdc) const { const double* visc = props_.viscosity(); double relperm[2]; double drelpermds[4]; props_.relperm(1, &s[0], &cell, relperm, drelpermds); - polyprops_.effectiveMobilities(c, visc, relperm, drelpermds, mob, dmobds, dmobwatdc); + polyprops_.effectiveMobilitiesWithDer(c, cmax, visc, relperm, drelpermds, mob, dmobds, dmobwatdc); } template - void mc(const PolyC& c, Mc& mcval, - DMcDc& dmcdcval) const + void computeMc(const PolyC& c, Mc& mc, + DMcDc& dmcdc) const { - polyprops_.computeMc(c, mcval, dmcdcval); + polyprops_.computeMcWithDer(c, mc, dmcdc); } private: @@ -512,7 +512,7 @@ main(int argc, char** argv) // reorder_model.initGravity(grav); } // Non-reordering solver. - NewtonPolymerTransportModel model (fluid, *grid->c_grid(), porevol, grav, guess_old_solution); + NewtonPolymerTransportModel model(fluid, *grid->c_grid(), porevol, grav, guess_old_solution); if (use_gravity) { model.initGravityTrans(*grid->c_grid(), psolver.getHalfTrans()); } @@ -526,7 +526,6 @@ main(int argc, char** argv) Opm::GravityColumnSolverPolymer colsolver(model, *grid->c_grid(), nl_tolerance, nl_maxiter); - // // // Not implemented for polymer. // // Control init. // Opm::ImplicitTransportDetails::NRReport rpt; diff --git a/opm/polymer/PolymerProperties.cpp b/opm/polymer/PolymerProperties.cpp new file mode 100644 index 000000000..8cd839287 --- /dev/null +++ b/opm/polymer/PolymerProperties.cpp @@ -0,0 +1,308 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include +#include +#include + +namespace Opm +{ + double PolymerProperties::cMax() const + { + return c_max_; + } + + double PolymerProperties::mixParam() const + { + return mix_param_; + } + + double PolymerProperties::rockDensity() const + { + return rock_density_; + } + + double PolymerProperties::deadPoreVol() const + { + return dead_pore_vol_; + } + + double PolymerProperties::resFactor() const + { + return res_factor_; + } + + double PolymerProperties::cMaxAds() const + { + return c_max_ads_; + } + + int PolymerProperties::adsIndex() const + { + return ads_index_; + } + + double PolymerProperties::viscMult(double c) const + { + return Opm::linearInterpolation(c_vals_visc_, visc_mult_vals_, c); + } + + double PolymerProperties::viscMultWithDer(double c, double* der) const + { + *der = Opm::linearInterpolationDerivative(c_vals_visc_, visc_mult_vals_, c); + return Opm::linearInterpolation(c_vals_visc_, visc_mult_vals_, c); + } + + void PolymerProperties::simpleAdsorption(double c, double& c_ads) const + { + double dummy; + simpleAdsorptionBoth(c, c_ads, dummy, false); + } + + void PolymerProperties::simpleAdsorptionWithDer(double c, double& c_ads, + double& dc_ads_dc) const + { + simpleAdsorptionBoth(c, c_ads, dc_ads_dc, true); + } + + void PolymerProperties::simpleAdsorptionBoth(double c, double& c_ads, + double& dc_ads_dc, bool if_with_der) const + { + c_ads = Opm::linearInterpolation(c_vals_ads_, ads_vals_, c);; + if (if_with_der) { + dc_ads_dc = Opm::linearInterpolationDerivative(c_vals_ads_, ads_vals_, c); + } else { + dc_ads_dc = 0.; + } + } + + void PolymerProperties::adsorption(double c, double cmax, double& c_ads) const + { + double dummy; + adsorptionBoth(c, cmax, c_ads, dummy, false); + } + + void PolymerProperties::adsorptionWithDer(double c, double cmax, + double& c_ads, double& dc_ads_dc) const + { + adsorptionBoth(c, cmax, c_ads, dc_ads_dc, true); + } + + void PolymerProperties::adsorptionBoth(double c, double cmax, + double& c_ads, double& dc_ads_dc, + bool if_with_der) const + { + if (ads_index_ == Desorption) { + simpleAdsorptionBoth(c, c_ads, dc_ads_dc, if_with_der); + } else if (ads_index_ == NoDesorption) { + if (c < cmax) { + simpleAdsorption(cmax, c_ads); + dc_ads_dc = 0.; + } else { + simpleAdsorptionBoth(c, c_ads, dc_ads_dc, if_with_der); + } + } else { + THROW("Invalid Adsoption index"); + } + } + + + void PolymerProperties::effectiveVisc(const double c, const double* visc, double* visc_eff) const { + effectiveViscBoth(c, visc, visc_eff, 0, false); + } + + void PolymerProperties::effectiveInvVisc(const double c, const double* visc, double* inv_visc_eff) const + { + effectiveViscBoth(c, visc, inv_visc_eff, 0, false); + inv_visc_eff[0] = 1./inv_visc_eff[0]; + inv_visc_eff[1] = 1./inv_visc_eff[1]; + } + + void PolymerProperties::effectiveViscWithDer(const double c, const double* visc, double* visc_eff, + double* dvisc_eff_dc) const { + effectiveViscBoth(c, visc, visc_eff, dvisc_eff_dc, true); + } + + void PolymerProperties::effectiveViscBoth(const double c, const double* visc, double* visc_eff, + double* dvisc_eff_dc, bool if_with_der) const + { + double cbar = c/c_max_; + double mu_w = visc[0]; + double mu_m; + double omega = mix_param_; + double mu_m_dc; + if (if_with_der) { + mu_m = viscMultWithDer(c, &mu_m_dc)*mu_w; + mu_m_dc *= mu_w; + } else { + mu_m = viscMult(c)*mu_w; + } + double mu_p = viscMult(c_max_)*mu_w; + double mu_m_omega = std::pow(mu_m, mix_param_); + double mu_w_e = mu_m_omega*std::pow(mu_w, 1.0 - mix_param_); + double mu_p_eff = mu_m_omega*std::pow(mu_p, 1.0 - mix_param_); + double mu_w_eff = 1./((1.0 - cbar)/mu_w_e + cbar/mu_p_eff); + visc_eff[0] = mu_w_eff; + visc_eff[1] = visc[1]; + if (if_with_der) { + double mu_w_e_dc = omega*mu_m_dc*std::pow(mu_m, omega - 1)*std::pow(mu_w, 1 - omega); + double mu_p_eff_dc = omega*mu_m_dc*std::pow(mu_m, omega - 1)*std::pow(mu_p, 1 - omega); + dvisc_eff_dc[0] = -1./c_max_*mu_w_eff*mu_w_eff*(1./mu_p_eff - 1./mu_w_e) + (1-cbar)*(mu_w_eff*mu_w_eff/(mu_w_e*mu_w_e))*mu_w_e_dc + cbar*(mu_w_eff*mu_w_eff/(mu_p_eff*mu_p_eff))*mu_p_eff_dc; + dvisc_eff_dc[1] = 0.; + } else { + dvisc_eff_dc = 0; + } + } + + void PolymerProperties::effectiveRelperm(const double c, + const double cmax, + const double* relperm, + double& eff_relperm_wat) const { + double dummy; + effectiveRelpermBoth(c, cmax, relperm, 0, eff_relperm_wat, + dummy, dummy, false); + } + + void PolymerProperties::effectiveRelpermWithDer (const double c, + const double cmax, + const double* relperm, + const double* drelperm_ds, + double& eff_relperm_wat, + double& deff_relperm_wat_ds, + double& deff_relperm_wat_dc) const { + effectiveRelpermBoth(c, cmax, relperm, + drelperm_ds, eff_relperm_wat, + deff_relperm_wat_ds, deff_relperm_wat_dc, + true); + } + + void PolymerProperties::effectiveRelpermBoth(const double c, + const double cmax, + const double* relperm, + const double* drelperm_ds, + double& eff_relperm_wat, + double& deff_relperm_wat_ds, + double& deff_relperm_wat_dc, + bool if_with_der) const { + double c_ads; + double dc_ads_dc; + adsorptionBoth(c, cmax, c_ads, dc_ads_dc, if_with_der); + double rk = 1 + (res_factor_ - 1)*c_ads/c_max_ads_; + eff_relperm_wat = relperm[0]/rk; + if (if_with_der) { + deff_relperm_wat_ds = drelperm_ds[0]/rk; + deff_relperm_wat_dc = dc_ads_dc*relperm[0]/(rk*rk*c_max_ads_); + } else { + deff_relperm_wat_ds = -1.0; + deff_relperm_wat_dc = -1.0; + } + } + + void PolymerProperties::effectiveMobilities(const double c, + const double cmax, + const double* visc, + const double* relperm, + std::vector& mob) const + { + double dummy1; + double dummy2[2]; + std::vector dummy3; + effectiveMobilitiesBoth(c, cmax, visc, relperm, + dummy2, mob, dummy3, dummy1, false); + } + + void PolymerProperties::effectiveMobilitiesWithDer(const double c, + const double cmax, + const double* visc, + const double* relperm, + const double* drelpermds, + std::vector& mob, + std::vector& dmobds, + double& dmobwatdc) const + { + effectiveMobilitiesBoth(c, cmax, visc, + relperm, drelpermds, mob,dmobds, + dmobwatdc, true); + } + + void PolymerProperties::effectiveMobilitiesBoth(const double c, + const double cmax, + const double* visc, + const double* relperm, + const double* drelperm_ds, + std::vector& mob, + std::vector& dmob_ds, + double& dmobwat_dc, + bool if_with_der) const + { + double visc_eff[2]; + double dvisc_eff_dc[2]; + effectiveViscBoth(c, visc, visc_eff, dvisc_eff_dc, if_with_der); + double mu_w_eff = visc_eff[0]; + double mu_w_eff_dc = dvisc_eff_dc[0]; + double eff_relperm_wat; + double deff_relperm_wat_ds; + double deff_relperm_wat_dc; + effectiveRelpermBoth(c, cmax, relperm, + drelperm_ds, eff_relperm_wat, + deff_relperm_wat_ds, deff_relperm_wat_dc, + if_with_der); + mob[0] = eff_relperm_wat/visc_eff[0]; + mob[1] = relperm[1]/visc_eff[1]; + + if (if_with_der) { + dmobwat_dc = - mob[0]*mu_w_eff_dc/(mu_w_eff*mu_w_eff) + + deff_relperm_wat_dc/mu_w_eff; + dmob_ds[0*2 + 0] = deff_relperm_wat_ds/visc_eff[0]; + dmob_ds[0*2 + 1] = -dmob_ds[0*2 + 0]; + dmob_ds[1*2 + 0] = drelperm_ds[1*2 + 0]/visc_eff[1]; + dmob_ds[1*2 + 1] = drelperm_ds[1*2 + 1]/visc_eff[1]; + } else { + mob.clear(); + dmob_ds.clear(); + } + } + + void PolymerProperties::computeMcWithDer(const double& c, double& mc) const + { + double dummy; + computeMcBoth(c, mc, dummy, false); + } + + void PolymerProperties::computeMcWithDer(const double& c, double& mc, + double& dmc_dc) const + { + computeMcBoth(c, mc, dmc_dc, true); + } + + void PolymerProperties::computeMcBoth(const double& c, double& mc, + double& dmc_dc, bool if_with_der) const + { + double cbar = c/c_max_; + double omega = mix_param_; + double r = std::pow(viscMult(c_max_), 1 - omega); // viscMult(c_max_)=mu_p/mu_w + mc = c/(cbar + (1 - cbar)*r); + if (if_with_der) { + dmc_dc = r/std::pow(cbar + (1 - cbar)*r, 2); + } else { + dmc_dc = 0.; + } + } +} diff --git a/opm/polymer/PolymerProperties.hpp b/opm/polymer/PolymerProperties.hpp index 91c185e98..27ee67ca6 100644 --- a/opm/polymer/PolymerProperties.hpp +++ b/opm/polymer/PolymerProperties.hpp @@ -23,7 +23,6 @@ #include #include -#include #include @@ -126,152 +125,87 @@ namespace Opm } - double cMax() const - { - return c_max_; - } + double cMax() const; - double mixParam() const - { - return mix_param_; - } + double mixParam() const; - double rockDensity() const - { - return rock_density_; - }; + double rockDensity() const; - double deadPoreVol() const - { - return dead_pore_vol_; - } + double deadPoreVol() const; - double resFactor() const - { - return res_factor_; - } + double resFactor() const; - double cMaxAds() const - { - return c_max_ads_; - } + double cMaxAds() const; - int adsIndex() const - { - return ads_index_; - } + int adsIndex() const; - double viscMult(double c) const - { - return Opm::linearInterpolation(c_vals_visc_, visc_mult_vals_, c); - } + double viscMult(double c) const; - double viscMultWithDer(double c, double* der) const - { - *der = Opm::linearInterpolationDerivative(c_vals_visc_, visc_mult_vals_, c); - return Opm::linearInterpolation(c_vals_visc_, visc_mult_vals_, c); - } + double viscMultWithDer(double c, double* der) const; - double simpleAdsorption(double c) const - { - return Opm::linearInterpolation(c_vals_ads_, ads_vals_, c); - } + void simpleAdsorption(double c, double& c_ads) const; - double simpleAdsorptionWithDer(double c, double* der) const - { - *der = Opm::linearInterpolationDerivative(c_vals_ads_, ads_vals_, c); - return Opm::linearInterpolation(c_vals_ads_, ads_vals_, c); - } + void simpleAdsorptionWithDer(double c, double& c_ads, + double& dc_ads_dc) const; - double adsorption(double c, double cmax) const - { - if (ads_index_ == Desorption) { - return simpleAdsorption(c); - } else if (ads_index_ == NoDesorption) { - return simpleAdsorption(std::max(c, cmax)); - } else { - THROW("Invalid Adsoption index"); - } - } + void adsorption(double c, double cmax, double& c_ads) const; - double adsorptionWithDer(double c, double cmax, double* der) const - { - if (ads_index_ == Desorption) { - return simpleAdsorptionWithDer(c, der); - } else if (ads_index_ == NoDesorption) { - if (c < cmax) { - *der = 0; - return simpleAdsorption(cmax); - } else { - return simpleAdsorptionWithDer(c, der); - } - } else { - THROW("Invalid Adsoption index"); - } - } + void adsorptionWithDer(double c, double cmax, + double& c_ads, double& dc_ads_dc) const; + void effectiveVisc(const double c, const double* visc, double* visc_eff) const; + + void effectiveInvVisc(const double c, const double* visc, + double* inv_visc_eff) const; + + void effectiveViscWithDer(const double c, const double* visc, + double* visc_eff, double* dvisc_eff_dc) const; + + void effectiveRelperm(const double c, + const double cmax, + const double* relperm, + double& eff_relperm_wat) const; + + void effectiveRelpermWithDer (const double c, + const double cmax, + const double* relperm, + const double* drelperm_ds, + double& eff_relperm_wat, + double& deff_relperm_wat_ds, + double& deff_relperm_wat_dc) const; - void effectiveInvVisc(const double c, const double* visc, double* inv_visc_eff) const - { - double cbar = c/c_max_; - double mu_w = visc[0]; - double mu_m = viscMult(c)*mu_w; - double mu_p = viscMult(c_max_)*mu_w; - double mu_m_omega = std::pow(mu_m, mix_param_); - double mu_w_e = mu_m_omega*std::pow(mu_w, 1.0 - mix_param_); - double mu_p_eff = mu_m_omega*std::pow(mu_p, 1.0 - mix_param_); - double inv_mu_w_eff = (1.0 - cbar)/mu_w_e + cbar/mu_p_eff; - inv_visc_eff[0] = inv_mu_w_eff; - inv_visc_eff[1] = 1.0/visc[1]; - } - void effectiveMobilities(const double c, + const double cmax, const double* visc, const double* relperm, - const double* drelpermds, - std::vector& mob, - std::vector& dmobds, - double& dmobwatdc) const - { - double cbar = c/c_max_; - double mu_w = visc[0]; - double mu_m_dc; // derivative of mu_m with respect to c - double mu_m = viscMultWithDer(c, &mu_m_dc)*mu_w; - mu_m_dc *= mu_w; - double mu_p = viscMult(c_max_)*mu_w; - double omega = mix_param_; - double mu_w_e = std::pow(mu_m, omega)*std::pow(mu_w, 1 - omega); - double mu_w_e_dc = omega*mu_m_dc*std::pow(mu_m, omega - 1)*std::pow(mu_w, 1 - omega); - double mu_p_eff = std::pow(mu_m, omega)*std::pow(mu_p, 1 - omega); - double mu_p_eff_dc = omega*mu_m_dc*std::pow(mu_m, omega - 1)*std::pow(mu_p, 1 - omega); - double mu_w_eff = 1./((1 - cbar)/mu_w_e + cbar/mu_p_eff); - double mu_w_eff_dc = -1./c_max_*mu_w_eff*mu_w_eff*(1./mu_p_eff - 1./mu_w_e) - + (1-cbar)*(mu_w_eff*mu_w_eff/(mu_w_e*mu_w_e))*mu_w_e_dc - + cbar*(mu_w_eff*mu_w_eff/(mu_p_eff*mu_p_eff))*mu_p_eff_dc; - double visc_eff[2] = { mu_w_eff, visc[1] }; + std::vector& mob) const; - dmobwatdc = - mob[0]*mu_w_eff_dc/(mu_w_eff*mu_w_eff); + void effectiveMobilitiesWithDer(const double c, + const double cmax, + const double* visc, + const double* relperm, + const double* drelpermds, + std::vector& mob, + std::vector& dmobds, + double& dmobwatdc) const; - mob[0] = relperm[0]/visc_eff[0]; - mob[1] = relperm[1]/visc_eff[1]; + void effectiveMobilitiesBoth(const double c, + const double cmax, + const double* visc, + const double* relperm, + const double* drelperm_ds, + std::vector& mob, + std::vector& dmob_ds, + double& dmobwat_dc, + bool if_with_der) const; - dmobds[0*2 + 0] = drelpermds[0*2 + 0]/visc_eff[0]; - dmobds[0*2 + 1] = drelpermds[0*2 + 1]/visc_eff[1]; - dmobds[1*2 + 0] = drelpermds[1*2 + 0]/visc_eff[0]; - dmobds[1*2 + 1] = drelpermds[1*2 + 1]/visc_eff[1]; - } + void computeMcWithDer(const double& c, double& mc) const; - void computeMc(const double& c, - double& mc, - double& dmcdc) const - { - double cbar = c/c_max_; - double omega = mix_param_; - double r = std::pow(viscMult(c_max_), 1 - omega); // viscMult(c_max_)=mu_p/mu_w - mc = c/(cbar + (1 - cbar)*r); - dmcdc = r/std::pow(cbar + (1 - cbar)*r, 2); - } + void computeMcWithDer(const double& c, double& mc, + double& dmc_dc) const; + void computeMcBoth(const double& c, double& mc, + double& dmc_dc, bool if_with_der) const; private: double c_max_; @@ -285,6 +219,21 @@ namespace Opm std::vector visc_mult_vals_; std::vector c_vals_ads_; std::vector ads_vals_; + void simpleAdsorptionBoth(double c, double& c_ads, + double& dc_ads_dc, bool if_with_der) const; + void adsorptionBoth(double c, double cmax, + double& c_ads, double& dc_ads_dc, + bool if_with_der) const; + void effectiveViscBoth(const double c, const double* visc, double* visc_eff, + double* dvisc_eff_dc, bool if_with_der) const; + void effectiveRelpermBoth(const double c, + const double cmax, + const double* relperm, + const double* drelperm_ds, + double& eff_relperm_wat, + double& deff_relperm_wat_ds, + double& deff_relperm_wat_dc, + bool if_with_der) const; }; } // namespace Opm diff --git a/opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp b/opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp index e006f82b3..00c3853c2 100644 --- a/opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp +++ b/opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp @@ -48,9 +48,11 @@ namespace Opm { class ModelParameterStorage { public: ModelParameterStorage(int nc, int totconn) - : drho_(0.0), deadporespace_(0.0), rockdensity_(0.0), mob_(0), dmobds_(0), dmobwatdc_(0), mc_(0), - dmcdc_(0), porevol_(0), porosity_(0), dg_(0), sw_(0), c_(0), - ds_(0), dsc_(0), dcads_(0), dcadsdc_(0), pc_(0), dpc_(0), trans_(0), data_() + : drho_(0.0), deadporespace_(0.0), rockdensity_(0.0), mob_(0), + dmobds_(0), dmobwatdc_(0), mc_(0), + dmcdc_(0), porevol_(0), porosity_(0), dg_(0), sw_(0), c_(0), cmax_(0), + ds_(0), dsc_(0), dcads_(0), dcadsdc_(0), pc_(0), dpc_(0), + trans_(0), data_() { size_t alloc_sz; @@ -65,6 +67,7 @@ namespace Opm { alloc_sz += 1 * totconn; // dg_ alloc_sz += 1 * nc; // sw_ alloc_sz += 1 * nc; // c_ + alloc_sz += 1 * nc; // cmax_ alloc_sz += 1 * nc; // dc_ alloc_sz += 1 * nc; // dsc_ alloc_sz += 1 * nc; // dcads_ @@ -84,7 +87,8 @@ namespace Opm { dg_ = porosity_ + (1 * nc ); sw_ = dg_ + (1 * totconn); c_ = sw_ + (1 * nc ); - ds_ = c_ + (1 * nc ); + cmax_ = c_ + (1 * nc ); + ds_ = cmax_ + (1 * nc ); dsc_ = ds_ + (1 * nc ); dcads_ = dsc_ + (1 * nc ); dcadsdc_ = dcads_ + (1 * nc ); @@ -133,6 +137,9 @@ namespace Opm { double& c(int cell) { return c_[cell] ; } double c(int cell) const { return c_[cell] ; } + double& cmax(int cell) { return cmax_[cell] ; } + double cmax(int cell) const { return cmax_[cell] ; } + double& ds(int cell) { return ds_[cell] ; } double ds(int cell) const { return ds_[cell] ; } @@ -168,6 +175,7 @@ namespace Opm { double *dg_ ; double *sw_ ; double *c_ ; + double *cmax_ ; double *ds_ ; double *dsc_ ; double *dcads_ ; // difference of cads to compute residual @@ -495,7 +503,7 @@ namespace Opm { std::vector mob(2, 0.); std::vector dmobds(4, 0.); double dmobwatdc; - double c; + double c, cmax; double mc, dmcdc; double pc, dpc; @@ -503,7 +511,7 @@ namespace Opm { sys.vector().solution(); const ::std::vector& sat = state.saturation(); const ::std::vector& cpoly = state.concentration(); - const ::std::vector& cmax = state.maxconcentration(); + const ::std::vector& cmaxpoly = state.maxconcentration(); bool in_range = true; for (int cell = 0; cell < g.number_of_cells; ++cell) { @@ -511,14 +519,16 @@ namespace Opm { store_.ds(cell) = x[2*cell + 0]; s[0] = sat[cell*2 + 0] + x[2*cell + 0]; c = cpoly[cell] + x[2*cell + 1]; + cmax = std::max(cpoly[cell] + x[2*cell + 1], cmaxpoly[cell]); store_.sw(cell) = s[0]; store_.c(cell) = c; + store_.cmax(cell) = cmax; store_.dsc(cell) = s[0]*c - sat[cell*2 + 0]*cpoly[cell]; double dcadsdc; double cads; - fluid_.adsorption(cpoly[cell], cmax[cell], cads, dcadsdc); + fluid_.adsorption(cpoly[cell], cmax, cads, dcadsdc); store_.dcads(cell) = -cads; - fluid_.adsorption(c, cmax[cell], cads, dcadsdc); + fluid_.adsorption(c, cmax, cads, dcadsdc); store_.dcads(cell) += cads; store_.dcadsdc(cell) = dcadsdc; double s_min = fluid_.s_min(cell); @@ -537,8 +547,8 @@ namespace Opm { s[0] = std::min(s_max, s[0]); s[1] = 1 - s[0]; - fluid_.mobility(cell, s, c, mob, dmobds, dmobwatdc); - fluid_.mc(c, mc, dmcdc); + fluid_.mobility(cell, s, c, cmax, mob, dmobds, dmobwatdc); + fluid_.computeMc(c, mc, dmcdc); fluid_.pc(cell, s, pc, dpc); store_.mob (cell)[0] = mob [0]; diff --git a/opm/polymer/TransportModelPolymer.cpp b/opm/polymer/TransportModelPolymer.cpp index 0292b665e..85438c4be 100644 --- a/opm/polymer/TransportModelPolymer.cpp +++ b/opm/polymer/TransportModelPolymer.cpp @@ -52,8 +52,13 @@ public: double computeResidualC(const double* x) const; void computeGradientResS(const double* x, double* res, double* gradient) const; void computeGradientResC(const double* x, double* res, double* gradient) const; - void computeJacobiRes(const double* x, double* res_s_ds_dc, double* res_c_ds_dc) const; - + void computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const; + +private: + void computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c, + const bool if_dres_s_dsdc, const bool if_dres_c_dsdc, + double* res, double* dres_s_dsdc, + double* dres_c_dsdc, double& mc, double& ff) const; }; @@ -229,7 +234,9 @@ namespace Opm } double operator()(double s) const { - return s - s0_ + dtpv_*(outflux_*tm_.fracFlow(s, c_, cell_) + influx_ + s*comp_term_); + double ff; + tm_.fracFlow(s, c_, cmax0_, cell_, ff); + return s - s0_ + dtpv_*(outflux_*ff + influx_ + s*comp_term_); } }; @@ -264,7 +271,9 @@ namespace Opm double dflux = -tm.source_[cell]; bool src_is_inflow = dflux < 0.0; influx = src_is_inflow ? dflux : 0.0; - influx_polymer = src_is_inflow ? dflux*tm.computeMc(tm.inflow_c_) : 0.0; + double mc; + tm.computeMc(tm.inflow_c_, mc); + influx_polymer = src_is_inflow ? dflux*mc : 0.0; outflux = !src_is_inflow ? dflux : 0.0; comp_term = tm.source_[cell]; // Note: this assumes that all source flux is water. dtpv = tm.dt_/tm.porevolume_[cell]; @@ -299,16 +308,18 @@ namespace Opm void computeBothResiduals(const double s_arg, const double c_arg, double& res_s, double& res_c, double& mc, double& ff) const { double dps = tm.polyprops_.deadPoreVol(); - ff = tm.fracFlow(s_arg, c_arg, cell); - mc = tm.computeMc(c_arg); + tm.fracFlow(s_arg, c_arg, cmax0, cell, ff); + tm.computeMc(c_arg, mc); double rhor = tm.polyprops_.rockDensity(); - double ads0 = tm.polyprops_.adsorption(c0, cmax0); - double ads = tm.polyprops_.adsorption(c_arg, cmax0); + double c_ads0; + tm.polyprops_.adsorption(c0, cmax0, c_ads0); + double c_ads; + tm.polyprops_.adsorption(c_arg, cmax0, c_ads); res_s = s_arg - s0 + dtpv*(outflux*ff + influx + s*comp_term); res_c = s_arg*(1 - dps)*c_arg - (s0 - dps)*c0 - + rhor*((1.0 - porosity)/porosity)*(ads - ads0) + + rhor*((1.0 - porosity)/porosity)*(c_ads - c_ads0) + dtpv*(outflux*ff*mc + influx_polymer) - + dtpv*(s_arg*c_arg*(1.0 - dps) - rhor*ads)*comp_term; + + dtpv*(s_arg*c_arg*(1.0 - dps) - rhor*c_ads)*comp_term; } @@ -322,15 +333,19 @@ namespace Opm // tm.maxit_, tm.tol_, iters_used); s = modifiedRegulaFalsi(res_s, s0, 0.0, 1.0, tm.maxit_, tm.tol_, iters_used); - double ff = tm.fracFlow(s, c, cell); - double mc = tm.computeMc(c); + double ff; + tm.fracFlow(s, c, cmax0, cell, ff); + double mc; + tm.computeMc(c, mc); double rhor = tm.polyprops_.rockDensity(); - double ads0 = tm.polyprops_.adsorption(c0, cmax0); - double ads = tm.polyprops_.adsorption(c, cmax0); + double c_ads0; + tm.polyprops_.adsorption(c0, cmax0, c_ads0); + double c_ads; + tm.polyprops_.adsorption(c, cmax0, c_ads); double res = (1 - dps)*s*c - (1 - dps)*s0*c0 - + rhor*((1.0 - porosity)/porosity)*(ads - ads0) + + rhor*((1.0 - porosity)/porosity)*(c_ads - c_ads0) + dtpv*(outflux*ff*mc + influx_polymer) - + dtpv*(s*c*(1.0 - dps) - rhor*ads)*comp_term; + + dtpv*(s*c*(1.0 - dps) - rhor*c_ads)*comp_term; #ifdef EXTRA_DEBUG_OUTPUT std::cout << "c = " << c << " s = " << s << " c-residual = " << res << std::endl; #endif @@ -357,13 +372,15 @@ namespace Opm cmax0 = tm.cmax_[cell]; dps = tm.polyprops_.deadPoreVol(); rhor = tm.polyprops_.rockDensity(); - ads0 = tm.polyprops_.adsorption(c0, cmax0); + tm.polyprops_.adsorption(c0, cmax0, ads0); res_factor = tm.polyprops_.resFactor(); c_max_ads = tm.polyprops_.cMaxAds(); double dflux = -tm.source_[cell]; bool src_is_inflow = dflux < 0.0; influx = src_is_inflow ? dflux : 0.0; - influx_polymer = src_is_inflow ? dflux*tm.computeMc(tm.inflow_c_) : 0.0; + double mc; + tm.computeMc(tm.inflow_c_, mc); + influx_polymer = src_is_inflow ? dflux*mc : 0.0; outflux = !src_is_inflow ? dflux : 0.0; dtpv = tm.dt_/tm.porevolume_[cell]; porosity = tm.porosity_[cell]; @@ -391,176 +408,141 @@ namespace Opm } } + void TransportModelPolymer::ResidualEquation::computeResidual(const double* x, double* res) const { - double s = x[0]; - double c = x[1]; - double ff = tm.fracFlow(s, c, cell); - double mc = tm.computeMc(c); - // double dps = tm.polyprops_.deadPoreVol(); - // double rhor = tm.polyprops_.rockDensity(); - // double ads0 = tm.polyprops_.adsorption(c0, cmax0); - double ads = tm.polyprops_.adsorption(c, cmax0); - res[0] = s - s0 + dtpv*(outflux*ff + influx); - res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0 - + rhor*((1.0 - porosity)/porosity)*(ads - ads0) - + dtpv*(outflux*ff*mc + influx_polymer); + double dummy; + computeResAndJacobi(x, true, true, false, false, res, 0, 0, dummy, dummy); } void TransportModelPolymer::ResidualEquation::computeResidual(const double* x, double* res, double& mc, double& ff) const { - double s = x[0]; - double c = x[1]; - ff = tm.fracFlow(s, c, cell); - mc = tm.computeMc(c); - // double dps = tm.polyprops_.deadPoreVol(); - // double rhor = tm.polyprops_.rockDensity(); - // double ads0 = tm.polyprops_.adsorption(c0, cmax0); - double ads = tm.polyprops_.adsorption(c, cmax0); - res[0] = s - s0 + dtpv*(outflux*ff + influx); - res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0 - + rhor*((1.0 - porosity)/porosity)*(ads - ads0) - + dtpv*(outflux*ff*mc + influx_polymer); + computeResAndJacobi(x, true, true, false, false, res, 0, 0, mc, ff); } double TransportModelPolymer::ResidualEquation::computeResidualS(const double* x) const { - double s = x[0]; - double c = x[1]; - double ff = tm.fracFlow(s, c, cell); - return s - s0 + dtpv*(outflux*ff + influx); + double res[2]; + double dummy; + computeResAndJacobi(x, true, false, false, false, res, 0, 0, dummy, dummy); + return res[0]; } double TransportModelPolymer::ResidualEquation::computeResidualC(const double* x) const { - double s = x[0]; - double c = x[1]; - double ff = tm.fracFlow(s, c, cell); - double mc = tm.computeMc(c); - double ads = tm.polyprops_.adsorption(c, cmax0); - return (1 - dps)*s*c - (1 - dps)*s0*c0 - + rhor*((1.0 - porosity)/porosity)*(ads - ads0) - + dtpv*(outflux*ff*mc + influx_polymer); + double res[2]; + double dummy; + computeResAndJacobi(x, false, true, false, false, res, 0, 0, dummy, dummy); + return res[1]; } void TransportModelPolymer::ResidualEquation::computeGradientResS(const double* x, double* res, double* gradient) const // If gradient_method == FinDif, use finite difference // If gradient_method == Analytic, use analytic expresions { - if (gradient_method == FinDif) { - double epsi = 1e-8; - double res_epsi[2]; - double x_epsi[2]; - computeResidual(x, res); - x_epsi[0] = x[0] + epsi; - x_epsi[1] = x[1]; - computeResidual(x_epsi, res_epsi); - gradient[0] = (res_epsi[0] - res[0])/epsi; - x_epsi[0] = x[0]; - x_epsi[1] = x[1] + epsi; - computeResidual(x_epsi, res_epsi); - gradient[1] = (res_epsi[0] - res[0])/epsi; - } else if (gradient_method == Analytic) { - double s = x[0]; - double c = x[1]; - double ff_ds_dc[2]; - double ff = tm.fracFlowWithDer(s, c, cell, ff_ds_dc); - double mc_dc; - double mc = tm.computeMcWithDer(c, &mc_dc); - double ads_dc; - double ads = tm.polyprops_.adsorptionWithDer(c, cmax0, &ads_dc); - res[0] = s - s0 + dtpv*(outflux*ff + influx); - res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0 - + rhor*((1.0 - porosity)/porosity)*(ads - ads0) - + dtpv*(outflux*ff*mc + influx_polymer); - gradient[0] = 1 + dtpv*outflux*ff_ds_dc[0]; - gradient[1] = dtpv*outflux*ff_ds_dc[1]; - } + double dummy; + computeResAndJacobi(x, true, true, true, false, res, gradient, 0, dummy, dummy); } void TransportModelPolymer::ResidualEquation::computeGradientResC(const double* x, double* res, double* gradient) const // If gradient_method == FinDif, use finite difference // If gradient_method == Analytic, use analytic expresions { - if (gradient_method == FinDif) { - double epsi = 1e-8; - double res_epsi[2]; - double x_epsi[2]; - computeResidual(x, res); - x_epsi[0] = x[0] + epsi; - x_epsi[1] = x[1]; - computeResidual(x_epsi, res_epsi); - gradient[0] = (res_epsi[1] - res[1])/epsi; - x_epsi[0] = x[0]; - x_epsi[1] = x[1] + epsi; - computeResidual(x_epsi, res_epsi); - gradient[1] = (res_epsi[1] - res[1])/epsi; - } else if (gradient_method == Analytic) { - double s = x[0]; - double c = x[1]; - double ff_ds_dc[2]; - double ff = tm.fracFlowWithDer(s, c, cell, ff_ds_dc); - double mc_dc; - double mc = tm.computeMcWithDer(c, &mc_dc); - // double dps = tm.polyprops_.deadPoreVol(); - // double rhor = tm.polyprops_.rockDensity(); - // double ads0 = tm.polyprops_.adsorption(c0, cmax0); - double ads_dc; - double ads = tm.polyprops_.adsorptionWithDer(c, cmax0, &ads_dc); - res[0] = s - s0 + dtpv*(outflux*ff + influx); - res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0 - + rhor*((1.0 - porosity)/porosity)*(ads - ads0) - + dtpv*(outflux*ff*mc + influx_polymer); - gradient[0] = (1 - dps)*c + dtpv*outflux*(ff_ds_dc[0])*mc; - gradient[1] = (1 - dps)*s + rhor*((1.0 - porosity)/porosity)*ads_dc - + dtpv*outflux*(ff_ds_dc[1]*mc + ff*mc_dc); - } + double dummy; + computeResAndJacobi(x, true, true, false, true, res, 0, gradient, dummy, dummy); } // Compute the Jacobian of the residual equations. - void TransportModelPolymer::ResidualEquation::computeJacobiRes(const double* x, double* res_s_ds_dc, double* res_c_ds_dc) const + void TransportModelPolymer::ResidualEquation::computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const { - if (gradient_method == FinDif) { + double dummy; + computeResAndJacobi(x, false, false, true, true, 0, dres_s_dsdc, dres_c_dsdc, dummy, dummy); + } + + void TransportModelPolymer::ResidualEquation::computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c, + const bool if_dres_s_dsdc, const bool if_dres_c_dsdc, + double* res, double* dres_s_dsdc, + double* dres_c_dsdc, double& mc, double& ff) const + { + if ((if_dres_s_dsdc || if_dres_c_dsdc) && gradient_method == Analytic) { + double s = x[0]; + double c = x[1]; + std::vector dff_dsdc(2); + double mc_dc; + double ads_dc; + double ads; + tm.fracFlowWithDer(s, c, cmax0, cell, ff, dff_dsdc); + if (if_dres_c_dsdc) { + tm.polyprops_.adsorptionWithDer(c, cmax0, ads, ads_dc); + tm.computeMcWithDer(c, mc, mc_dc); + } else { + tm.polyprops_.adsorption(c, cmax0, ads); + tm.computeMc(c, mc); + } + if (if_res_s) { + res[0] = s - s0 + dtpv*(outflux*ff + influx); + } + if (if_res_c) { + res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0 + + rhor*((1.0 - porosity)/porosity)*(ads - ads0) + + dtpv*(outflux*ff*mc + influx_polymer); + } + if (if_dres_s_dsdc) { + dres_s_dsdc[0] = 1 + dtpv*outflux*dff_dsdc[0]; + dres_s_dsdc[1] = dtpv*outflux*dff_dsdc[1]; + } + if (if_dres_c_dsdc) { + dres_c_dsdc[0] = (1 - dps)*c + dtpv*outflux*(dff_dsdc[0])*mc; + dres_c_dsdc[1] = (1 - dps)*s + rhor*((1.0 - porosity)/porosity)*ads_dc + + dtpv*outflux*(dff_dsdc[1]*mc + ff*mc_dc); + } + + } else if (if_res_c || if_res_s) { + double s = x[0]; + double c = x[1]; + tm.fracFlow(s, c, cmax0, cell, ff); + tm.computeMc(c, mc); + double ads; + tm.polyprops_.adsorption(c, cmax0, ads); + if (if_res_s) { + res[0] = s - s0 + dtpv*(outflux*ff + influx); + } + if (if_res_c) { + res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0 + + rhor*((1.0 - porosity)/porosity)*(ads - ads0) + + dtpv*(outflux*ff*mc + influx_polymer); + } + } + + if ((if_dres_c_dsdc || if_dres_s_dsdc) && gradient_method == FinDif) { double epsi = 1e-8; - double res[2]; double res_epsi[2]; double x_epsi[2]; computeResidual(x, res); - x_epsi[0] = x[0] + epsi; - x_epsi[1] = x[1]; - computeResidual(x_epsi, res_epsi); - res_s_ds_dc[0] = (res_epsi[0] - res[0])/epsi; - x_epsi[0] = x[0]; - x_epsi[1] = x[1] + epsi; - computeResidual(x_epsi, res_epsi); - res_s_ds_dc[1] = (res_epsi[0] - res[0])/epsi; - x_epsi[0] = x[0] + epsi; - x_epsi[1] = x[1]; - computeResidual(x_epsi, res_epsi); - res_c_ds_dc[0] = (res_epsi[1] - res[1])/epsi; - x_epsi[0] = x[0]; - x_epsi[1] = x[1] + epsi; - computeResidual(x_epsi, res_epsi); - res_c_ds_dc[1] = (res_epsi[1] - res[1])/epsi; - } else if (gradient_method == Analytic) { - double s = x[0]; - double c = x[1]; - double ff_ds_dc[2]; - double ff = tm.fracFlowWithDer(s, c, cell, ff_ds_dc); - double mc_dc; - double mc = tm.computeMcWithDer(c, &mc_dc); - double ads_dc; - tm.polyprops_.adsorptionWithDer(c, cmax0, &ads_dc); - res_s_ds_dc[0] = 1 + dtpv*outflux*ff_ds_dc[0]; - res_s_ds_dc[1] = dtpv*outflux*ff_ds_dc[1]; - res_c_ds_dc[0] = (1 - dps)*c + dtpv*outflux*(ff_ds_dc[0])*mc; - res_c_ds_dc[1] = (1 - dps)*s + rhor*((1.0 - porosity)/porosity)*ads_dc - + dtpv*outflux*(ff_ds_dc[1]*mc + ff*mc_dc); - } + if (if_dres_s_dsdc) { + x_epsi[0] = x[0] + epsi; + x_epsi[1] = x[1]; + computeResidual(x_epsi, res_epsi); + dres_s_dsdc[0] = (res_epsi[0] - res[0])/epsi; + x_epsi[0] = x[0]; + x_epsi[1] = x[1] + epsi; + computeResidual(x_epsi, res_epsi); + dres_s_dsdc[1] = (res_epsi[0] - res[0])/epsi; + } + if (if_dres_c_dsdc) { + x_epsi[0] = x[0] + epsi; + x_epsi[1] = x[1]; + computeResidual(x_epsi, res_epsi); + dres_c_dsdc[0] = (res_epsi[1] - res[1])/epsi; + x_epsi[0] = x[0]; + x_epsi[1] = x[1] + epsi; + computeResidual(x_epsi, res_epsi); + dres_c_dsdc[1] = (res_epsi[1] - res[1])/epsi; + } + } } - void TransportModelPolymer::solveSingleCell(const int cell) { switch (method_) { @@ -596,8 +578,9 @@ namespace Opm concentration_[cell] = modifiedRegulaFalsi(res, a, b, maxit_, tol_, iters_used); cmax_[cell] = std::max(cmax_[cell], concentration_[cell]); saturation_[cell] = res.lastSaturation(); - fractionalflow_[cell] = fracFlow(saturation_[cell], concentration_[cell], cell); - mc_[cell] = computeMc(concentration_[cell]); + fracFlow(saturation_[cell], concentration_[cell], cmax_[cell], cell, + fractionalflow_[cell]); + computeMc(concentration_[cell], mc_[cell]); } @@ -797,8 +780,9 @@ namespace Opm // Must set initial fractional flows etc. before we start. for (int i = 0; i < num_cells; ++i) { const int cell = cells[i]; - fractionalflow_[cell] = fracFlow(saturation_[cell], concentration_[cell], cell); - mc_[cell] = computeMc(concentration_[cell]); + fracFlow(saturation_[cell], concentration_[cell], cmax_[cell], + cell, fractionalflow_[cell]); + computeMc(concentration_[cell], mc_[cell]); s0[i] = saturation_[cell]; c0[i] = concentration_[cell]; cmax0[i] = cmax_[i]; @@ -842,115 +826,55 @@ namespace Opm } - - - double TransportModelPolymer::fracFlow(double s, double c, int cell) const + void TransportModelPolymer::fracFlow(double s, double c, double cmax, + int cell, double& ff) const { - double c_max_limit = polyprops_.cMax(); - double cbar = c/c_max_limit; - double c_ads = polyprops_.adsorption(c, cmax_[cell]); - double c_max_ads = polyprops_.cMaxAds(); - double res_factor = polyprops_.resFactor(); - double res_k = 1 + (res_factor - 1)*c_ads/c_max_ads; - double mu_w = visc_[0]; - double mu_m = polyprops_.viscMult(c)*mu_w; - double mu_p = polyprops_.viscMult(polyprops_.cMax())*mu_w; - double omega = polyprops_.mixParam(); - double mu_m_omega = std::pow(mu_m, omega); - double mu_w_e = mu_m_omega*std::pow(mu_w, 1.0 - omega); - double mu_p_eff = mu_m_omega*std::pow(mu_p, 1.0 - omega); - double inv_mu_w_eff = (1.0 - cbar)/mu_w_e + cbar/mu_p_eff; - double inv_visc_eff[2] = { inv_mu_w_eff, 1.0/visc_[1] }; - double sat[2] = { s, 1.0 - s }; - double mob[2]; - props_.relperm(1, sat, &cell, mob, 0); - mob[0] *= inv_visc_eff[0]/res_k; - mob[1] *= inv_visc_eff[1]; - return mob[0]/(mob[0] + mob[1]); + std::vector dummy; + fracFlowBoth(s, c, cmax, cell, ff, dummy, false); + } + + void TransportModelPolymer::fracFlowWithDer(double s, double c, double cmax, + int cell, double& ff, + std::vector& dff_dsdc) const + { + fracFlowBoth(s, c, cmax, cell, ff, dff_dsdc, true); } - double TransportModelPolymer::fracFlowWithDer(double s, double c, int cell, double* der) const + void TransportModelPolymer::fracFlowBoth(double s, double c, double cmax, int cell, + double& ff, std::vector& dff_dsdc, + bool if_with_der) const { - double c_max_limit = polyprops_.cMax(); - double cbar = c/c_max_limit; - double c_ads_dc; - double c_max = cmax_[cell]; - double c_ads = polyprops_.adsorptionWithDer(c, c_max, &c_ads_dc); - double c_max_ads = polyprops_.cMaxAds(); - double res_factor = polyprops_.resFactor(); - double res_k = 1 + (res_factor - 1)*c_ads/c_max_ads; - double res_k_dc = (res_factor - 1)*c_ads_dc/c_max_ads; - double mu_w = visc_[0]; - double mu_m_dc; // derivative of mu_m with respect to c - double mu_m = polyprops_.viscMultWithDer(c, &mu_m_dc)*mu_w; - mu_m_dc *= mu_w; - double mu_p = polyprops_.viscMult(c_max_limit)*mu_w; - double omega = polyprops_.mixParam(); - double mu_w_e = std::pow(mu_m, omega)*std::pow(mu_w, 1 - omega); - double mu_w_e_dc = omega*mu_m_dc*std::pow(mu_m, omega - 1)*std::pow(mu_w, 1 - omega); - double mu_p_eff = std::pow(mu_m, omega)*std::pow(mu_p, 1 - omega); - double mu_p_eff_dc = omega*mu_m_dc*std::pow(mu_m, omega - 1)*std::pow(mu_p, 1 - omega); - double mu_w_eff = 1./((1 - cbar)/mu_w_e + cbar/mu_p_eff); - double mu_w_eff_dc = -1./c_max_limit*mu_w_eff*mu_w_eff*(1./mu_p_eff - 1./mu_w_e) - + (1-cbar)*(mu_w_eff*mu_w_eff/(mu_w_e*mu_w_e))*mu_w_e_dc - + cbar*(mu_w_eff*mu_w_eff/(mu_p_eff*mu_p_eff))*mu_p_eff_dc; - double visc_eff[2] = { mu_w_eff, visc_[1] }; - double sat[2] = { s, 1.0 - s }; - double mob[2]; - double mob_ds[2]; - double mob_dc[2]; - double perm[2]; - double perm_ds[4]; - props_.relperm(1, sat, &cell, perm, perm_ds); - mob[0] = perm[0]/visc_eff[0]/res_k; - mob[1] = perm[1]/visc_eff[1]; - mob_ds[0] = perm_ds[0]/visc_eff[0]/res_k; - mob_ds[1] = perm_ds[1]/visc_eff[1]; - mob_dc[0] = - perm[0]*(mu_w_eff_dc/(mu_w_eff*mu_w_eff*res_k) + res_k_dc/(res_k*res_k*mu_w_eff)); - mob_dc[1] = 0.; - der[0] = (mob_ds[0]*mob[1] - mob_ds[1]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); - der[1] = (mob_dc[0]*mob[1] - mob_dc[1]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); - return mob[0]/(mob[0] + mob[1]); + double relperm[2]; + double drelperm_ds[4]; + double sat[2] = {s, 1 - s}; + props_.relperm(1, sat, &cell, relperm, drelperm_ds); + std::vector mob(2); + std::vector dmob_ds(2); + std::vector dmob_dc(2); + double dmobwat_dc; + polyprops_.effectiveMobilitiesBoth(c, cmax, visc_, relperm, drelperm_ds, + mob, dmob_ds, dmobwat_dc, if_with_der); + dmob_dc[0] = dmobwat_dc; + dmob_dc[1] = 0.; + ff = mob[0]/(mob[0] + mob[1]); + if (if_with_der) { + dff_dsdc[0] = (dmob_ds[0]*mob[1] - dmob_ds[1]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); + dff_dsdc[1] = (dmob_dc[0]*mob[1] - dmob_dc[1]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); + } else { + dff_dsdc.clear(); + } } - double TransportModelPolymer::computeMc(double c) const + void TransportModelPolymer::computeMc(double c, double& mc) const { - double c_max_limit = polyprops_.cMax(); - double cbar = c/c_max_limit; - double mu_w = visc_[0]; - double mu_m = polyprops_.viscMult(c)*mu_w; - double mu_p = polyprops_.viscMult(polyprops_.cMax())*mu_w; - double omega = polyprops_.mixParam(); - double mu_m_omega = std::pow(mu_m, omega); - double mu_w_e = mu_m_omega*std::pow(mu_w, 1.0 - omega); - double mu_p_eff = mu_m_omega*std::pow(mu_p, 1.0 - omega); - double inv_mu_w_eff = (1.0 - cbar)/mu_w_e + cbar/mu_p_eff; - return c/(inv_mu_w_eff*mu_p_eff); + double dummy; + polyprops_.computeMcBoth(c, mc, dummy, false); } - double TransportModelPolymer::computeMcWithDer(double c, double* der) const + void TransportModelPolymer::computeMcWithDer(double c, double& mc, + double &dmc_dc) const { - double c_max_limit = polyprops_.cMax(); - double cbar = c/c_max_limit; - double mu_w = visc_[0]; - double mu_m_dc; // derivative of mu_m with respect to c - double mu_m = polyprops_.viscMultWithDer(c, &mu_m_dc)*mu_w; - mu_m_dc *= mu_w; - double mu_p = polyprops_.viscMult(polyprops_.cMax())*mu_w; - double omega = polyprops_.mixParam(); - double mu_m_omega = std::pow(mu_m, omega); - double mu_m_omega_minus1 = std::pow(mu_m, omega-1); - double mu_w_omega = std::pow(mu_w, 1.0 - omega); - double mu_w_e = mu_m_omega*mu_w_omega; - double mu_w_e_dc = omega*mu_m_dc*mu_m_omega_minus1*mu_w_omega; - double mu_p_omega = std::pow(mu_p, 1.0 - omega); - double mu_p_eff = mu_m_omega*mu_p_omega; - double mu_p_eff_dc = omega*mu_m_dc*mu_m_omega_minus1*mu_p_omega; - double mu_w_eff = 1./((1 - cbar)/mu_w_e + cbar/mu_p_eff); - double inv_mu_w_eff_dc = -mu_w_e_dc/(mu_w_e*mu_w_e)*(1. - cbar) - mu_p_eff_dc/(mu_p_eff*mu_p_eff)*cbar + (1./mu_p_eff - 1./mu_w_e); - double mu_w_eff_dc = -mu_w_eff*mu_w_eff*inv_mu_w_eff_dc; - *der = mu_w_eff/mu_p_eff + c*mu_w_eff_dc/mu_p_eff - c*mu_p_eff_dc*mu_w_eff/(mu_p_eff*mu_p_eff); - return c*mu_w_eff/mu_p_eff; + polyprops_.computeMcBoth(c, mc, dmc_dc, true); } } // namespace Opm @@ -1084,17 +1008,17 @@ namespace bool solveNewtonStep(const double* x, const Opm::TransportModelPolymer::ResidualEquation& res_eq, const double* res, double* x_new) { - double res_s_ds_dc[2]; - double res_c_ds_dc[2]; + double dres_s_dsdc[2]; + double dres_c_dsdc[2]; - res_eq.computeJacobiRes(x, res_s_ds_dc, res_c_ds_dc); + res_eq.computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc); - double det = res_s_ds_dc[0]*res_c_ds_dc[1] - res_c_ds_dc[0]*res_s_ds_dc[1]; + double det = dres_s_dsdc[0]*dres_c_dsdc[1] - dres_c_dsdc[0]*dres_s_dsdc[1]; if (std::abs(det) < 1e-8) { return false; } else { - x_new[0] = x[0] - (res[0]*res_c_ds_dc[1] - res[1]*res_s_ds_dc[1])/det; - x_new[1] = x[1] - (res[1]*res_s_ds_dc[0] - res[0]*res_c_ds_dc[0])/det; + x_new[0] = x[0] - (res[0]*dres_c_dsdc[1] - res[1]*dres_s_dsdc[1])/det; + x_new[1] = x[1] - (res[1]*dres_s_dsdc[0] - res[0]*dres_c_dsdc[0])/det; return true; } } diff --git a/opm/polymer/TransportModelPolymer.hpp b/opm/polymer/TransportModelPolymer.hpp index 4e9261aa5..d7ed11a51 100644 --- a/opm/polymer/TransportModelPolymer.hpp +++ b/opm/polymer/TransportModelPolymer.hpp @@ -97,10 +97,13 @@ namespace Opm struct ResidualS; - double fracFlow(double s, double c, int cell) const; - double fracFlowWithDer(double s, double c, int cell, double* der) const; - double computeMc(double c) const; - double computeMcWithDer(double c, double* der) const; + void fracFlow(double s, double c, double cmax, int cell, double& ff) const; + void fracFlowWithDer(double s, double cmax, double c, int cell, double& ff, + std::vector& dff_dsdc) const; + void fracFlowBoth(double s, double c, double cmax, int cell, double& ff, + std::vector& dff_dsdc, bool if_with_der) const; + void computeMc(double c, double& mc) const; + void computeMcWithDer(double c, double& mc, double& dmc_dc) const; }; } // namespace Opm diff --git a/opm/polymer/polymerUtilities.cpp b/opm/polymer/polymerUtilities.cpp index 92b54aaf8..8cf367ced 100644 --- a/opm/polymer/polymerUtilities.cpp +++ b/opm/polymer/polymerUtilities.cpp @@ -204,7 +204,9 @@ namespace Opm const double* poro = props.porosity(); double abs_mass = 0.0; for (int cell = 0; cell < num_cells; ++cell) { - abs_mass += polyprops.simpleAdsorption(cmax[cell])*pv[cell]*((1.0 - poro[cell])/poro[cell])*rhor; + double c_ads; + polyprops.simpleAdsorption(cmax[cell], c_ads); + abs_mass += c_ads*pv[cell]*((1.0 - poro[cell])/poro[cell])*rhor; } return abs_mass; }