Merge pull request #227 from totto82/fix_updateState3

Fixes in updateState
This commit is contained in:
Atgeirr Flø Rasmussen 2014-11-10 09:24:23 +01:00
commit 70f390f705
2 changed files with 17 additions and 23 deletions

View File

@ -60,7 +60,7 @@ namespace Opm {
{
double dp_max_rel_;
double ds_max_;
double drs_max_rel_;
double dr_max_rel_;
enum RelaxType relax_type_;
double relax_max_;
double relax_increment_;
@ -338,7 +338,7 @@ namespace Opm {
double dpMaxRel() const { return param_.dp_max_rel_; }
double dsMax() const { return param_.ds_max_; }
double drsMaxRel() const { return param_.drs_max_rel_; }
double drMaxRel() const { return param_.dr_max_rel_; }
enum RelaxType relaxType() const { return param_.relax_type_; }
double relaxMax() const { return param_.relax_max_; };
double relaxIncrement() const { return param_.relax_increment_; };

View File

@ -163,7 +163,7 @@ namespace {
// default values for the solver parameters
dp_max_rel_ = 1.0e9;
ds_max_ = 0.2;
drs_max_rel_ = 1.0e9;
dr_max_rel_ = 1.0e9;
relax_type_ = DAMPEN;
relax_max_ = 0.5;
relax_increment_ = 0.1;
@ -189,7 +189,7 @@ namespace {
// overload with given parameters
dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_);
ds_max_ = param.getDefault("ds_max", ds_max_);
drs_max_rel_ = param.getDefault("drs_max_rel", drs_max_rel_);
dr_max_rel_ = param.getDefault("dr_max_rel", dr_max_rel_);
relax_max_ = param.getDefault("relax_max", relax_max_);
max_iter_ = param.getDefault("max_iter", max_iter_);
@ -1403,29 +1403,22 @@ namespace {
}
// Update rs and rv
const double drsmaxrel = drsMaxRel();
const double drvmax = 1e9;//% same as in Mrst
const double drmaxrel = drMaxRel();
V rs;
if (has_disgas_) {
const V rs_old = Eigen::Map<const V>(&state.gasoilratio()[0], nc);
const V drs = isRs * dxvar;
const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drsmaxrel);
const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drmaxrel);
rs = rs_old - drs_limited;
}
V rv;
if (has_vapoil_) {
const V rv_old = Eigen::Map<const V>(&state.rv()[0], nc);
const V drv = isRv * dxvar;
const V drv_limited = sign(drv) * drv.abs().min(drvmax);
const V drv_limited = sign(drv) * drv.abs().min(rv_old.abs()*drmaxrel);
rv = rv_old - drv_limited;
}
// Update the state
if (has_disgas_)
std::copy(&rs[0], &rs[0] + nc, state.gasoilratio().begin());
if (has_vapoil_)
std::copy(&rv[0], &rv[0] + nc, state.rv().begin());
// Sg is used as primal variable for water only cells.
const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
@ -1441,14 +1434,10 @@ namespace {
// The obvious case
auto hasGas = (sg > 0 && isRs == 0);
// keep oil saturated if previous sg is sufficient large:
const int pos = pu.phase_pos[ Gas ];
auto hadGas = (sg <= 0 && s_old.col(pos) > epsilon);
// Set oil saturated if previous rs is sufficiently large
const V rs_old = Eigen::Map<const V>(&state.gasoilratio()[0], nc);
auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
auto useSg = watOnly || hasGas || hadGas || gasVaporized;
auto useSg = watOnly || hasGas || gasVaporized;
for (int c = 0; c < nc; ++c) {
if (useSg[c]) { rs[c] = rsSat[c];}
else { primalVariable_[c] = PrimalVariables::RS; }
@ -1464,13 +1453,10 @@ namespace {
// The obvious case
auto hasOil = (so > 0 && isRv == 0);
// keep oil saturated if previous so is sufficient large:
const int pos = pu.phase_pos[ Oil ];
auto hadOil = (so <= 0 && s_old.col(pos) > epsilon );
// Set oil saturated if previous rv is sufficiently large
const V rv_old = Eigen::Map<const V>(&state.rv()[0], nc);
auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv == 1) && (rv_old > rvSat0 * (1-epsilon)) );
auto useSg = watOnly || hasOil || hadOil || oilCondensed;
auto useSg = watOnly || hasOil || oilCondensed;
for (int c = 0; c < nc; ++c) {
if (useSg[c]) { rv[c] = rvSat[c]; }
else {primalVariable_[c] = PrimalVariables::RV; }
@ -1479,6 +1465,14 @@ namespace {
}
// Update the state
if (has_disgas_) {
std::copy(&rs[0], &rs[0] + nc, state.gasoilratio().begin());
}
if (has_vapoil_) {
std::copy(&rv[0], &rv[0] + nc, state.rv().begin());
}
// Qs update.
// Since we need to update the wellrates, that are ordered by wells,