correcting buggy visc pointer usage in effective viscosity calculation

This commit is contained in:
Kai Bao 2016-04-14 14:41:54 +02:00
parent c64593a023
commit 30774d8309
7 changed files with 47 additions and 42 deletions

View File

@ -140,7 +140,7 @@ namespace Opm
props_.viscosity(nc, cell_p, cell_T, cell_z, &allcells_[0], &cell_viscosity_[0], 0);
cell_phasemob_.resize(nc*np);
for (int cell = 0; cell < nc; ++cell) {
poly_props_.effectiveVisc((*c_)[cell], &cell_viscosity_[np*cell + 0], cell_eff_viscosity_[np*cell + 0]);
poly_props_.effectiveVisc((*c_)[cell], cell_viscosity_[np*cell + 0], cell_eff_viscosity_[np*cell + 0]);
poly_props_.effectiveMobilities((*c_)[cell], (*cmax_)[cell], &cell_viscosity_[np*cell + 0], &cell_relperm_[np*cell + 0], &cell_phasemob_[np*cell + 0]);
}

View File

@ -192,37 +192,37 @@ namespace Opm
}
void PolymerProperties::effectiveVisc(const double c, const double* visc, double& mu_w_eff) const {
void PolymerProperties::effectiveVisc(const double c, const double visc, double& mu_w_eff) const {
effectiveInvVisc(c, visc, mu_w_eff);
mu_w_eff = 1./mu_w_eff;
}
void PolymerProperties::effectiveViscWithDer(const double c, const double* visc, double& mu_w_eff, double dmu_w_eff_dc) const {
void PolymerProperties::effectiveViscWithDer(const double c, const double visc, double& mu_w_eff, double dmu_w_eff_dc) const {
effectiveInvViscWithDer(c, visc, mu_w_eff, dmu_w_eff_dc);
mu_w_eff = 1./mu_w_eff;
dmu_w_eff_dc = -dmu_w_eff_dc*mu_w_eff*mu_w_eff;
}
void PolymerProperties::effectiveInvVisc(const double c, const double* visc, double& inv_mu_w_eff) const
void PolymerProperties::effectiveInvVisc(const double c, const double visc, double& inv_mu_w_eff) const
{
double dummy;
effectiveInvViscBoth(c, visc, inv_mu_w_eff, dummy, false);
}
void PolymerProperties::effectiveInvViscWithDer(const double c, const double* visc,
void PolymerProperties::effectiveInvViscWithDer(const double c, const double visc,
double& inv_mu_w_eff,
double& dinv_mu_w_eff_dc) const {
effectiveInvViscBoth(c, visc, inv_mu_w_eff, dinv_mu_w_eff_dc, true);
}
void PolymerProperties::effectiveInvViscBoth(const double c, const double* visc,
void PolymerProperties::effectiveInvViscBoth(const double c, const double visc,
double& inv_mu_w_eff,
double& dinv_mu_w_eff_dc,
bool if_with_der) const {
double cbar = c/c_max_;
double mu_w = visc[0];
const double cbar = c/c_max_;
const double mu_w = visc;
double mu_m;
double omega = mix_param_;
const double omega = mix_param_;
double dmu_m_dc;
if (if_with_der) {
mu_m = viscMultWithDer(c, &dmu_m_dc)*mu_w;
@ -244,7 +244,7 @@ namespace Opm
}
void PolymerProperties::effectiveInvPolyVisc(const double c,
const double* visc,
const double visc,
double& inv_mu_p_eff) const
{
double dummy;
@ -253,7 +253,7 @@ namespace Opm
}
void PolymerProperties::effectiveInvPolyViscWithDer(const double c,
const double* visc,
const double visc,
double& inv_mu_p_eff,
double& d_inv_mu_p_eff_dc) const
{
@ -262,13 +262,13 @@ namespace Opm
}
void PolymerProperties::effectiveInvPolyViscBoth(const double c,
const double* visc,
const double visc,
double& inv_mu_p_eff,
double& dinv_mu_p_eff_dc,
const bool if_with_der) const
{
const double omega = mix_param_;
const double mu_w = visc[0];
const double mu_w = visc;
double mu_m = 0.0;
double dmu_m_dc = 0.0;
@ -373,7 +373,7 @@ namespace Opm
{
double inv_mu_w_eff;
double dinv_mu_w_eff_dc;
effectiveInvViscBoth(c, visc, inv_mu_w_eff, dinv_mu_w_eff_dc, if_with_der);
effectiveInvViscBoth(c, visc[0], inv_mu_w_eff, dinv_mu_w_eff_dc, if_with_der);
double eff_relperm_wat;
double deff_relperm_wat_ds;
double deff_relperm_wat_dc;

View File

@ -278,26 +278,27 @@ namespace Opm
void adsorptionWithDer(double c, double cmax,
double& c_ads, double& dc_ads_dc) const;
void effectiveVisc(const double c, const double* visc,
void effectiveVisc(const double c, const double visc,
double& mu_w_eff) const;
void effectiveViscWithDer(const double c, const double* visc
void effectiveViscWithDer(const double c, const double visc
, double& mu_w_eff
, double dmu_w_eff_dc) const;
void effectiveInvVisc(const double c, const double* visc,
void effectiveInvVisc(const double c, const double visc,
double& inv_mu_w_eff) const;
void effectiveInvViscWithDer(const double c,
const double* visc,
double& inv_mu_w_eff,
double& dinv_mu_w_eff_dc) const;
const double visc,
double& inv_mu_w_eff,
double& dinv_mu_w_eff_dc) const;
void effectiveInvPolyVisc(const double c,
const double* visc,
const double visc,
double& inv_mu_p_eff) const;
void effectiveInvPolyViscWithDer(const double c,
const double* visc,
const double visc,
double& inv_mu_p_eff,
double& d_inv_mu_p_eff_dc) const;
@ -407,12 +408,13 @@ namespace Opm
void adsorptionBoth(double c, double cmax,
double& c_ads, double& dc_ads_dc,
bool if_with_der) const;
void effectiveInvViscBoth(const double c, const double* visc,
void effectiveInvViscBoth(const double c, const double visc,
double& inv_mu_w_eff,
double& dinv_mu_w_eff_dc, bool if_with_der) const;
void effectiveInvPolyViscBoth(const double c,
const double* visc,
const double visc,
double& inv_mu_p_eff,
double& dinv_mu_p_eff_dc,
const bool if_with_der) const;

View File

@ -460,11 +460,11 @@ namespace Opm {
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
const ADB mc = computeMc(state);
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr);
const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value().data());
const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
// Reduce mobility of water phase by relperm reduction and effective viscosity increase.
rq_[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc;
// Compute polymer mobility.
const ADB inv_poly_eff_visc = polymer_props_ad_.effectiveInvPolymerVisc(state.concentration, mu.value().data());
const ADB inv_poly_eff_visc = polymer_props_ad_.effectiveInvPolymerVisc(state.concentration, mu.value());
rq_[poly_pos_].mob = tr_mult * mc * krw_eff * inv_poly_eff_visc;
rq_[poly_pos_].b = rq_[actph].b;
rq_[poly_pos_].dh = rq_[actph].dh;
@ -624,7 +624,7 @@ namespace Opm {
ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration,
cmax,
kr[canonicalPhaseIdx]);
ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value().data());
ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
rq_[ phase ].mob = tr_mult * krw_eff * inv_wat_eff_visc;
const V& polymer_conc = state.concentration.value();

View File

@ -871,7 +871,7 @@ namespace {
const ADB& temp = state.temperature;
const ADB mu_w = fluidViscosity(0, press[0], temp, cond, cells_);
ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu_w.value().data());
ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu_w.value());
rq_[0].mob = tr_mult * krw_eff * inv_wat_eff_vis;
rq_[2].mob = tr_mult * mc * krw_eff * inv_wat_eff_vis;
const ADB mu_o = fluidViscosity(1, press[1], temp, cond, cells_);

View File

@ -150,13 +150,14 @@ namespace Opm {
V PolymerPropsAd::effectiveInvWaterVisc(const V& c,
const double* visc) const
const V& visc) const
{
assert(c.size() == visc.size());
const int nc = c.size();
V inv_mu_w_eff(nc);
for (int i = 0; i < nc; ++i) {
double im = 0;
polymer_props_.effectiveInvVisc(c(i), visc, im);
polymer_props_.effectiveInvVisc(c(i), visc(i), im);
inv_mu_w_eff(i) = im;
}
@ -168,14 +169,15 @@ namespace Opm {
ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c,
const double* visc) const
const V& visc) const
{
assert(c.size() == visc.size());
const int nc = c.size();
V inv_mu_w_eff(nc);
V dinv_mu_w_eff(nc);
for (int i = 0; i < nc; ++i) {
double im = 0, dim = 0;
polymer_props_.effectiveInvViscWithDer(c.value()(i), visc, im, dim);
polymer_props_.effectiveInvViscWithDer(c.value()(i), visc(i), im, dim);
inv_mu_w_eff(i) = im;
dinv_mu_w_eff(i) = dim;
}
@ -192,8 +194,9 @@ namespace Opm {
ADB PolymerPropsAd::effectiveInvPolymerVisc(const ADB& c, const double* visc) const
ADB PolymerPropsAd::effectiveInvPolymerVisc(const ADB& c, const V& visc) const
{
assert(c.size() == visc.size());
const int nc = c.size();
V inv_mu_p_eff(nc);
V dinv_mu_p_eff(nc);
@ -201,7 +204,7 @@ namespace Opm {
double im = 0;
double dim = 0;
// TODO: the usage of visc can be likely wrong, while more investigation will be requried.
polymer_props_.effectiveInvPolyViscWithDer(c.value()(i), visc, im, dim);
polymer_props_.effectiveInvPolyViscWithDer(c.value()(i), visc(i), im, dim);
inv_mu_p_eff(i) = im;
dinv_mu_p_eff(i) = dim;
}

View File

@ -81,22 +81,22 @@ namespace Opm {
~PolymerPropsAd();
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] visc Array of 2 viscosity value.
/// \return value of inverse effective water viscosity.
/// \param[in] visc Array of n water viscosity values.
/// \return Array of inverse effective water viscosity.
V
effectiveInvWaterVisc(const V& c,const double* visc) const;
effectiveInvWaterVisc(const V& c, const V& visc) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] visc Array of 2 viscosity value
/// \return value of inverse effective water viscosity.
/// \param[in] c ADB of polymer concentraion.
/// \param[in] visc Array of water viscosity value.
/// \return ADB of inverse effective water viscosity.
ADB
effectiveInvWaterVisc(const ADB& c,const double* visc) const;
effectiveInvWaterVisc(const ADB& c,const V& visc) const;
/// \param[in] c ADB of polymer concentraion values.
/// \param[in] visc Array of water viscosity values
/// \return ADB of inverse effective polymer viscosity.
ADB
effectiveInvPolymerVisc(const ADB& c, const double* visc) const;
effectiveInvPolymerVisc(const ADB& c, const V& visc) const;
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n mc values, here mc means m(c) * c.