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

Refactor updatState() to use numWellVars().

Browse Source 
This allows us to use the base version of updateState().
This commit is contained in:
Atgeirr Flø Rasmussen 2015-11-20 11:45:35 +01:00 committed by Kai Bao
parent df30546841
commit 61cccfebf8
4 changed files with 29 additions and 233 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
opm/autodiff/BlackoilModelBase.hpp
View File

@ -345,6 +345,8 @@ namespace Opm {
// return wells object
const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; }
int numWellVars() const;
void
makeConstantState(SolutionState& state) const;

16
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -442,6 +442,19 @@ namespace detail {
template <class Grid, class Implementation>
int
BlackoilModelBase<Grid, Implementation>::numWellVars() const
{
// For each well, we have a bhp variable, and one flux per phase.
const int nw = localWellsActive() ? wells().number_of_wells : 0;
return (numPhases() + 1) * nw;
}
template <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::makeConstantState(SolutionState& state) const
@ -1958,7 +1971,6 @@ namespace detail {
using namespace Opm::AutoDiffGrid;
const int np = fluid_.numPhases();
const int nc = numCells(grid_);
const int nw = localWellsActive() ? wells().number_of_wells : 0;
const V null;
assert(null.size() == 0);
const V zero = V::Zero(nc);
@ -1973,7 +1985,7 @@ namespace detail {
varstart += dxvar.size();
// Extract well parts np phase rates + bhp
const V dwells = subset(dx, Span((np+1)*nw, 1, varstart));
const V dwells = subset(dx, Span(asImpl().numWellVars(), 1, varstart));
varstart += dwells.size();
assert(varstart == dx.size());

10
opm/autodiff/BlackoilMultiSegmentModel.hpp
View File

@ -106,14 +106,6 @@ namespace Opm {
WellState& well_state,
const bool initial_assembly);
/// Apply an update to the primary variables, chopped if appropriate.
/// \param[in] dx updates to apply to primary variables
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
void updateState(const V& dx,
ReservoirState& reservoir_state,
WellState& well_state);
using Base::numPhases;
using Base::numMaterials;
using Base::materialName;
@ -300,6 +292,8 @@ namespace Opm {
const WellState& xw,
const V& aliveWells);
int numWellVars() const;
void
makeConstantState(SolutionState& state) const;

234
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
View File

@ -274,6 +274,17 @@ namespace Opm {
template <class Grid>
int
BlackoilMultiSegmentModel<Grid>::numWellVars() const
{
// For each segment, we have a pressure variable, and one flux per phase.
const int nseg = wops_ms_.p2s.rows();
return (numPhases() + 1) * nseg;
}
template <class Grid>
@ -1484,229 +1495,6 @@ namespace Opm {
template <class Grid>
void BlackoilMultiSegmentModel<Grid>::updateState(const V& dx,
ReservoirState& reservoir_state,
WellState& well_state)
{
using namespace Opm::AutoDiffGrid;
const int np = fluid_.numPhases();
const int nc = numCells(grid_);
const V null;
assert(null.size() == 0);
const V zero = V::Zero(nc);
// Extract parts of dx corresponding to each part.
const V dp = subset(dx, Span(nc));
int varstart = nc;
const V dsw = active_[Water] ? subset(dx, Span(nc, 1, varstart)) : null;
varstart += dsw.size();
const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null;
varstart += dxvar.size();
// Extract well parts np phase rates + bhp
const int nseg_total = well_state.numSegments();
const V dwells = subset(dx, Span((np+1)*nseg_total, 1, varstart));
varstart += dwells.size();
assert(varstart == dx.size());
// Pressure update.
const double dpmaxrel = dpMaxRel();
const V p_old = Eigen::Map<const V>(&reservoir_state.pressure()[0], nc, 1);
const V absdpmax = dpmaxrel*p_old.abs();
const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
const V p = (p_old - dp_limited).max(zero);
std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
// Saturation updates.
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
const double dsmax = dsMax();
V so;
V sw;
V sg;
{
V maxVal = zero;
V dso = zero;
if (active_[Water]){
maxVal = dsw.abs().max(maxVal);
dso = dso - dsw;
}
V dsg;
if (active_[Gas]){
dsg = isSg_ * dxvar - isRv_ * dsw;
maxVal = dsg.abs().max(maxVal);
dso = dso - dsg;
}
maxVal = dso.abs().max(maxVal);
V step = dsmax/maxVal;
step = step.min(1.);
if (active_[Water]) {
const int pos = pu.phase_pos[ Water ];
const V sw_old = s_old.col(pos);
sw = sw_old - step * dsw;
}
if (active_[Gas]) {
const int pos = pu.phase_pos[ Gas ];
const V sg_old = s_old.col(pos);
sg = sg_old - step * dsg;
}
const int pos = pu.phase_pos[ Oil ];
const V so_old = s_old.col(pos);
so = so_old - step * dso;
}
// Appleyard chop process.
auto ixg = sg < 0;
for (int c = 0; c < nc; ++c) {
if (ixg[c]) {
sw[c] = sw[c] / (1-sg[c]);
so[c] = so[c] / (1-sg[c]);
sg[c] = 0;
}
}
auto ixo = so < 0;
for (int c = 0; c < nc; ++c) {
if (ixo[c]) {
sw[c] = sw[c] / (1-so[c]);
sg[c] = sg[c] / (1-so[c]);
so[c] = 0;
}
}
auto ixw = sw < 0;
for (int c = 0; c < nc; ++c) {
if (ixw[c]) {
so[c] = so[c] / (1-sw[c]);
sg[c] = sg[c] / (1-sw[c]);
sw[c] = 0;
}
}
//const V sumSat = sw + so + sg;
//sw = sw / sumSat;
//so = so / sumSat;
//sg = sg / sumSat;
// Update the reservoir_state
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
}
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
}
if (active_[ Oil ]) {
const int pos = pu.phase_pos[ Oil ];
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pos] = so[c];
}
}
// Update rs and rv
const double drmaxrel = drMaxRel();
V rs;
if (has_disgas_) {
const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
const V drs = isRs_ * dxvar;
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>(&reservoir_state.rv()[0], nc);
const V drv = isRv_ * dxvar;
const V drv_limited = sign(drv) * drv.abs().min(rv_old.abs()*drmaxrel);
rv = rv_old - drv_limited;
}
// Sg is used as primal variable for water only cells.
const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
auto watOnly = sw > (1 - epsilon);
// phase translation sg <-> rs
std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg);
if (has_disgas_) {
const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
const V rsSat = fluidRsSat(p, so, cells_);
// The obvious case
auto hasGas = (sg > 0 && isRs_ == 0);
// Set oil saturated if previous rs is sufficiently large
const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs_ == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
auto useSg = watOnly || hasGas || gasVaporized;
for (int c = 0; c < nc; ++c) {
if (useSg[c]) {
rs[c] = rsSat[c];
} else {
primalVariable_[c] = PrimalVariables::RS;
}
}
}
// phase transitions so <-> rv
if (has_vapoil_) {
// The gas pressure is needed for the rvSat calculations
const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas));
const V gaspress = computeGasPressure(p, sw, so, sg);
const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_);
const V rvSat = fluidRvSat(gaspress, so, cells_);
// The obvious case
auto hasOil = (so > 0 && isRv_ == 0);
// Set oil saturated if previous rv is sufficiently large
const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv_ == 1) && (rv_old > rvSat0 * (1-epsilon)) );
auto useSg = watOnly || hasOil || oilCondensed;
for (int c = 0; c < nc; ++c) {
if (useSg[c]) {
rv[c] = rvSat[c];
} else {
primalVariable_[c] = PrimalVariables::RV;
}
}
}
// Update the reservoir_state
if (has_disgas_) {
std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
}
if (has_vapoil_) {
std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
}
updateWellState(dwells,well_state);
// Update phase conditions used for property calculations.
updatePhaseCondFromPrimalVariable();
}
template <class Grid>
void
BlackoilMultiSegmentModel<Grid>::updateWellState(const V& dwells,