OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #652 from GitPaean/buggy_pointer_usage

Browse Source 
correcting buggy pointer usage for visc parameters in effective viscosity calcuation
This commit is contained in:
Atgeirr Flø Rasmussen
2016-04-14 15:42:00 +02:00
parent c64593a023 c120313f09
commit fd99a9e6e9
7 changed files with 115 additions and 113 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
opm/polymer/CompressibleTpfaPolymer.cpp
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]);
}

36
opm/polymer/PolymerProperties.cpp
View File

@@ -192,37 +192,36 @@ namespace Opm
}
void PolymerProperties::effectiveVisc(const double c, const double* visc, double& mu_w_eff) const {
effectiveInvVisc(c, visc, mu_w_eff);
void PolymerProperties::effectiveVisc(const double c, const double mu_w, double& mu_w_eff) const {
effectiveInvVisc(c, mu_w, 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 {
effectiveInvViscWithDer(c, visc, mu_w_eff, dmu_w_eff_dc);
void PolymerProperties::effectiveViscWithDer(const double c, const double mu_w, double& mu_w_eff, double dmu_w_eff_dc) const {
effectiveInvViscWithDer(c, mu_w, 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 mu_w, double& inv_mu_w_eff) const
{
double dummy;
effectiveInvViscBoth(c, visc, inv_mu_w_eff, dummy, false);
effectiveInvViscBoth(c, mu_w, inv_mu_w_eff, dummy, false);
}
void PolymerProperties::effectiveInvViscWithDer(const double c, const double* visc,
void PolymerProperties::effectiveInvViscWithDer(const double c, const double mu_w,
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);
effectiveInvViscBoth(c, mu_w, 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 mu_w,
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_;
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,31 +243,30 @@ namespace Opm
}
void PolymerProperties::effectiveInvPolyVisc(const double c,
const double* visc,
const double mu_w,
double& inv_mu_p_eff) const
{
double dummy;
effectiveInvPolyViscBoth(c, visc, inv_mu_p_eff, dummy, false);
effectiveInvPolyViscBoth(c, mu_w, inv_mu_p_eff, dummy, false);
}
void PolymerProperties::effectiveInvPolyViscWithDer(const double c,
const double* visc,
const double mu_w,
double& inv_mu_p_eff,
double& d_inv_mu_p_eff_dc) const
{
effectiveInvPolyViscBoth(c, visc, inv_mu_p_eff, d_inv_mu_p_eff_dc, true);
effectiveInvPolyViscBoth(c, mu_w, inv_mu_p_eff, d_inv_mu_p_eff_dc, true);
}
void PolymerProperties::effectiveInvPolyViscBoth(const double c,
const double* visc,
const double mu_w,
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];
double mu_m = 0.0;
double dmu_m_dc = 0.0;
@@ -373,7 +371,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;

22
opm/polymer/PolymerProperties.hpp
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 mu_w,
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 mu_w,
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 mu_w,
double& inv_mu_w_eff,
double& dinv_mu_w_eff_dc) const;
void effectiveInvPolyVisc(const double c,
const double* visc,
const double mu_w,
double& inv_mu_p_eff) const;
void effectiveInvPolyViscWithDer(const double c,
const double* visc,
const double mu_w,
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 mu_w,
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 mu_w,
double& inv_mu_p_eff,
double& dinv_mu_p_eff_dc,
const bool if_with_der) const;

6
opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
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();

2
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
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_);

56
opm/polymer/fullyimplicit/PolymerPropsAd.cpp
View File

@@ -54,11 +54,11 @@ namespace Opm {
double
PolymerPropsAd::cMax() const
{
return polymer_props_.cMax();
}
double
PolymerPropsAd::cMax() const
{
return polymer_props_.cMax();
}
const std::vector<double>&
PolymerPropsAd::shearWaterVelocity() const
@@ -130,7 +130,7 @@ namespace Opm {
}
PolymerPropsAd::PolymerPropsAd(const PolymerProperties& polymer_props)
: polymer_props_ (polymer_props)
@@ -150,13 +150,14 @@ namespace Opm {
V PolymerPropsAd::effectiveInvWaterVisc(const V& c,
const double* visc) const
const V& mu_w) const
{
assert(c.size() == mu_w.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), mu_w(i), im);
inv_mu_w_eff(i) = im;
}
@@ -168,17 +169,18 @@ namespace Opm {
ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c,
const double* visc) const
const V& mu_w) const
{
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);
inv_mu_w_eff(i) = im;
dinv_mu_w_eff(i) = dim;
}
assert(c.size() == mu_w.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), mu_w(i), im, dim);
inv_mu_w_eff(i) = im;
dinv_mu_w_eff(i) = dim;
}
ADB::M dim_diag(dinv_mu_w_eff.matrix().asDiagonal());
const int num_blocks = c.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
@@ -192,16 +194,16 @@ namespace Opm {
ADB PolymerPropsAd::effectiveInvPolymerVisc(const ADB& c, const double* visc) const
ADB PolymerPropsAd::effectiveInvPolymerVisc(const ADB& c, const V& mu_w) const
{
assert(c.size() == mu_w.size());
const int nc = c.size();
V inv_mu_p_eff(nc);
V dinv_mu_p_eff(nc);
for (int i = 0; i < nc; ++i) {
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), mu_w(i), im, dim);
inv_mu_p_eff(i) = im;
dinv_mu_p_eff(i) = dim;
}
@@ -230,7 +232,7 @@ namespace Opm {
polymer_props_.computeMc(c(i), m);
mc(i) = m;
}
return mc;
}
@@ -240,11 +242,11 @@ namespace Opm {
ADB PolymerPropsAd::polymerWaterVelocityRatio(const ADB& c) const
{
const int nc = c.size();
V mc(nc);
V dmc(nc);
for (int i = 0; i < nc; ++i) {
double m = 0;
double dm = 0;
@@ -317,7 +319,7 @@ namespace Opm {
V
PolymerPropsAd::effectiveRelPerm(const V& c,
PolymerPropsAd::effectiveRelPerm(const V& c,
const V& cmax_cells,
const V& krw) const
{
@@ -350,8 +352,8 @@ namespace Opm {
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 rk = one + ads * factor;
return krw / rk;
}

104
opm/polymer/fullyimplicit/PolymerPropsAd.hpp
View File

@@ -28,18 +28,18 @@
#include <opm/polymer/PolymerProperties.hpp>
namespace Opm {
class PolymerPropsAd
class PolymerPropsAd
{
public:
/// \return Reference rock density.
/// \return Reference rock density.
double rockDensity() const;
/// \return The value of dead pore volume.
/// \return The value of dead pore volume.
double deadPoreVol() const;
/// \return The max concentration injected.
double cMax() const;
/// \return The max concentration injected.
double cMax() const;
/// \ return The water velcoity or shear rate in the PLYSHLOG table
const std::vector<double>& shearWaterVelocity() const;
@@ -67,75 +67,75 @@ namespace Opm {
double viscMult(double c) const; // multipler interpolated from PLYVISC table
typedef AutoDiffBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
V viscMult(const V& c) const;
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n viscosity multiplier from PLVISC table.
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n viscosity multiplier from PLVISC table.
/// Constructor wrapping a polymer props.
/// Constructor wrapping a polymer props.
PolymerPropsAd(const PolymerProperties& polymer_props);
/// Destructor.
/// Destructor.
~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.
V
effectiveInvWaterVisc(const V& c,const double* 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.
ADB
effectiveInvWaterVisc(const ADB& c,const double* visc) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] mu_w Array of n water viscosity values.
/// \return Array of inverse effective water viscosity.
V
effectiveInvWaterVisc(const V& c, const V& mu_w) const;
/// \param[in] c ADB of polymer concentraion values.
/// \param[in] visc Array of water viscosity values
/// \return ADB of inverse effective polymer viscosity.
/// \param[in] c ADB of polymer concentraion.
/// \param[in] mu_w Array of water viscosity value.
/// \return ADB of inverse effective water viscosity.
ADB
effectiveInvPolymerVisc(const ADB& c, const double* visc) const;
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n mc values, here mc means m(c) * c.
V
effectiveInvWaterVisc(const ADB& c,const V& mu_w) const;
/// \param[in] c ADB of polymer concentraion values.
/// \param[in] mu_w Array of water viscosity values
/// \return ADB of inverse effective polymer viscosity.
ADB
effectiveInvPolymerVisc(const ADB& c, const V& mu_w) const;
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n mc values, here mc means m(c) * c.
V
polymerWaterVelocityRatio(const V& c) const;
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n mc values, here mc means m(c) * c.
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n mc values, here mc means m(c) * c.
ADB
polymerWaterVelocityRatio(const ADB& c) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \return Array of n adsorption values.
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \return Array of n adsorption values.
V
adsorption(const V& c, const V& cmax_cells) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \return Array of n adsorption values.
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \return Array of n adsorption values.
ADB
adsorption(const ADB& c, const ADB& cmax_cells) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \param[in] relperm Array of n relative water relperm values.
/// \return Array of n adsorption values.
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \param[in] relperm Array of n relative water relperm values.
/// \return Array of n adsorption values.
V
effectiveRelPerm(const V& c, const V& cmax_cells, const V& relperm) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \param[in] relperm Array of n relative water relperm values.
/// \return Array of n adsorption values.
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \param[in] relperm Array of n relative water relperm values.
/// \return Array of n adsorption values.
ADB
effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw) const;
@@ -150,7 +150,7 @@ namespace Opm {
private:
const PolymerProperties& polymer_props_;
};
} //namespace Opm
#endif// OPM_POLYMERPROPSAD_HEADED_INLCUDED