OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-15 11:23:24 -06:00
Code Issues Packages Projects Releases Wiki Activity

Remove functions that are identical in BlackoilModelBase.

Browse Source 
Also refactor some functions that are different to call the
base version and then do additional processing. However this
process has not been carried out on all methods at this point.
This commit is contained in:
Atgeirr Flø Rasmussen 2015-05-26 01:46:34 +02:00
parent 26484e91a5
commit f2e5177594
2 changed files with 19 additions and 958 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

132
opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
View File

@ -49,6 +49,9 @@ namespace Opm {
typedef BlackoilModelBase<Grid, BlackoilPolymerModel<Grid> > Base;
typedef typename Base::ReservoirState ReservoirState;
typedef typename Base::WellState WellState;
// The next line requires C++11 support available in g++ 4.7.
// friend Base;
friend BlackoilModelBase<Grid, BlackoilPolymerModel<Grid> >;
/// Construct the model. It will retain references to the
/// arguments of this functions, and they are expected to
@ -93,26 +96,12 @@ namespace Opm {
ReservoirState& reservoir_state,
WellState& well_state);
/// Assemble the residual and Jacobian of the nonlinear system.
/// \param[in] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
/// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep
void assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly);
// void assemble(const PolymerBlackoilState& reservoir_state,
// WellStateFullyImplicitBlackoilPolymer& well_state,
// const bool initial_assembly);
/// \brief Compute the residual norms of the mass balance for each phase,
/// the well flux, and the well equation.
/// \return a vector that contains for each phase the norm of the mass balance
/// and afterwards the norm of the residual of the well flux and the well equation.
std::vector<double> computeResidualNorms() const;
/// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual.
V solveJacobianSystem() const;
/// 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
@ -171,10 +160,18 @@ namespace Opm {
// Need to declare Base members we want to use here.
using Base::wellsActive;
using Base::wells;
SolutionState
constantState(const ReservoirState& x,
const WellState& xw) const;
using Base::computePressures;
using Base::computeGasPressure;
using Base::applyThresholdPressures;
using Base::fluidViscosity;
using Base::fluidReciprocFVF;
using Base::fluidDensity;
using Base::fluidRsSat;
using Base::fluidRvSat;
using Base::poroMult;
using Base::transMult;
using Base::updatePrimalVariableFromState;
using Base::updatePhaseCondFromPrimalVariable;
void
makeConstantState(SolutionState& state) const;
@ -187,39 +184,14 @@ namespace Opm {
computeAccum(const SolutionState& state,
const int aix );
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw);
void
addWellControlEq(const SolutionState& state,
const WellState& xw,
const V& aliveWells);
assembleMassBalanceEq(const SolutionState& state);
void
addWellEq(const SolutionState& state,
WellState& xw,
V& aliveWells);
void updateWellControls(WellState& xw) const;
std::vector<ADB>
computePressures(const SolutionState& state) const;
std::vector<ADB>
computePressures(const ADB& po,
const ADB& sw,
const ADB& so,
const ADB& sg) const;
V
computeGasPressure(const V& po,
const V& sw,
const V& so,
const V& sg) const;
std::vector<ADB>
computeRelPerm(const SolutionState& state) const;
void
computeMassFlux(const int actph ,
const V& transi,
@ -233,81 +205,9 @@ namespace Opm {
ADB
computeMc(const SolutionState& state) const;
void applyThresholdPressures(ADB& dp);
ADB
fluidViscosity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;
ADB
fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;
ADB
fluidDensity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;
V
fluidRsSat(const V& p,
const V& so,
const std::vector<int>& cells) const;
ADB
fluidRsSat(const ADB& p,
const ADB& so,
const std::vector<int>& cells) const;
V
fluidRvSat(const V& p,
const V& so,
const std::vector<int>& cells) const;
ADB
fluidRvSat(const ADB& p,
const ADB& so,
const std::vector<int>& cells) const;
ADB
poroMult(const ADB& p) const;
ADB
transMult(const ADB& p) const;
void
classifyCondition(const SolutionState& state,
std::vector<PhasePresence>& cond ) const;
const std::vector<PhasePresence>
phaseCondition() const {return this->phaseCondition_;}
void
classifyCondition(const ReservoirState& state);
/// update the primal variable for Sg, Rv or Rs. The Gas phase must
/// be active to call this method.
void
updatePrimalVariableFromState(const ReservoirState& state);
/// Update the phaseCondition_ member based on the primalVariable_ member.
void
updatePhaseCondFromPrimalVariable();
/// \brief Compute the reduction within the convergence check.
/// \param[in] B A matrix with MaxNumPhases columns and the same number rows
/// as the number of cells of the grid. B.col(i) contains the values

845
opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
View File

@ -134,44 +134,12 @@ namespace Opm {
template <class Grid>
typename BlackoilPolymerModel<Grid>::SolutionState
BlackoilPolymerModel<Grid>::constantState(const ReservoirState& x,
const WellState& xw) const
{
auto state = variableState(x, xw);
makeConstantState(state);
return state;
}
template <class Grid>
void
BlackoilPolymerModel<Grid>::makeConstantState(SolutionState& state) const
{
// HACK: throw away the derivatives. this may not be the most
// performant way to do things, but it will make the state
// automatically consistent with variableState() (and doing
// things automatically is all the rage in this module ;)
state.pressure = ADB::constant(state.pressure.value());
state.temperature = ADB::constant(state.temperature.value());
state.rs = ADB::constant(state.rs.value());
state.rv = ADB::constant(state.rv.value());
Base::makeConstantState(state);
state.concentration = ADB::constant(state.concentration.value());
const int num_phases = state.saturation.size();
for (int phaseIdx = 0; phaseIdx < num_phases; ++ phaseIdx) {
state.saturation[phaseIdx] = ADB::constant(state.saturation[phaseIdx].value());
}
state.qs = ADB::constant(state.qs.value());
state.bhp = ADB::constant(state.bhp.value());
assert(state.canonical_phase_pressures.size() == static_cast<std::size_t>(Opm::BlackoilPhases::MaxNumPhases));
for (int canphase = 0; canphase < Opm::BlackoilPhases::MaxNumPhases; ++canphase) {
ADB& pp = state.canonical_phase_pressures[canphase];
pp = ADB::constant(pp.value());
}
}
@ -421,189 +389,17 @@ namespace Opm {
template <class Grid>
void BlackoilPolymerModel<Grid>::computeWellConnectionPressures(const SolutionState& state,
const WellState& xw)
{
if( ! wellsActive() ) return ;
using namespace Opm::AutoDiffGrid;
// 1. Compute properties required by computeConnectionPressureDelta().
// Note that some of the complexity of this part is due to the function
// taking std::vector<double> arguments, and not Eigen objects.
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
// Compute the average pressure in each well block
const V perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
V avg_press = perf_press*0;
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
}
// Use cell values for the temperature as the wells don't knows its temperature yet.
const ADB perf_temp = subset(state.temperature, well_cells);
// Compute b, rsmax, rvmax values for perforations.
// Evaluate the properties using average well block pressures
// and cell values for rs, rv, phase condition and temperature.
const ADB avg_press_ad = ADB::constant(avg_press);
std::vector<PhasePresence> perf_cond(nperf);
const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf; ++perf) {
perf_cond[perf] = pc[well_cells[perf]];
}
const PhaseUsage& pu = fluid_.phaseUsage();
DataBlock b(nperf, pu.num_phases);
std::vector<double> rsmax_perf(nperf, 0.0);
std::vector<double> rvmax_perf(nperf, 0.0);
if (pu.phase_used[BlackoilPhases::Aqua]) {
const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
}
assert(active_[Oil]);
const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const ADB perf_rs = subset(state.rs, well_cells);
const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = subset(state.rv, well_cells);
const V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
}
// b is row major, so can just copy data.
std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases);
// Extract well connection depths.
const V depth = cellCentroidsZToEigen(grid_);
const V pdepth = subset(depth, well_cells);
std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
// Surface density.
std::vector<double> surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases);
// Gravity
double grav = 0.0;
const double* g = geo_.gravity();
const int dim = dimensions(grid_);
if (g) {
// Guard against gravity in anything but last dimension.
for (int dd = 0; dd < dim - 1; ++dd) {
assert(g[dd] == 0.0);
}
grav = g[dim - 1];
}
// 2. Compute pressure deltas, and store the results.
std::vector<double> cdp = WellDensitySegmented
::computeConnectionPressureDelta(wells(), xw, fluid_.phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, perf_depth,
surf_dens, grav);
well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf);
}
template <class Grid>
void
BlackoilPolymerModel<Grid>::
assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly)
assembleMassBalanceEq(const SolutionState& state)
{
using namespace Opm::AutoDiffGrid;
// Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells
updateWellControls(well_state);
// Create the primary variables.
SolutionState state = variableState(reservoir_state, well_state);
if (initial_assembly) {
// Create the (constant, derivativeless) initial state.
SolutionState state0 = state;
makeConstantState(state0);
// Compute initial accumulation contributions
// and well connection pressures.
computeAccum(state0, 0);
computeWellConnectionPressures(state0, well_state);
}
// OPM_AD_DISKVAL(state.pressure);
// OPM_AD_DISKVAL(state.saturation[0]);
// OPM_AD_DISKVAL(state.saturation[1]);
// OPM_AD_DISKVAL(state.saturation[2]);
// OPM_AD_DISKVAL(state.rs);
// OPM_AD_DISKVAL(state.rv);
// OPM_AD_DISKVAL(state.qs);
// OPM_AD_DISKVAL(state.bhp);
// -------- Mass balance equations --------
// Compute b_p and the accumulation term b_p*s_p for each phase,
// except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
// These quantities are stored in rq_[phase].accum[1].
// The corresponding accumulation terms from the start of
// the timestep (b^0_p*s^0_p etc.) were already computed
// on the initial call to assemble() and stored in rq_[phase].accum[0].
computeAccum(state, 1);
// Set up the common parts of the mass balance equations
// for each active phase.
const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
const std::vector<ADB> kr = computeRelPerm(state);
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
residual_.material_balance_eq[ phaseIdx ] =
pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
+ ops_.div*rq_[phaseIdx].mflux;
}
// -------- Extra (optional) rs and rv contributions to the mass balance equations --------
// Add the extra (flux) terms to the mass balance equations
// From gas dissolved in the oil phase (rs) and oil vaporized in the gas phase (rv)
// The extra terms in the accumulation part of the equation are already handled.
if (active_[ Oil ] && active_[ Gas ]) {
const int po = fluid_.phaseUsage().phase_pos[ Oil ];
const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
const UpwindSelector<double> upwindOil(grid_, ops_,
rq_[po].head.value());
const ADB rs_face = upwindOil.select(state.rs);
const UpwindSelector<double> upwindGas(grid_, ops_,
rq_[pg].head.value());
const ADB rv_face = upwindGas.select(state.rv);
residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * rq_[po].mflux);
residual_.material_balance_eq[ po ] += ops_.div * (rv_face * rq_[pg].mflux);
// OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]);
}
Base::assembleMassBalanceEq(state);
// Add polymer equation.
if (has_polymer_) {
residual_.material_balance_eq[poly_pos_] = pvdt_ * (rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0])
+ ops_.div*rq_[poly_pos_].mflux;
}
// Add contribution from wells and set up the well equations.
V aliveWells;
addWellEq(state, well_state, aliveWells);
addWellControlEq(state, well_state, aliveWells);
}
@ -787,161 +583,6 @@ namespace Opm {
template <class Grid>
void BlackoilPolymerModel<Grid>::updateWellControls(WellState& xw) const
{
if( ! wellsActive() ) return ;
std::string modestring[3] = { "BHP", "RESERVOIR_RATE", "SURFACE_RATE" };
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
int current = xw.currentControls()[w];
// Loop over all controls except the current one, and also
// skip any RESERVOIR_RATE controls, since we cannot
// handle those.
const int nwc = well_controls_get_num(wc);
int ctrl_index = 0;
for (; ctrl_index < nwc; ++ctrl_index) {
if (ctrl_index == current) {
// This is the currently used control, so it is
// used as an equation. So this is not used as an
// inequality constraint, and therefore skipped.
continue;
}
if (detail::constraintBroken(xw.bhp(), xw.wellRates(), w, np, wells().type[w], wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
}
}
if (ctrl_index != nwc) {
// Constraint number ctrl_index was broken, switch to it.
if (terminal_output_)
{
std::cout << "Switching control mode for well " << wells().name[w]
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
}
xw.currentControls()[w] = ctrl_index;
current = xw.currentControls()[w];
}
// Updating well state and primary variables.
// Target values are used as initial conditions for BHP and SURFACE_RATE
const double target = well_controls_iget_target(wc, current);
const double* distr = well_controls_iget_distr(wc, current);
switch (well_controls_iget_type(wc, current)) {
case BHP:
xw.bhp()[w] = target;
break;
case RESERVOIR_RATE:
// No direct change to any observable quantity at
// surface condition. In this case, use existing
// flow rates as initial conditions as reservoir
// rate acts only in aggregate.
break;
case SURFACE_RATE:
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0) {
xw.wellRates()[np*w + phase] = target * distr[phase];
}
}
break;
}
}
}
template <class Grid>
void BlackoilPolymerModel<Grid>::addWellControlEq(const SolutionState& state,
const WellState& xw,
const V& aliveWells)
{
if( ! wellsActive() ) return;
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
V bhp_targets = V::Zero(nw);
V rate_targets = V::Zero(nw);
M rate_distr(nw, np*nw);
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
const int current = xw.currentControls()[w];
switch (well_controls_iget_type(wc, current)) {
case BHP:
{
bhp_targets (w) = well_controls_iget_target(wc, current);
rate_targets(w) = -1e100;
}
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
{
// RESERVOIR and SURFACE rates look the same, from a
// high-level point of view, in the system of
// simultaneous linear equations.
const double* const distr =
well_controls_iget_distr(wc, current);
for (int p = 0; p < np; ++p) {
rate_distr.insert(w, p*nw + w) = distr[p];
}
bhp_targets (w) = -1.0e100;
rate_targets(w) = well_controls_iget_target(wc, current);
}
break;
}
}
const ADB bhp_residual = state.bhp - bhp_targets;
const ADB rate_residual = rate_distr * state.qs - rate_targets;
// Choose bhp residual for positive bhp targets.
Selector<double> bhp_selector(bhp_targets);
residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
// For wells that are dead (not flowing), and therefore not communicating
// with the reservoir, we set the equation to be equal to the well's total
// flow. This will be a solution only if the target rate is also zero.
M rate_summer(nw, np*nw);
for (int w = 0; w < nw; ++w) {
for (int phase = 0; phase < np; ++phase) {
rate_summer.insert(w, phase*nw + w) = 1.0;
}
}
Selector<double> alive_selector(aliveWells, Selector<double>::NotEqualZero);
residual_.well_eq = alive_selector.select(residual_.well_eq, rate_summer * state.qs);
// OPM_AD_DUMP(residual_.well_eq);
}
template <class Grid>
V BlackoilPolymerModel<Grid>::solveJacobianSystem() const
{
return linsolver_.computeNewtonIncrement(residual_);
}
template <class Grid>
void BlackoilPolymerModel<Grid>::updateState(const V& dx,
@ -1215,117 +856,6 @@ namespace Opm {
template <class Grid>
std::vector<ADB>
BlackoilPolymerModel<Grid>::computeRelPerm(const SolutionState& state) const
{
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const ADB zero = ADB::constant(V::Zero(nc));
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB& sw = (active_[ Water ]
? state.saturation[ pu.phase_pos[ Water ] ]
: zero);
const ADB& so = (active_[ Oil ]
? state.saturation[ pu.phase_pos[ Oil ] ]
: zero);
const ADB& sg = (active_[ Gas ]
? state.saturation[ pu.phase_pos[ Gas ] ]
: zero);
return fluid_.relperm(sw, so, sg, cells_);
}
template <class Grid>
std::vector<ADB>
BlackoilPolymerModel<Grid>::computePressures(const SolutionState& state) const
{
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const ADB null = ADB::constant(V::Zero(nc));
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB& sw = (active_[ Water ]
? state.saturation[ pu.phase_pos[ Water ] ]
: null);
const ADB& so = (active_[ Oil ]
? state.saturation[ pu.phase_pos[ Oil ] ]
: null);
const ADB& sg = (active_[ Gas ]
? state.saturation[ pu.phase_pos[ Gas ] ]
: null);
return computePressures(state.pressure, sw, so, sg);
}
template <class Grid>
std::vector<ADB>
BlackoilPolymerModel<Grid>::
computePressures(const ADB& po,
const ADB& sw,
const ADB& so,
const ADB& sg) const
{
// convert the pressure offsets to the capillary pressures
std::vector<ADB> pressure = fluid_.capPress(sw, so, sg, cells_);
for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
// The reference pressure is always the liquid phase (oil) pressure.
if (phaseIdx == BlackoilPhases::Liquid)
continue;
pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
}
// Since pcow = po - pw, but pcog = pg - po,
// we have
// pw = po - pcow
// pg = po + pcgo
// This is an unfortunate inconsistency, but a convention we must handle.
for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
if (phaseIdx == BlackoilPhases::Aqua) {
pressure[phaseIdx] = po - pressure[phaseIdx];
} else {
pressure[phaseIdx] += po;
}
}
return pressure;
}
template <class Grid>
V
BlackoilPolymerModel<Grid>::computeGasPressure(const V& po,
const V& sw,
const V& so,
const V& sg) const
{
assert (active_[Gas]);
std::vector<ADB> cp = fluid_.capPress(ADB::constant(sw),
ADB::constant(so),
ADB::constant(sg),
cells_);
return cp[Gas].value() + po;
}
template <class Grid>
void
BlackoilPolymerModel<Grid>::computeMassFlux(const int actph ,
@ -1389,36 +919,6 @@ namespace Opm {
template <class Grid>
void
BlackoilPolymerModel<Grid>::applyThresholdPressures(ADB& dp)
{
// We support reversible threshold pressures only.
// Method: if the potential difference is lower (in absolute
// value) than the threshold for any face, then the potential
// (and derivatives) is set to zero. If it is above the
// threshold, the threshold pressure is subtracted from the
// absolute potential (the potential is moved towards zero).
// Identify the set of faces where the potential is under the
// threshold, that shall have zero flow. Storing the bool
// Array as a V (a double Array) with 1 and 0 elements, a
// 1 where flow is allowed, a 0 where it is not.
const V high_potential = (dp.value().abs() >= threshold_pressures_by_interior_face_).template cast<double>();
// Create a sparse vector that nullifies the low potential elements.
const M keep_high_potential = spdiag(high_potential);
// Find the current sign for the threshold modification
const V sign_dp = sign(dp.value());
const V threshold_modification = sign_dp * threshold_pressures_by_interior_face_;
// Modify potential and nullify where appropriate.
dp = keep_high_potential * (dp - threshold_modification);
}
template <class Grid>
@ -1666,133 +1166,9 @@ namespace Opm {
}
template <class Grid>
ADB
BlackoilPolymerModel<Grid>::fluidViscosity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
switch (phase) {
case Water:
return fluid_.muWat(p, temp, cells);
case Oil: {
return fluid_.muOil(p, temp, rs, cond, cells);
}
case Gas:
return fluid_.muGas(p, temp, rv, cond, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
}
template <class Grid>
ADB
BlackoilPolymerModel<Grid>::fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
switch (phase) {
case Water:
return fluid_.bWat(p, temp, cells);
case Oil: {
return fluid_.bOil(p, temp, rs, cond, cells);
}
case Gas:
return fluid_.bGas(p, temp, rv, cond, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
}
template <class Grid>
ADB
BlackoilPolymerModel<Grid>::fluidDensity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
ADB b = fluidReciprocFVF(phase, p, temp, rs, rv, cond, cells);
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
if (phase == Oil && active_[Gas]) {
// It is correct to index into rhos with canonical phase indices.
rho += V::Constant(p.size(), 1, rhos[Gas]) * rs * b;
}
if (phase == Gas && active_[Oil]) {
// It is correct to index into rhos with canonical phase indices.
rho += V::Constant(p.size(), 1, rhos[Oil]) * rv * b;
}
return rho;
}
template <class Grid>
V
BlackoilPolymerModel<Grid>::fluidRsSat(const V& p,
const V& satOil,
const std::vector<int>& cells) const
{
return fluid_.rsSat(ADB::constant(p), ADB::constant(satOil), cells).value();
}
template <class Grid>
ADB
BlackoilPolymerModel<Grid>::fluidRsSat(const ADB& p,
const ADB& satOil,
const std::vector<int>& cells) const
{
return fluid_.rsSat(p, satOil, cells);
}
template <class Grid>
V
BlackoilPolymerModel<Grid>::fluidRvSat(const V& p,
const V& satOil,
const std::vector<int>& cells) const
{
return fluid_.rvSat(ADB::constant(p), ADB::constant(satOil), cells).value();
}
template <class Grid>
ADB
BlackoilPolymerModel<Grid>::fluidRvSat(const ADB& p,
const ADB& satOil,
const std::vector<int>& cells) const
{
return fluid_.rvSat(p, satOil, cells);
}
template <class Grid>
ADB
BlackoilPolymerModel<Grid>::computeMc(const SolutionState& state) const
@ -1802,221 +1178,6 @@ namespace Opm {
template <class Grid>
ADB
BlackoilPolymerModel<Grid>::poroMult(const ADB& p) const
{
const int n = p.size();
if (rock_comp_props_ && rock_comp_props_->isActive()) {
V pm(n);
V dpm(n);
for (int i = 0; i < n; ++i) {
pm[i] = rock_comp_props_->poroMult(p.value()[i]);
dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]);
}
ADB::M dpm_diag = spdiag(dpm);
const int num_blocks = p.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dpm_diag, p.derivative()[block], jacs[block]);
}
return ADB::function(std::move(pm), std::move(jacs));
} else {
return ADB::constant(V::Constant(n, 1.0));
}
}
template <class Grid>
ADB
BlackoilPolymerModel<Grid>::transMult(const ADB& p) const
{
const int n = p.size();
if (rock_comp_props_ && rock_comp_props_->isActive()) {
V tm(n);
V dtm(n);
for (int i = 0; i < n; ++i) {
tm[i] = rock_comp_props_->transMult(p.value()[i]);
dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]);
}
ADB::M dtm_diag = spdiag(dtm);
const int num_blocks = p.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dtm_diag, p.derivative()[block], jacs[block]);
}
return ADB::function(std::move(tm), std::move(jacs));
} else {
return ADB::constant(V::Constant(n, 1.0));
}
}
/*
template <class Grid>
void
FullyImplicitBlackoilSolver<T>::
classifyCondition(const SolutionState& state,
std::vector<PhasePresence>& cond ) const
{
const PhaseUsage& pu = fluid_.phaseUsage();
if (active_[ Gas ]) {
// Oil/Gas or Water/Oil/Gas system
const int po = pu.phase_pos[ Oil ];
const int pg = pu.phase_pos[ Gas ];
const V& so = state.saturation[ po ].value();
const V& sg = state.saturation[ pg ].value();
cond.resize(sg.size());
for (V::Index c = 0, e = sg.size(); c != e; ++c) {
if (so[c] > 0) { cond[c].setFreeOil (); }
if (sg[c] > 0) { cond[c].setFreeGas (); }
if (active_[ Water ]) { cond[c].setFreeWater(); }
}
}
else {
// Water/Oil system
assert (active_[ Water ]);
const int po = pu.phase_pos[ Oil ];
const V& so = state.saturation[ po ].value();
cond.resize(so.size());
for (V::Index c = 0, e = so.size(); c != e; ++c) {
cond[c].setFreeWater();
if (so[c] > 0) { cond[c].setFreeOil(); }
}
}
} */
template <class Grid>
void
BlackoilPolymerModel<Grid>::classifyCondition(const ReservoirState& state)
{
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const int np = state.numPhases();
const PhaseUsage& pu = fluid_.phaseUsage();
const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
if (active_[ Gas ]) {
// Oil/Gas or Water/Oil/Gas system
const V so = s.col(pu.phase_pos[ Oil ]);
const V sg = s.col(pu.phase_pos[ Gas ]);
for (V::Index c = 0, e = sg.size(); c != e; ++c) {
if (so[c] > 0) { phaseCondition_[c].setFreeOil (); }
if (sg[c] > 0) { phaseCondition_[c].setFreeGas (); }
if (active_[ Water ]) { phaseCondition_[c].setFreeWater(); }
}
}
else {
// Water/Oil system
assert (active_[ Water ]);
const V so = s.col(pu.phase_pos[ Oil ]);
for (V::Index c = 0, e = so.size(); c != e; ++c) {
phaseCondition_[c].setFreeWater();
if (so[c] > 0) { phaseCondition_[c].setFreeOil(); }
}
}
}
template <class Grid>
void
BlackoilPolymerModel<Grid>::updatePrimalVariableFromState(const ReservoirState& state)
{
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const int np = state.numPhases();
const PhaseUsage& pu = fluid_.phaseUsage();
const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
// Water/Oil/Gas system
assert (active_[ Gas ]);
// reset the primary variables if RV and RS is not set Sg is used as primary variable.
primalVariable_.resize(nc);
std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg);
const V sg = s.col(pu.phase_pos[ Gas ]);
const V so = s.col(pu.phase_pos[ Oil ]);
const V sw = s.col(pu.phase_pos[ Water ]);
const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
auto watOnly = sw > (1 - epsilon);
auto hasOil = so > 0;
auto hasGas = sg > 0;
// For oil only cells Rs is used as primal variable. For cells almost full of water
// the default primal variable (Sg) is used.
if (has_disgas_) {
for (V::Index c = 0, e = sg.size(); c != e; ++c) {
if ( !watOnly[c] && hasOil[c] && !hasGas[c] ) {primalVariable_[c] = PrimalVariables::RS; }
}
}
// For gas only cells Rv is used as primal variable. For cells almost full of water
// the default primal variable (Sg) is used.
if (has_vapoil_) {
for (V::Index c = 0, e = so.size(); c != e; ++c) {
if ( !watOnly[c] && hasGas[c] && !hasOil[c] ) {primalVariable_[c] = PrimalVariables::RV; }
}
}
updatePhaseCondFromPrimalVariable();
}
/// Update the phaseCondition_ member based on the primalVariable_ member.
template <class Grid>
void
BlackoilPolymerModel<Grid>::updatePhaseCondFromPrimalVariable()
{
if (! active_[Gas]) {
OPM_THROW(std::logic_error, "updatePhaseCondFromPrimarVariable() logic requires active gas phase.");
}
const int nc = primalVariable_.size();
for (int c = 0; c < nc; ++c) {
phaseCondition_[c] = PhasePresence(); // No free phases.
phaseCondition_[c].setFreeWater(); // Not necessary for property calculation usage.
switch (primalVariable_[c]) {
case PrimalVariables::Sg:
phaseCondition_[c].setFreeOil();
phaseCondition_[c].setFreeGas();
break;
case PrimalVariables::RS:
phaseCondition_[c].setFreeOil();
break;
case PrimalVariables::RV:
phaseCondition_[c].setFreeGas();
break;
default:
OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << c << ": " << primalVariable_[c]);
}
}
}
} // namespace Opm
#endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED