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 #685 from GitPaean/well_refactoring_interface_refactoring

Browse Source 
Well interface refactoring
This commit is contained in:
Atgeirr Flø Rasmussen 2016-05-24 10:58:57 +02:00
commit ad604bca89
29 changed files with 1120 additions and 1235 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
opm/autodiff/BlackoilModel.hpp
View File

@ -65,13 +65,13 @@ namespace Opm {
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo,
const RockCompressibility* rock_comp_props,
const Wells* wells_arg,
const StandardWells& std_wells,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output)
: Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver,
: Base(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
eclState, has_disgas, has_vapoil, terminal_output)
{
}

42
opm/autodiff/BlackoilModelBase.hpp
View File

@ -142,7 +142,7 @@ namespace Opm {
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const Wells* wells,
const WellModel& well_model,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
@ -276,7 +276,6 @@ namespace Opm {
const BlackoilPropsAdInterface& fluid_;
const DerivedGeology& geo_;
const RockCompressibility* rock_comp_props_;
WellModel std_wells_;
VFPProperties vfp_properties_;
const NewtonIterationBlackoilInterface& linsolver_;
// For each canonical phase -> true if active
@ -294,6 +293,10 @@ namespace Opm {
std::vector<ReservoirResidualQuant> rq_;
std::vector<PhasePresence> phaseCondition_;
// Well Model
WellModel well_model_;
V isRs_;
V isRv_;
V isSg_;
@ -328,19 +331,17 @@ namespace Opm {
}
/// return the WellModel object
WellModel& stdWells() { return std_wells_; }
const WellModel& stdWells() const { return std_wells_; }
WellModel& wellModel() { return well_model_; }
const WellModel& wellModel() const { return well_model_; }
/// return the Well struct in the WellModel
const Wells& wells() const { return std_wells_.wells(); }
const Wells& wells() const { return well_model_.wells(); }
/// return true if wells are available in the reservoir
bool wellsActive() const { return std_wells_.wellsActive(); }
bool wellsActive() const { return well_model_.wellsActive(); }
/// return true if wells are available on this process
bool localWellsActive() const { return std_wells_.localWellsActive(); }
int numWellVars() const;
bool localWellsActive() const { return well_model_.localWellsActive(); }
void
makeConstantState(SolutionState& state) const;
@ -364,11 +365,6 @@ namespace Opm {
const std::vector<int>& indices,
std::vector<ADB>& vars) const;
void
variableStateExtractWellsVars(const std::vector<int>& indices,
std::vector<ADB>& vars,
SolutionState& state) const;
void
computeAccum(const SolutionState& state,
const int aix );
@ -376,18 +372,6 @@ namespace Opm {
void
assembleMassBalanceEq(const SolutionState& state);
// TODO: only kept for now due to flow_multisegment
// will be removed soon
void
extractWellPerfProperties(const SolutionState& state,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const;
// TODO: only kept for now due to flow_multisegment
// will be removed soon
void updateWellState(const V& dwells,
WellState& well_state);
bool
solveWellEq(const std::vector<ADB>& mob_perfcells,
@ -495,6 +479,12 @@ namespace Opm {
void
updatePhaseCondFromPrimalVariable(const ReservoirState& state);
// TODO: added since the interfaces of the function are different
// TODO: for StandardWells and MultisegmentWells
void
computeWellConnectionPressures(const SolutionState& state,
const WellState& well_state);
/// \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

205
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -142,6 +142,21 @@ namespace detail {
return act2can;
}
inline
double getGravity(const double* g, const int dim) {
double grav = 0.0;
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];
}
return grav;
}
} // namespace detail
@ -152,7 +167,7 @@ namespace detail {
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const Wells* wells_arg,
const WellModel& well_model,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
@ -162,7 +177,6 @@ namespace detail {
, fluid_ (fluid)
, geo_ (geo)
, rock_comp_props_(rock_comp_props)
, std_wells_ (wells_arg)
, vfp_properties_(
eclState->getTableManager().getVFPInjTables(),
eclState->getTableManager().getVFPProdTables())
@ -177,6 +191,7 @@ namespace detail {
, use_threshold_pressure_(false)
, rq_ (fluid.numPhases())
, phaseCondition_(AutoDiffGrid::numCells(grid))
, well_model_ (well_model)
, isRs_(V::Zero(AutoDiffGrid::numCells(grid)))
, isRv_(V::Zero(AutoDiffGrid::numCells(grid)))
, isSg_(V::Zero(AutoDiffGrid::numCells(grid)))
@ -201,6 +216,15 @@ namespace detail {
assert(numMaterials() == std::accumulate(active_.begin(), active_.end(), 0)); // Due to the material_name_ init above.
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_);
well_model_.init(&fluid_, &active_, &phaseCondition_, &vfp_properties_, gravity, depth);
// TODO: put this for now to avoid modify the following code.
// TODO: this code can be fragile.
const Wells* wells_arg = &(asImpl().well_model_.wells());
#if HAVE_MPI
if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
{
@ -213,7 +237,7 @@ namespace detail {
int local_number_of_wells = localWellsActive() ? wells().number_of_wells : 0;
int global_number_of_wells = info.communicator().sum(local_number_of_wells);
const bool wells_active = ( wells_arg && global_number_of_wells > 0 );
stdWells().setWellsActive(wells_active);
wellModel().setWellsActive(wells_active);
// Compute the global number of cells
std::vector<int> v( Opm::AutoDiffGrid::numCells(grid_), 1);
global_nc_ = 0;
@ -221,7 +245,7 @@ namespace detail {
}else
#endif
{
stdWells().setWellsActive( localWellsActive() );
wellModel().setWellsActive( localWellsActive() );
global_nc_ = Opm::AutoDiffGrid::numCells(grid_);
}
}
@ -436,19 +460,6 @@ namespace detail {
template <class Grid, class WellModel, class Implementation>
int
BlackoilModelBase<Grid, WellModel, Implementation>::numWellVars() const
{
// For each well, we have a bhp variable, and one flux per phase.
const int nw = stdWells().localWellsActive() ? wells().number_of_wells : 0;
return (numPhases() + 1) * nw;
}
template <class Grid, class WellModel, class Implementation>
void
BlackoilModelBase<Grid, WellModel, Implementation>::
@ -509,7 +520,7 @@ namespace detail {
// and bhp and Q for the wells
vars0.reserve(np + 1);
variableReservoirStateInitials(x, vars0);
asImpl().stdWells().variableWellStateInitials(xw, vars0);
asImpl().wellModel().variableWellStateInitials(xw, vars0);
return vars0;
}
@ -571,7 +582,7 @@ namespace detail {
if (active_[Gas]) {
indices[Xvar] = next++;
}
asImpl().stdWells().variableStateWellIndices(indices, next);
asImpl().wellModel().variableStateWellIndices(indices, next);
assert(next == fluid_.numPhases() + 2);
return indices;
}
@ -657,7 +668,7 @@ namespace detail {
}
}
// wells
asImpl().variableStateExtractWellsVars(indices, vars, state);
asImpl().wellModel().variableStateExtractWellsVars(indices, vars, state);
return state;
}
@ -665,24 +676,6 @@ namespace detail {
template <class Grid, class WellModel, class Implementation>
void
BlackoilModelBase<Grid, WellModel, Implementation>::
variableStateExtractWellsVars(const std::vector<int>& indices,
std::vector<ADB>& vars,
SolutionState& state) const
{
// Qs.
state.qs = std::move(vars[indices[Qs]]);
// Bhp.
state.bhp = std::move(vars[indices[Bhp]]);
}
template <class Grid, class WellModel, class Implementation>
void
BlackoilModelBase<Grid, WellModel, Implementation>::
@ -727,21 +720,6 @@ namespace detail {
}
}
namespace detail {
inline
double getGravity(const double* g, const int dim) {
double grav = 0.0;
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];
}
return grav;
}
}
@ -755,9 +733,6 @@ namespace detail {
{
using namespace Opm::AutoDiffGrid;
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const V depth = cellCentroidsZToEigen(grid_);
// If we have VFP tables, we need the well connection
// pressures for the "simple" hydrostatic correction
// between well depth and vfp table depth.
@ -765,16 +740,12 @@ namespace detail {
SolutionState state = asImpl().variableState(reservoir_state, well_state);
SolutionState state0 = state;
asImpl().makeConstantState(state0);
// asImpl().computeWellConnectionPressures(state0, well_state);
// Extract well connection depths.
asImpl().stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity);
asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
}
// Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells
// asImpl().updateWellControls(well_state);
// asImpl().stdWells().updateWellControls(well_state);
asImpl().stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state);
asImpl().wellModel().updateWellControls(terminal_output_, well_state);
// Create the primary variables.
SolutionState state = asImpl().variableState(reservoir_state, well_state);
@ -786,8 +757,7 @@ namespace detail {
// Compute initial accumulation contributions
// and well connection pressures.
asImpl().computeAccum(state0, 0);
// asImpl().computeWellConnectionPressures(state0, well_state);
asImpl().stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity);
asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
}
// OPM_AD_DISKVAL(state.pressure);
@ -810,25 +780,23 @@ namespace detail {
std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells;
asImpl().stdWells().extractWellPerfProperties(state, rq_, fluid_.numPhases(), fluid_, active_, mob_perfcells, b_perfcells);
asImpl().wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
}
V aliveWells;
std::vector<ADB> cq_s;
asImpl().stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s);
asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
asImpl().stdWells().addWellFluxEq(cq_s, state, residual_);
asImpl().wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
asImpl().wellModel().addWellFluxEq(cq_s, state, residual_);
asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
asImpl().stdWells().addWellControlEq(state, well_state, aliveWells, active_, vfp_properties_, gravity, residual_);
// asImpl().computeWellPotentials(state, mob_perfcells, b_perfcells, well_state);
{
asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
if (param_.compute_well_potentials_) {
SolutionState state0 = state;
asImpl().makeConstantState(state0);
asImpl().stdWells().computeWellPotentials(state0, mob_perfcells, b_perfcells,
fluid_.phaseUsage(), active_, vfp_properties_,
param_.compute_well_potentials_, gravity, well_state);
asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state);
}
}
@ -956,37 +924,7 @@ namespace detail {
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = asImpl().numPhases();
for (int phase = 0; phase < np; ++phase) {
residual_.material_balance_eq[phase] -= superset(cq_s[phase], stdWells().wellOps().well_cells, nc);
}
}
template <class Grid, class WellModel, class Implementation>
void
BlackoilModelBase<Grid, WellModel, Implementation>::
extractWellPerfProperties(const SolutionState&,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const
{
// If we have wells, extract the mobilities and b-factors for
// the well-perforated cells.
if (!asImpl().localWellsActive()) {
mob_perfcells.clear();
b_perfcells.clear();
return;
} else {
const int np = asImpl().numPhases();
const std::vector<int>& well_cells = stdWells().wellOps().well_cells;
mob_perfcells.resize(np, ADB::null());
b_perfcells.resize(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
mob_perfcells[phase] = subset(rq_[phase].mob, well_cells);
b_perfcells[phase] = subset(rq_[phase].b, well_cells);
}
residual_.material_balance_eq[phase] -= superset(cq_s[phase], wellModel().wellOps().well_cells, nc);
}
}
@ -1041,7 +979,7 @@ namespace detail {
V aliveWells;
const int np = wells().number_of_phases;
std::vector<ADB> cq_s(np, ADB::null());
std::vector<int> indices = asImpl().stdWells().variableWellStateIndices();
std::vector<int> indices = asImpl().wellModel().variableWellStateIndices();
SolutionState state0 = state;
WellState well_state0 = well_state;
asImpl().makeConstantState(state0);
@ -1058,25 +996,21 @@ namespace detail {
}
}
// gravity
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
int it = 0;
bool converged;
do {
// bhp and Q for the wells
std::vector<V> vars0;
vars0.reserve(2);
asImpl().stdWells().variableWellStateInitials(well_state, vars0);
asImpl().wellModel().variableWellStateInitials(well_state, vars0);
std::vector<ADB> vars = ADB::variables(vars0);
SolutionState wellSolutionState = state0;
asImpl().variableStateExtractWellsVars(indices, vars, wellSolutionState);
asImpl().stdWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
asImpl().stdWells().addWellFluxEq(cq_s, wellSolutionState, residual_);
asImpl().stdWells().addWellControlEq(wellSolutionState, well_state, aliveWells,
active_, vfp_properties_, gravity, residual_);
asImpl().wellModel().variableStateExtractWellsVars(indices, vars, wellSolutionState);
asImpl().wellModel().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
asImpl().wellModel().addWellFluxEq(cq_s, wellSolutionState, residual_);
asImpl().wellModel().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_);
converged = getWellConvergence(it);
if (converged) {
@ -1099,9 +1033,8 @@ namespace detail {
ADB::V total_residual_v = total_residual.value();
const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
assert(dx.size() == total_residual_v.size());
// asImpl().updateWellState(dx.array(), well_state);
asImpl().stdWells().updateWellState(dx.array(), gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state);
asImpl().stdWells(). updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state);
asImpl().wellModel().updateWellState(dx.array(), dpMaxRel(), well_state);
asImpl().wellModel().updateWellControls(terminal_output_, well_state);
}
} while (it < 15);
@ -1128,9 +1061,7 @@ namespace detail {
std::vector<ADB::M> old_derivs = state.qs.derivative();
state.qs = ADB::function(std::move(new_qs), std::move(old_derivs));
}
// asImpl().computeWellConnectionPressures(state, well_state);
const ADB::V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_);
asImpl().stdWells().computeWellConnectionPressures(state, well_state, fluid_, active_, phaseCondition(), depth, gravity);
asImpl().computeWellConnectionPressures(state, well_state);
}
if (!converged) {
@ -1281,7 +1212,7 @@ namespace detail {
varstart += dxvar.size();
// Extract well parts np phase rates + bhp
const V dwells = subset(dx, Span(asImpl().numWellVars(), 1, varstart));
const V dwells = subset(dx, Span(asImpl().wellModel().numWellVars(), 1, varstart));
varstart += dwells.size();
assert(varstart == dx.size());
@ -1488,10 +1419,8 @@ namespace detail {
std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
}
// TODO: gravity should be stored as a member
// const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
// asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state);
asImpl().updateWellState(dwells,well_state);
asImpl().wellModel().updateWellState(dwells, dpMaxRel(), well_state);
// Update phase conditions used for property calculations.
updatePhaseCondFromPrimalVariable(reservoir_state);
@ -2324,20 +2253,14 @@ namespace detail {
// TODO: only kept for now due to flow_multisegment
// will be removed soon
template <class Grid, class WellModel, class Implementation>
void
BlackoilModelBase<Grid, WellModel, Implementation>::
updateWellState(const V& dwells,
WellState& well_state)
{
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(),
active_, vfp_properties_, well_state);
}
template <class Grid, class WellModel, class Implementation>
void
BlackoilModelBase<Grid, WellModel, Implementation>::
computeWellConnectionPressures(const SolutionState& state,
const WellState& well_state)
{
asImpl().wellModel().computeWellConnectionPressures(state, well_state);
}

62
opm/autodiff/BlackoilMultiSegmentModel.hpp
View File

@ -51,11 +51,11 @@ namespace Opm {
/// \tparam Grid UnstructuredGrid or CpGrid.
/// \tparam Implementation Provides concrete state types.
template<class Grid>
class BlackoilMultiSegmentModel : public BlackoilModelBase<Grid, StandardWells, BlackoilMultiSegmentModel<Grid> >
class BlackoilMultiSegmentModel : public BlackoilModelBase<Grid, MultisegmentWells, BlackoilMultiSegmentModel<Grid> >
{
public:
typedef BlackoilModelBase<Grid, StandardWells, BlackoilMultiSegmentModel<Grid> > Base; // base class
typedef BlackoilModelBase<Grid, MultisegmentWells, BlackoilMultiSegmentModel<Grid> > Base; // base class
typedef typename Base::ReservoirState ReservoirState;
typedef typename Base::WellState WellState;
typedef BlackoilMultiSegmentSolutionState SolutionState;
@ -85,13 +85,12 @@ namespace Opm {
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const Wells* wells,
const MultisegmentWells& well_model,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output,
const std::vector<WellMultiSegmentConstPtr>& wells_multisegment);
const bool terminal_output);
/// Called once before each time step.
/// \param[in] dt time step size
@ -136,11 +135,12 @@ namespace Opm {
using Base::cells_;
using Base::param_;
using Base::linsolver_;
using Base::phaseCondition_;
using Base::vfp_properties_;
using Base::well_model_;
MultisegmentWells ms_wells_;
using Base::stdWells;
using Base::wells;
using Base::wellModel;
// using Base::wells;
using Base::wellsActive;
using Base::updatePrimalVariableFromState;
using Base::phaseCondition;
@ -159,35 +159,10 @@ namespace Opm {
using Base::asImpl;
using Base::variableReservoirStateInitials;
// TODO: fixing the confusing naming
const MultisegmentWells& msWells() const { return ms_wells_; }
MultisegmentWells& msWells() { return ms_wells_; }
const std::vector<WellMultiSegmentConstPtr>& wellsMultiSegment() const { return msWells().wells(); }
const std::vector<WellMultiSegmentConstPtr>& wellsMultiSegment() const { return well_model_.msWells(); }
const MultisegmentWells::MultisegmentWellOps& msWellOps() const { return msWells().wellOps(); }
// TODO: kept for now. to be removed soon.
void updateWellState(const V& dwells,
WellState& well_state);
std::vector<V>
variableStateInitials(const ReservoirState& x,
const WellState& xw) const;
void
variableWellStateInitials(const WellState& xw,
std::vector<V>& vars0) const;
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw);
/// added to fixing the flow_multisegment running
bool
baseSolveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
SolutionState& state,
WellState& well_state);
const MultisegmentWells::MultisegmentWellOps& msWellOps() const { return well_model_.wellOps(); }
bool
solveWellEq(const std::vector<ADB>& mob_perfcells,
@ -195,21 +170,14 @@ namespace Opm {
SolutionState& state,
WellState& well_state);
void
updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw) const;
int numWellVars() const;
void
makeConstantState(SolutionState& state) const;
// TODO: added since the interfaces of the function are different
// TODO: for StandardWells and MultisegmentWells
void
variableStateExtractWellsVars(const std::vector<int>& indices,
std::vector<ADB>& vars,
SolutionState& state) const;
computeWellConnectionPressures(const SolutionState& state,
const WellState& well_state);
};

537
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
View File

@ -62,16 +62,14 @@ namespace Opm {
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const Wells* wells_arg,
const MultisegmentWells& well_model,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output,
const std::vector<WellMultiSegmentConstPtr>& wells_multisegment)
: Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver,
const bool terminal_output)
: Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver,
eclState, has_disgas, has_vapoil, terminal_output)
, ms_wells_(wells_multisegment, fluid.numPhases())
{
}
@ -94,7 +92,7 @@ namespace Opm {
const int nw = wellsMultiSegment().size();
if ( !msWellOps().has_multisegment_wells ) {
msWells().segVDt() = V::Zero(nw);
wellModel().segVDt() = V::Zero(nw);
return;
}
@ -107,20 +105,7 @@ namespace Opm {
segment_volume.insert(segment_volume.end(), segment_volume_well.begin(), segment_volume_well.end());
}
assert(int(segment_volume.size()) == nseg_total);
msWells().segVDt() = Eigen::Map<V>(segment_volume.data(), nseg_total) / dt;
}
template <class Grid>
int
BlackoilMultiSegmentModel<Grid>::numWellVars() const
{
// For each segment, we have a pressure variable, and one flux per phase.
const int nseg = msWellOps().p2s.rows();
return (numPhases() + 1) * nseg;
wellModel().segVDt() = Eigen::Map<V>(segment_volume.data(), nseg_total) / dt;
}
@ -140,286 +125,6 @@ namespace Opm {
template <class Grid>
void
BlackoilMultiSegmentModel<Grid>::variableWellStateInitials(const WellState& xw, std::vector<V>& vars0) const
{
// Initial well rates
if ( wellsMultiSegment().size() > 0 )
{
// Need to reshuffle well segment rates, from phase running fastest
const int nseg = xw.numSegments();
const int np = xw.numPhases();
// The transpose() below switches the ordering of the segment rates
const DataBlock segrates = Eigen::Map<const DataBlock>(& xw.segPhaseRates()[0], nseg, np).transpose();
// segment phase rates in surface volume
const V segqs = Eigen::Map<const V>(segrates.data(), nseg * np);
vars0.push_back(segqs);
// for the pressure of the segments
const V segp = Eigen::Map<const V>(& xw.segPress()[0], xw.segPress().size());
vars0.push_back(segp);
}
else
{
// push null sates for segqs and segp
vars0.push_back(V());
vars0.push_back(V());
}
}
template <class Grid>
void
BlackoilMultiSegmentModel<Grid>::variableStateExtractWellsVars(const std::vector<int>& indices,
std::vector<ADB>& vars,
SolutionState& state) const
{
// TODO: using the original Qs for the segment rates for now, to be fixed eventually.
// TODO: using the original Bhp for the segment pressures for now, to be fixed eventually.
// segment phase rates in surface volume
state.segqs = std::move(vars[indices[Qs]]);
// segment pressures
state.segp = std::move(vars[indices[Bhp]]);
// The qs and bhp are no longer primary variables, but could
// still be used in computations. They are identical to the
// pressures and flows of the top segments.
const int np = numPhases();
const int ns = state.segp.size();
const int nw = msWells().topWellSegments().size();
state.qs = ADB::constant(ADB::V::Zero(np*nw));
for (int phase = 0; phase < np; ++phase) {
// Extract segment fluxes for this phase (ns consecutive elements).
ADB segqs_phase = subset(state.segqs, Span(ns, 1, ns*phase));
// Extract top segment fluxes (= well fluxes)
ADB wellqs_phase = subset(segqs_phase, msWells().topWellSegments());
// Expand to full size of qs (which contains all phases) and add.
state.qs += superset(wellqs_phase, Span(nw, 1, nw*phase), nw*np);
}
state.bhp = subset(state.segp, msWells().topWellSegments());
}
// TODO: This is just a preliminary version, remains to be improved later when we decide a better way
// TODO: to intergrate the usual wells and multi-segment wells.
template <class Grid>
void BlackoilMultiSegmentModel<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_total = xw.numPerforations();
const int nw = xw.numWells();
const std::vector<int>& well_cells = msWellOps().well_cells;
stdWells().wellPerforationDensities() = V::Zero(nperf_total);
const V perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf_total);
V avg_press = perf_press * 0.0;
// for the non-segmented/regular wells, calculated the average pressures.
// If it is the top perforation, then average with the bhp().
// If it is not the top perforation, then average with the perforation above it().
int start_segment = 0;
for (int w = 0; w < nw; ++w) {
const int nseg = wellsMultiSegment()[w]->numberOfSegments();
if (wellsMultiSegment()[w]->isMultiSegmented()) {
// maybe we should give some reasonable values to prevent the following calculations fail
start_segment += nseg;
continue;
}
std::string well_name(wellsMultiSegment()[w]->name());
typedef typename WellStateMultiSegment::SegmentedWellMapType::const_iterator const_iterator;
const_iterator it_well = xw.segmentedWellMap().find(well_name);
assert(it_well != xw.segmentedWellMap().end());
const int start_perforation = (*it_well).second.start_perforation;
const int end_perforation = start_perforation + (*it_well).second.number_of_perforations;
for (int perf = start_perforation; perf < end_perforation; ++perf) {
const double p_above = perf == start_perforation ? state.segp.value()[start_segment] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
start_segment += nseg;
}
assert(start_segment == xw.numSegments());
// 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_total);
const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf_total; ++perf) {
perf_cond[perf] = pc[well_cells[perf]];
}
const PhaseUsage& pu = fluid_.phaseUsage();
DataBlock b(nperf_total, pu.num_phases);
std::vector<double> rsmax_perf(nperf_total, 0.0);
std::vector<double> rvmax_perf(nperf_total, 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_total);
}
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_total);
}
// b is row major, so can just copy data.
std::vector<double> b_perf(b.data(), b.data() + nperf_total * pu.num_phases);
// Extract well connection depths.
const V depth = cellCentroidsZToEigen(grid_);
const V perfcelldepth = subset(depth, well_cells);
std::vector<double> perf_cell_depth(perfcelldepth.data(), perfcelldepth.data() + nperf_total);
// Surface density.
// The compute density segment wants the surface densities as
// an np * number of wells cells array
V rho = superset(fluid_.surfaceDensity(0 , well_cells), Span(nperf_total, pu.num_phases, 0), nperf_total * pu.num_phases);
for (int phase = 1; phase < pu.num_phases; ++phase) {
rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf_total, pu.num_phases, phase), nperf_total * pu.num_phases);
}
std::vector<double> surf_dens_perf(rho.data(), rho.data() + nperf_total * pu.num_phases);
// Gravity
double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
// 2. Compute densities
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_.phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// 3. Compute pressure deltas
std::vector<double> cdp =
WellDensitySegmented::computeConnectionPressureDelta(
wells(), perf_cell_depth, cd, grav);
// 4. Store the results
stdWells().wellPerforationDensities() = Eigen::Map<const V>(cd.data(), nperf_total); // This one is not useful for segmented wells at all
stdWells().wellPerforationPressureDiffs() = Eigen::Map<const V>(cdp.data(), nperf_total);
if ( !msWellOps().has_multisegment_wells ) {
msWells().wellPerforationCellDensities() = V::Zero(nperf_total);
msWells().wellPerforationCellPressureDiffs() = V::Zero(nperf_total);
return;
}
// compute the average of the fluid densites in the well blocks.
// the average is weighted according to the fluid relative permeabilities.
const std::vector<ADB> kr_adb = Base::computeRelPerm(state);
size_t temp_size = kr_adb.size();
std::vector<V> perf_kr;
for(size_t i = 0; i < temp_size; ++i) {
// const ADB kr_phase_adb = subset(kr_adb[i], well_cells);
const V kr_phase = (subset(kr_adb[i], well_cells)).value();
perf_kr.push_back(kr_phase);
}
// compute the averaged density for the well block
// TODO: for the non-segmented wells, they should be set to zero
// TODO: for the moment, they are still calculated, while not used later.
for (int i = 0; i < nperf_total; ++i) {
double sum_kr = 0.;
int np = perf_kr.size(); // make sure it is 3
for (int p = 0; p < np; ++p) {
sum_kr += perf_kr[p][i];
}
for (int p = 0; p < np; ++p) {
perf_kr[p][i] /= sum_kr;
}
}
V rho_avg_perf = V::Constant(nperf_total, 0.0);
// TODO: make sure the order of the density and the order of the kr are the same.
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
const int canonicalPhaseIdx = canph_[phaseIdx];
const ADB fluid_density = fluidDensity(canonicalPhaseIdx, rq_[phaseIdx].b, state.rs, state.rv);
const V rho_perf = subset(fluid_density, well_cells).value();
// TODO: phaseIdx or canonicalPhaseIdx ?
rho_avg_perf += rho_perf * perf_kr[phaseIdx];
}
msWells().wellPerforationCellDensities() = Eigen::Map<const V>(rho_avg_perf.data(), nperf_total);
// We should put this in a global class
std::vector<double> perf_depth_vec;
perf_depth_vec.reserve(nperf_total);
for (int w = 0; w < nw; ++w) {
WellMultiSegmentConstPtr well = wellsMultiSegment()[w];
const std::vector<double>& perf_depth_well = well->perfDepth();
perf_depth_vec.insert(perf_depth_vec.end(), perf_depth_well.begin(), perf_depth_well.end());
}
assert(int(perf_depth_vec.size()) == nperf_total);
const V perf_depth = Eigen::Map<V>(perf_depth_vec.data(), nperf_total);
const V perf_cell_depth_diffs = perf_depth - perfcelldepth;
msWells().wellPerforationCellPressureDiffs() = grav * msWells().wellPerforationCellDensities() * perf_cell_depth_diffs;
// Calculating the depth difference between segment nodes and perforations.
// TODO: should be put somewhere else for better clarity later
msWells().wellSegmentPerforationDepthDiffs() = V::Constant(nperf_total, -1e100);
int start_perforation = 0;
for (int w = 0; w < nw; ++w) {
WellMultiSegmentConstPtr well = wellsMultiSegment()[w];
const int nseg = well->numberOfSegments();
const int nperf = well->numberOfPerforations();
const std::vector<std::vector<int>>& segment_perforations = well->segmentPerforations();
for (int s = 0; s < nseg; ++s) {
const int nperf_seg = segment_perforations[s].size();
const double segment_depth = well->segmentDepth()[s];
for (int perf = 0; perf < nperf_seg; ++perf) {
const int perf_number = segment_perforations[s][perf] + start_perforation;
msWells().wellSegmentPerforationDepthDiffs()[perf_number] = segment_depth - perf_depth[perf_number];
}
}
start_perforation += nperf;
}
assert(start_perforation == nperf_total);
}
template <class Grid>
void
BlackoilMultiSegmentModel<Grid>::
@ -442,7 +147,7 @@ namespace Opm {
// Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells
msWells().updateWellControls(terminal_output_, well_state);
wellModel().updateWellControls(terminal_output_, well_state);
// Create the primary variables.
SolutionState state = asImpl().variableState(reservoir_state, well_state);
@ -454,12 +159,13 @@ namespace Opm {
// Compute initial accumulation contributions
// and well connection pressures.
asImpl().computeAccum(state0, 0);
msWells().computeSegmentFluidProperties(state0, phaseCondition(), active_, fluid_);
wellModel().computeSegmentFluidProperties(state0);
const int np = numPhases();
assert(np == int(msWells().segmentCompSurfVolumeInitial().size()));
assert(np == int(wellModel().segmentCompSurfVolumeInitial().size()));
for (int phase = 0; phase < np; ++phase) {
msWells().segmentCompSurfVolumeInitial()[phase] = msWells().segmentCompSurfVolumeCurrent()[phase].value();
wellModel().segmentCompSurfVolumeInitial()[phase] = wellModel().segmentCompSurfVolumeCurrent()[phase].value();
}
asImpl().computeWellConnectionPressures(state0, well_state);
}
@ -481,16 +187,14 @@ namespace Opm {
return;
}
// asImpl().computeSegmentFluidProperties(state);
msWells().computeSegmentFluidProperties(state, phaseCondition(), active_, fluid_);
wellModel().computeSegmentFluidProperties(state);
// asImpl().computeSegmentPressuresDelta(state);
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
msWells().computeSegmentPressuresDelta(gravity);
wellModel().computeSegmentPressuresDelta(gravity);
std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells;
asImpl().extractWellPerfProperties(state, mob_perfcells, b_perfcells);
wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
@ -500,65 +204,11 @@ namespace Opm {
// it is related to the segment location
V aliveWells;
std::vector<ADB> cq_s;
const int nw = wellsMultiSegment().size();
const int np = numPhases();
const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
const V perf_press_diffs = stdWells().wellPerforationPressureDiffs();
msWells().computeWellFlux(state, fluid_.phaseUsage(), active_,
perf_press_diffs, compi,
mob_perfcells, b_perfcells, aliveWells, cq_s);
asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
msWells().addWellFluxEq(cq_s, state, residual_);
wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
wellModel().addWellFluxEq(cq_s, state, residual_);
asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
// asImpl().addWellControlEq(state, well_state, aliveWells);
msWells().addWellControlEq(state, well_state, aliveWells, active_, residual_);
}
template <class Grid>
void BlackoilMultiSegmentModel<Grid>::updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw) const
{
// Update the perforation phase rates (used to calculate the pressure drop in the wellbore).
const int np = numPhases();
const int nw = wellsMultiSegment().size();
const int nperf_total = xw.perfPress().size();
V cq = superset(cq_s[0].value(), Span(nperf_total, np, 0), nperf_total * np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(cq_s[phase].value(), Span(nperf_total, np, phase), nperf_total * np);
}
xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf_total * np);
// Update the perforation pressures for usual wells first to recover the resutls
// without mutlti segment wells. For segment wells, it has not been decided if
// we need th concept of preforation pressures
xw.perfPress().resize(nperf_total, -1.e100);
const V& cdp = stdWells().wellPerforationPressureDiffs();
int start_segment = 0;
int start_perforation = 0;
for (int i = 0; i < nw; ++i) {
WellMultiSegmentConstPtr well = wellsMultiSegment()[i];
const int nperf = well->numberOfPerforations();
const int nseg = well->numberOfSegments();
if (well->isMultiSegmented()) {
start_segment += nseg;
start_perforation += nperf;
continue;
}
const V cdp_well = subset(cdp, Span(nperf, 1, start_perforation));
const ADB segp = subset(state.segp, Span(nseg, 1, start_segment));
const V perfpressure = (well->wellOps().s2p * segp.value().matrix()).array() + cdp_well;
std::copy(perfpressure.data(), perfpressure.data() + nperf, &xw.perfPress()[start_perforation]);
start_segment += nseg;
start_perforation += nperf;
}
wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
}
@ -571,7 +221,7 @@ namespace Opm {
SolutionState& state,
WellState& well_state)
{
const bool converged = baseSolveWellEq(mob_perfcells, b_perfcells, state, well_state);
const bool converged = Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
if (converged) {
// We must now update the state.segp and state.segqs members,
@ -608,150 +258,21 @@ namespace Opm {
template <class Grid>
void
BlackoilMultiSegmentModel<Grid>::updateWellState(const V& dwells,
WellState& well_state)
{
msWells().updateWellState(dwells, dpMaxRel(), well_state);
}
/// added to fixing the flow_multisegment running
template <class Grid>
bool
BlackoilMultiSegmentModel<Grid>::baseSolveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
SolutionState& state,
WellState& well_state) {
V aliveWells;
const int np = wells().number_of_phases;
std::vector<ADB> cq_s(np, ADB::null());
std::vector<int> indices = stdWells().variableWellStateIndices();
SolutionState state0 = state;
WellState well_state0 = well_state;
makeConstantState(state0);
std::vector<ADB> mob_perfcells_const(np, ADB::null());
std::vector<ADB> b_perfcells_const(np, ADB::null());
if ( Base::localWellsActive() ){
// If there are non well in the sudomain of the process
// thene mob_perfcells_const and b_perfcells_const would be empty
for (int phase = 0; phase < np; ++phase) {
mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value());
b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value());
}
}
int it = 0;
bool converged;
do {
// bhp and Q for the wells
std::vector<V> vars0;
vars0.reserve(2);
variableWellStateInitials(well_state, vars0);
std::vector<ADB> vars = ADB::variables(vars0);
SolutionState wellSolutionState = state0;
variableStateExtractWellsVars(indices, vars, wellSolutionState);
const int nw = wellsMultiSegment().size();
const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
const V perf_press_diffs = stdWells().wellPerforationPressureDiffs();
msWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_,
perf_press_diffs, compi,
mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
msWells().addWellFluxEq(cq_s, wellSolutionState, residual_);
// addWellControlEq(wellSolutionState, well_state, aliveWells);
msWells().addWellControlEq(wellSolutionState, well_state, aliveWells, active_, residual_);
converged = Base::getWellConvergence(it);
if (converged) {
break;
}
++it;
if( Base::localWellsActive() )
{
std::vector<ADB> eqs;
eqs.reserve(2);
eqs.push_back(residual_.well_flux_eq);
eqs.push_back(residual_.well_eq);
ADB total_residual = vertcatCollapseJacs(eqs);
const std::vector<M>& Jn = total_residual.derivative();
typedef Eigen::SparseMatrix<double> Sp;
Sp Jn0;
Jn[0].toSparse(Jn0);
const Eigen::SparseLU< Sp > solver(Jn0);
ADB::V total_residual_v = total_residual.value();
const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
assert(dx.size() == total_residual_v.size());
asImpl().updateWellState(dx.array(), well_state);
msWells().updateWellControls(terminal_output_, well_state);
}
} while (it < 15);
if (converged) {
if ( terminal_output_ ) {
std::cout << "well converged iter: " << it << std::endl;
}
const int nw = wells().number_of_wells;
{
// We will set the bhp primary variable to the new ones,
// but we do not change the derivatives here.
ADB::V new_bhp = Eigen::Map<ADB::V>(well_state.bhp().data(), nw);
// Avoiding the copy below would require a value setter method
// in AutoDiffBlock.
std::vector<ADB::M> old_derivs = state.bhp.derivative();
state.bhp = ADB::function(std::move(new_bhp), std::move(old_derivs));
}
{
// Need to reshuffle well rates, from phase running fastest
// to wells running fastest.
// The transpose() below switches the ordering.
const DataBlock wrates = Eigen::Map<const DataBlock>(well_state.wellRates().data(), nw, np).transpose();
ADB::V new_qs = Eigen::Map<const V>(wrates.data(), nw*np);
std::vector<ADB::M> old_derivs = state.qs.derivative();
state.qs = ADB::function(std::move(new_qs), std::move(old_derivs));
}
computeWellConnectionPressures(state, well_state);
}
if (!converged) {
well_state = well_state0;
}
return converged;
}
template <class Grid>
std::vector<V>
BlackoilMultiSegmentModel<Grid>::
variableStateInitials(const ReservoirState& x,
const WellState& xw) const
computeWellConnectionPressures(const SolutionState& state,
const WellState& well_state)
{
assert(active_[ Oil ]);
const int np = x.numPhases();
std::vector<V> vars0;
// p, Sw and Rs, Rv or Sg is used as primary depending on solution conditions
// and bhp and Q for the wells
vars0.reserve(np + 1);
variableReservoirStateInitials(x, vars0);
variableWellStateInitials(xw, vars0);
return vars0;
const int np = numPhases();
const std::vector<ADB> kr_adb = Base::computeRelPerm(state);
std::vector<ADB> fluid_density(np, ADB::null());
// TODO: make sure the order of the density and the order of the kr are the same.
for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
const int canonicalPhaseIdx = canph_[phaseIdx];
fluid_density[phaseIdx] = fluidDensity(canonicalPhaseIdx, rq_[phaseIdx].b, state.rs, state.rv);
}
wellModel().computeWellConnectionPressures(state, well_state, kr_adb, fluid_density);
}
} // namespace Opm

9
opm/autodiff/BlackoilSolventModel.hpp
View File

@ -71,7 +71,7 @@ namespace Opm {
const DerivedGeology& geo,
const RockCompressibility* rock_comp_props,
const SolventPropsAdFromDeck& solvent_props,
const Wells* wells,
const StandardWellsSolvent& well_model,
const NewtonIterationBlackoilInterface& linsolver,
const EclipseStateConstPtr eclState,
const bool has_disgas,
@ -129,7 +129,7 @@ namespace Opm {
// --------- Protected methods ---------
// Need to declare Base members we want to use here.
using Base::stdWells;
using Base::wellModel;
using Base::wells;
using Base::variableState;
using Base::computeGasPressure;
@ -215,11 +215,6 @@ namespace Opm {
const std::vector<PhasePresence>
phaseCondition() const {return this->phaseCondition_;}
void extractWellPerfProperties(const SolutionState& state,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells);
// compute effective viscosities (mu_eff_) and effective b factors (b_eff_) using the ToddLongstaff model
void computeEffectiveProperties(const SolutionState& state);

58
opm/autodiff/BlackoilSolventModel_impl.hpp
View File

@ -76,7 +76,7 @@ namespace Opm {
const DerivedGeology& geo,
const RockCompressibility* rock_comp_props,
const SolventPropsAdFromDeck& solvent_props,
const Wells* wells_arg,
const StandardWellsSolvent& well_model,
const NewtonIterationBlackoilInterface& linsolver,
const EclipseStateConstPtr eclState,
const bool has_disgas,
@ -84,7 +84,7 @@ namespace Opm {
const bool terminal_output,
const bool has_solvent,
const bool is_miscible)
: Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver,
: Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver,
eclState, has_disgas, has_vapoil, terminal_output),
has_solvent_(has_solvent),
solvent_pos_(detail::solventPos(fluid.phaseUsage())),
@ -100,7 +100,8 @@ namespace Opm {
Base::material_name_.push_back("Solvent");
assert(solvent_pos_ == fluid_.numPhases());
residual_.matbalscale.resize(fluid_.numPhases() + 1, 0.0031); // use the same as gas
stdWells().initSolvent(&solvent_props_, solvent_pos_, has_solvent_);
wellModel().initSolvent(&solvent_props_, solvent_pos_, has_solvent_);
}
if (is_miscible_) {
mu_eff_.resize(fluid_.numPhases() + 1, ADB::null());
@ -240,7 +241,7 @@ namespace Opm {
}
}
// wells
Base::variableStateExtractWellsVars(indices, vars, state);
wellModel().variableStateExtractWellsVars(indices, vars, state);
return state;
}
@ -377,6 +378,11 @@ namespace Opm {
}
}
template <class Grid>
void
BlackoilSolventModel<Grid>::
@ -409,7 +415,7 @@ namespace Opm {
varstart += dss.size();
// Extract well parts np phase rates + bhp
const V dwells = subset(dx, Span(Base::numWellVars(), 1, varstart));
const V dwells = subset(dx, Span(wellModel().numWellVars(), 1, varstart));
varstart += dwells.size();
assert(varstart == dx.size());
@ -615,10 +621,7 @@ namespace Opm {
std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
}
// TODO: gravity should be stored as a member
// const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
// asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state);
Base::updateWellState(dwells,well_state);
wellModel().updateWellState(dwells, dpMaxRel(), well_state);
// Update phase conditions used for property calculations.
updatePhaseCondFromPrimalVariable(reservoir_state);
@ -816,43 +819,6 @@ namespace Opm {
}
template <class Grid>
void
BlackoilSolventModel<Grid>::extractWellPerfProperties(const SolutionState& state,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells)
{
Base::extractWellPerfProperties(state, mob_perfcells, b_perfcells);
if (has_solvent_) {
int gas_pos = fluid_.phaseUsage().phase_pos[Gas];
const std::vector<int>& well_cells = stdWells().wellOps().well_cells;
const int nperf = well_cells.size();
// Gas and solvent is combinded and solved together
// The input in the well equation is then the
// total gas phase = hydro carbon gas + solvent gas
// The total mobility is the sum of the solvent and gas mobiliy
mob_perfcells[gas_pos] += subset(rq_[solvent_pos_].mob, well_cells);
// A weighted sum of the b-factors of gas and solvent are used.
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB zero = ADB::constant(V::Zero(nc));
const ADB& ss = state.solvent_saturation;
const ADB& sg = (active_[ Gas ]
? state.saturation[ pu.phase_pos[ Gas ] ]
: zero);
Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
V ones = V::Constant(nperf,1.0);
b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos];
b_perfcells[gas_pos] += (F_solvent * subset(rq_[solvent_pos_].b, well_cells));
}
}
template <class Grid>
void

173
opm/autodiff/MultisegmentWells.cpp
View File

@ -140,10 +140,17 @@ namespace Opm {
MultisegmentWells::
MultisegmentWells(const std::vector<WellMultiSegmentConstPtr>& wells_ms, const int np)
: wells_multisegment_(wells_ms)
, wops_ms_(wells_ms)
, num_phases_(np)
MultisegmentWells(const Wells* wells_arg,
const std::vector<WellConstPtr>& wells_ecl,
const int time_step)
: wells_multisegment_( createMSWellVector(wells_arg, wells_ecl, time_step) )
, wops_ms_(wells_multisegment_)
, num_phases_(wells_arg ? wells_arg->number_of_phases : 0)
, wells_(wells_arg)
, fluid_(nullptr)
, active_(nullptr)
, phase_condition_(nullptr)
, vfp_properties_(nullptr)
, well_segment_perforation_pressure_diffs_(ADB::null())
, well_segment_densities_(ADB::null())
, well_segment_pressures_delta_(ADB::null())
@ -152,16 +159,14 @@ namespace Opm {
, segment_mass_flow_rates_(ADB::null())
, segment_viscosities_(ADB::null())
{
// TODO: repeated with the wellOps's initialization, delete one soon.
// Count the total number of perforations and segments.
const int nw = wells_ms.size();
top_well_segments_.resize(nw);
const int nw = wells_multisegment_.size();
int nperf_total = 0;
int nseg_total = 0;
top_well_segments_.resize(nw);
for (int w = 0; w < nw; ++w) {
top_well_segments_[w] = nseg_total;
nperf_total += wells_ms[w]->numberOfPerforations();
nseg_total += wells_ms[w]->numberOfSegments();
nperf_total += wells_multisegment_[w]->numberOfPerforations();
nseg_total += wells_multisegment_[w]->numberOfSegments();
}
nperf_total_ = nperf_total;
@ -172,8 +177,107 @@ namespace Opm {
std::vector<WellMultiSegmentConstPtr>
MultisegmentWells::createMSWellVector(const Wells* wells_arg,
const std::vector<WellConstPtr>& wells_ecl,
const int time_step)
{
std::vector<WellMultiSegmentConstPtr> wells_multisegment;
if (wells_arg) {
// number of wells in wells_arg structure
const int nw = wells_arg->number_of_wells;
// number of wells in EclipseState
const int nw_ecl = wells_ecl.size();
wells_multisegment.reserve(nw);
for(int i = 0; i < nw_ecl; ++i) {
// not processing SHUT wells
if (wells_ecl[i]->getStatus(time_step) == WellCommon::SHUT) {
continue;
}
// checking if the well can be found in the wells
const std::string& well_name = wells_ecl[i]->name();
int index_well;
for (index_well = 0; index_well < nw; ++index_well) {
if (well_name == std::string(wells_arg->name[index_well])) {
break;
}
}
if (index_well != nw) { // found in the wells
wells_multisegment.push_back(std::make_shared<WellMultiSegment>(wells_ecl[i], time_step, wells_arg));
}
}
}
return wells_multisegment;
}
void
MultisegmentWells::init(const BlackoilPropsAdInterface* fluid_arg,
const std::vector<bool>* active_arg,
const std::vector<PhasePresence>* pc_arg,
const VFPProperties* vfp_properties_arg,
const double gravity_arg,
const Vector& depth_arg)
{
fluid_ = fluid_arg;
active_ = active_arg;
phase_condition_ = pc_arg;
vfp_properties_ = vfp_properties_arg;
gravity_ = gravity_arg;
perf_cell_depth_ = subset(depth_arg, wellOps().well_cells);;
const int nw = wells_multisegment_.size();
// Calculating the depth difference between perforation and the cell center in the peforated cells.
std::vector<double> perf_depth_vec;
perf_depth_vec.reserve(nperf_total_);
for (int w = 0; w < nw; ++w) {
WellMultiSegmentConstPtr well = wells_multisegment_[w];
const std::vector<double>& perf_depth_well = well->perfDepth();
perf_depth_vec.insert(perf_depth_vec.end(), perf_depth_well.begin(), perf_depth_well.end());
}
assert(int(perf_depth_vec.size()) == nperf_total_);
const Vector perf_depth = Eigen::Map<Vector>(perf_depth_vec.data(), nperf_total_);
perf_cell_depth_diffs_ = perf_depth - perf_cell_depth_;
// Calculating the depth difference between segment nodes and perforations.
well_segment_perforation_depth_diffs_ = Vector::Constant(nperf_total_, -1e100);
int start_perforation = 0;
for (int w = 0; w < nw; ++w) {
WellMultiSegmentConstPtr well = wells_multisegment_[w];
const int nseg = well->numberOfSegments();
const int nperf = well->numberOfPerforations();
const std::vector<std::vector<int>>& segment_perforations = well->segmentPerforations();
for (int s = 0; s < nseg; ++s) {
const int nperf_seg = segment_perforations[s].size();
const double segment_depth = well->segmentDepth()[s];
for (int perf = 0; perf < nperf_seg; ++perf) {
const int perf_number = segment_perforations[s][perf] + start_perforation;
well_segment_perforation_depth_diffs_[perf_number] = segment_depth - perf_depth[perf_number];
}
}
start_perforation += nperf;
}
assert(start_perforation == nperf_total_);
}
const std::vector<WellMultiSegmentConstPtr>&
MultisegmentWells::wells() const
MultisegmentWells::msWells() const
{
return wells_multisegment_;
}
@ -182,6 +286,17 @@ namespace Opm {
const Wells&
MultisegmentWells::wells() const
{
assert(wells_ != nullptr);
return *(wells_);
}
const MultisegmentWells::MultisegmentWellOps&
MultisegmentWells::wellOps() const
{
@ -196,7 +311,7 @@ namespace Opm {
MultisegmentWells::
computeSegmentPressuresDelta(const double grav)
{
const int nw = wells().size();
const int nw = msWells().size();
const int nseg_total = nseg_total_;
if ( !wellOps().has_multisegment_wells ) {
@ -210,7 +325,7 @@ namespace Opm {
Vector segment_depth_delta = Vector::Zero(nseg_total);
int start_segment = 0;
for (int w = 0; w < nw; ++w) {
WellMultiSegmentConstPtr well = wells()[w];
WellMultiSegmentConstPtr well = msWells()[w];
const int nseg = well->numberOfSegments();
for (int s = 1; s < nseg; ++s) {
const int s_outlet = well->outletSegment()[s];
@ -228,6 +343,38 @@ namespace Opm {
well_segment_perforation_pressure_diffs_ = grav * well_segment_perforation_depth_diffs_ * well_segment_perforation_densities;
}
void
MultisegmentWells::
variableStateWellIndices(std::vector<int>& indices,
int& next) const
{
indices[Qs] = next++;
indices[Bhp] = next++;
}
std::vector<int>
MultisegmentWells::
variableWellStateIndices() const
{
// Black oil model standard is 5 equation.
// For the pure well solve, only the well equations are picked.
std::vector<int> indices(5, -1);
int next = 0;
variableStateWellIndices(indices, next);
assert(next == 2);
return indices;
}
} // end of namespace Opm

112
opm/autodiff/MultisegmentWells.hpp
View File

@ -37,8 +37,10 @@
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/autodiff/WellHelpers.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/WellMultiSegment.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
@ -80,17 +82,48 @@ namespace Opm {
// --------- Public methods ---------
// TODO: using a vector of WellMultiSegmentConstPtr for now
// TODO: it should use const Wells or something else later.
MultisegmentWells(const std::vector<WellMultiSegmentConstPtr>& wells_multisegment, const int np);
MultisegmentWells(const Wells* wells_arg,
const std::vector<WellConstPtr>& wells_ecl,
const int time_step);
const std::vector<WellMultiSegmentConstPtr>& wells() const;
std::vector<WellMultiSegmentConstPtr> createMSWellVector(const Wells* wells_arg,
const std::vector<WellConstPtr>& wells_ecl,
const int time_step);
void init(const BlackoilPropsAdInterface* fluid_arg,
const std::vector<bool>* active_arg,
const std::vector<PhasePresence>* pc_arg,
const VFPProperties* vfp_properties_arg,
const double gravity_arg,
const Vector& depth_arg);
const std::vector<WellMultiSegmentConstPtr>& msWells() const;
const MultisegmentWellOps& wellOps() const;
const Wells& wells() const;
int numPhases() const { return num_phases_; };
int numWells() const { return wells_multisegment_.size(); }
int numSegment() const { return nseg_total_; };
int numPerf() const { return nperf_total_; };
bool wellsActive() const { return wells_active_; };
void setWellsActive(const bool wells_active) { wells_active_ = wells_active; };
bool localWellsActive() const { return ! wells_multisegment_.empty(); };
int numWellVars() const { return (num_phases_ + 1) * nseg_total_; };
template <class ReservoirResidualQuant, class SolutionState>
void extractWellPerfProperties(const SolutionState& state,
const std::vector<ReservoirResidualQuant>& rq,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const;
Vector& wellPerforationCellPressureDiffs() { return well_perforation_cell_pressure_diffs_; };
Vector& wellSegmentPerforationDepthDiffs() { return well_segment_perforation_depth_diffs_; };
@ -108,7 +141,11 @@ namespace Opm {
Vector& segVDt() { return segvdt_; };
const Vector& wellPerforationDensities() const { return well_perforation_densities_; };
Vector& wellPerforationDensities() { return well_perforation_densities_; };
const Vector& wellPerforationPressureDiffs() const { return well_perforation_pressure_diffs_; };
Vector& wellPerforationPressureDiffs() { return well_perforation_pressure_diffs_; };
template <class WellState>
void
@ -121,24 +158,22 @@ namespace Opm {
template <class SolutionState>
void
computeWellFlux(const SolutionState& state,
const Opm::PhaseUsage& pu,
const std::vector<bool>& active,
const Vector& well_perforation_pressure_diffs,
const DataBlock& compi,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
Vector& aliveWells,
std::vector<ADB>& cq_s) const;
template <class SolutionState, class WellState>
void updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw) const;
// Calculate the density of the mixture in the segments
// And the surface volume of the components in the segments by dt
template <class SolutionState>
void
computeSegmentFluidProperties(const SolutionState& state,
const std::vector<PhasePresence>& pc,
const std::vector<bool>& active,
const BlackoilPropsAdInterface& fluid);
computeSegmentFluidProperties(const SolutionState& state);
void
computeSegmentPressuresDelta(const double grav);
@ -154,7 +189,6 @@ namespace Opm {
addWellControlEq(const SolutionState& state,
const WellState& xw,
const Vector& aliveWells,
const std::vector<bool>& active,
LinearisedBlackoilResidual& residual);
template <class WellState>
@ -162,14 +196,56 @@ namespace Opm {
updateWellControls(const bool terminal_output,
WellState& xw) const;
// TODO: these code are same with the StandardWells
// to find a better solution later.
void
variableStateWellIndices(std::vector<int>& indices,
int& next) const;
template <class SolutionState>
void
variableStateExtractWellsVars(const std::vector<int>& indices,
std::vector<ADB>& vars,
SolutionState& state) const;
std::vector<int>
variableWellStateIndices() const;
template <class WellState>
void
variableWellStateInitials(const WellState& xw,
std::vector<Vector>& vars0) const;
template <class SolutionState, class WellState>
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const std::vector<ADB>& kr_adb,
const std::vector<ADB>& fluid_density);
protected:
// TODO: probably a wells_active_ will be required here.
const std::vector<WellMultiSegmentConstPtr> wells_multisegment_;
const MultisegmentWellOps wops_ms_;
bool wells_active_;
std::vector<WellMultiSegmentConstPtr> wells_multisegment_;
MultisegmentWellOps wops_ms_;
const int num_phases_;
int nseg_total_;
int nperf_total_;
// TODO: put the Wells here to simplify the interface
// TODO: at the moment, both wells_ and wells_multisegment_
// TODO: acutally contain all the wells
// TODO: they should be split eventually.
const Wells* wells_;
const BlackoilPropsAdInterface* fluid_;
const std::vector<bool>* active_;
const std::vector<PhasePresence>* phase_condition_;
const VFPProperties* vfp_properties_;
double gravity_;
// The depth of the all the cell centers
// It is different from the perforation depth in MultisegmentWells
Vector perf_cell_depth_;
// Pressure correction due to the different depth of the perforation
// and the cell center of the grid block
// For the non-segmented wells, since the perforation are forced to be
@ -183,6 +259,9 @@ namespace Opm {
// The depth difference between segment nodes and perforations
Vector well_segment_perforation_depth_diffs_;
// The depth difference between the perforations and the perforation cells.
Vector perf_cell_depth_diffs_;
// the average of the fluid densities in the grid block
// which is used to calculate the hydrostatic head correction due to the depth difference of the perforation
// and the cell center of the grid block
@ -218,6 +297,13 @@ namespace Opm {
// to handle the volume effects of the segment
Vector segvdt_;
// technically, they are only useful for standard wells
// since at the moment, we are handling both the standard
// wells and the multi-segment wells under the MultisegmentWells
// we need them to remove the dependency on StandardWells
Vector well_perforation_densities_;
Vector well_perforation_pressure_diffs_;
};
} // namespace Opm

432
opm/autodiff/MultisegmentWells_impl.hpp
View File

@ -47,6 +47,37 @@ namespace Opm
template <class ReservoirResidualQuant, class SolutionState>
void
MultisegmentWells::
extractWellPerfProperties(const SolutionState& /* state */,
const std::vector<ReservoirResidualQuant>& rq,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const
{
// If we have wells, extract the mobilities and b-factors for
// the well-perforated cells.
if ( !localWellsActive() ) {
mob_perfcells.clear();
b_perfcells.clear();
return;
} else {
const std::vector<int>& well_cells = wellOps().well_cells;
mob_perfcells.resize(num_phases_, ADB::null());
b_perfcells.resize(num_phases_, ADB::null());
for (int phase = 0; phase < num_phases_; ++phase) {
mob_perfcells[phase] = subset(rq[phase].mob, well_cells);
b_perfcells[phase] = subset(rq[phase].b, well_cells);
}
}
}
template <class WellState>
void
MultisegmentWells::
@ -54,9 +85,9 @@ namespace Opm
const double dpmaxrel,
WellState& well_state) const
{
if (!wells().empty())
if (!msWells().empty())
{
const int nw = wells().size();
const int nw = msWells().size();
const int nseg_total = nseg_total_;
const int np = numPhases();
@ -102,7 +133,7 @@ namespace Opm
wr[p + np * w] = well_state.segPhaseRates()[p + np * start_segment];
}
const int nseg = wells()[w]->numberOfSegments();
const int nseg = msWells()[w]->numberOfSegments();
start_segment += nseg;
}
@ -122,25 +153,23 @@ namespace Opm
void
MultisegmentWells::
computeWellFlux(const SolutionState& state,
const Opm::PhaseUsage& pu,
const std::vector<bool>& active,
const Vector& well_perforation_pressure_diffs,
const DataBlock& compi,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
Vector& aliveWells,
std::vector<ADB>& cq_s) const
{
if (wells().size() == 0) return;
if (msWells().size() == 0) return;
const int np = numPhases();
const int nw = wells().size();
const int nw = msWells().size();
aliveWells = Vector::Constant(nw, 1.0);
const int nseg = nseg_total_;
const int nperf = nperf_total_;
const Opm::PhaseUsage& pu = fluid_->phaseUsage();
cq_s.resize(np, ADB::null());
{
@ -160,7 +189,7 @@ namespace Opm
// Create selector for perforations of multi-segment vs. regular wells.
Vector is_multisegment_well(nw);
for (int w = 0; w < nw; ++w) {
is_multisegment_well[w] = double(wells()[w]->isMultiSegmented());
is_multisegment_well[w] = double(msWells()[w]->isMultiSegmented());
}
// Take one flag per well and expand to one flag per perforation.
Vector is_multisegment_perf = wellOps().w2p * is_multisegment_well.matrix();
@ -168,7 +197,7 @@ namespace Opm
// Compute drawdown.
ADB h_nc = msperf_selector.select(well_segment_perforation_pressure_diffs_,
ADB::constant(well_perforation_pressure_diffs));
ADB::constant(well_perforation_pressure_diffs_));
const Vector h_cj = msperf_selector.select(well_perforation_cell_pressure_diffs_, Vector::Zero(nperf));
// Special handling for when we are called from solveWellEq().
@ -201,7 +230,7 @@ namespace Opm
cq_ps[phase] = b_perfcells[phase] * cq_p;
}
if (active[Oil] && active[Gas]) {
if ((*active_)[Oil] && (*active_)[Gas]) {
const int oilpos = pu.phase_pos[Oil];
const int gaspos = pu.phase_pos[Gas];
const ADB cq_psOil = cq_ps[oilpos];
@ -234,10 +263,12 @@ namespace Opm
// TODO: involves one operations that are not valid now. (i.e. how to transverse from the leaves to the root,
// TODO: although we can begin from the brutal force way)
// TODO: stop using wells() here.
// TODO: stop using msWells() here.
std::vector<ADB> wbq(np, ADB::null());
ADB wbqt = ADB::constant(Vector::Zero(nseg));
const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
for (int phase = 0; phase < np; ++phase) {
const ADB& q_ps = wellOps().p2s * cq_ps[phase];
const ADB& q_s = subset(state.segqs, Span(nseg, 1, phase * nseg));
@ -263,7 +294,7 @@ namespace Opm
if (wbqt.value()[topseg] == 0.0) { // yes we really mean == here, no fuzzyness
aliveWells[w] = 0.0;
}
topseg += wells()[w]->numberOfSegments();
topseg += msWells()[w]->numberOfSegments();
}
}
@ -285,13 +316,13 @@ namespace Opm
for (int phase = 0; phase < np; ++phase) {
ADB tmp = cmix_s[phase];
if (phase == Oil && active[Gas]) {
if (phase == Oil && (*active_)[Gas]) {
const int gaspos = pu.phase_pos[Gas];
tmp = tmp - rv_perfcells * cmix_s[gaspos] / d;
tmp = (tmp - rv_perfcells * cmix_s[gaspos]) / d;
}
if (phase == Gas && active[Oil]) {
if (phase == Gas && (*active_)[Oil]) {
const int oilpos = pu.phase_pos[Oil];
tmp = tmp - rs_perfcells * cmix_s[oilpos] / d;
tmp = (tmp - rs_perfcells * cmix_s[oilpos]) / d;
}
volumeRatio += tmp / b_perfcells[phase];
}
@ -310,16 +341,63 @@ namespace Opm
template <class SolutionState, class WellState>
void
MultisegmentWells::
updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw) const
{
// Update the perforation phase rates (used to calculate the pressure drop in the wellbore).
const int np = numPhases();
const int nw = numWells();
Vector cq = superset(cq_s[0].value(), Span(nperf_total_, np, 0), nperf_total_ * np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(cq_s[phase].value(), Span(nperf_total_, np, phase), nperf_total_ * np);
}
xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf_total_ * np);
// Update the perforation pressures for usual wells first to recover the resutls
// without mutlti segment wells. For segment wells, it has not been decided if
// we need th concept of preforation pressures
xw.perfPress().resize(nperf_total_, -1.e100);
const Vector& cdp = well_perforation_pressure_diffs_;
int start_segment = 0;
int start_perforation = 0;
for (int i = 0; i < nw; ++i) {
WellMultiSegmentConstPtr well = wells_multisegment_[i];
const int nperf = well->numberOfPerforations();
const int nseg = well->numberOfSegments();
if (well->isMultiSegmented()) {
start_segment += nseg;
start_perforation += nperf;
continue;
}
const Vector cdp_well = subset(cdp, Span(nperf, 1, start_perforation));
const ADB segp = subset(state.segp, Span(nseg, 1, start_segment));
const Vector perfpressure = (well->wellOps().s2p * segp.value().matrix()).array() + cdp_well;
std::copy(perfpressure.data(), perfpressure.data() + nperf, &xw.perfPress()[start_perforation]);
start_segment += nseg;
start_perforation += nperf;
}
assert(start_segment == nseg_total_);
assert(start_perforation == nperf_total_);
}
template <class SolutionState>
void
MultisegmentWells::
computeSegmentFluidProperties(const SolutionState& state,
const std::vector<PhasePresence>& pc,
const std::vector<bool>& active,
const BlackoilPropsAdInterface& fluid)
computeSegmentFluidProperties(const SolutionState& state)
{
const int np = numPhases();
const int nw = wells().size();
const int nw = msWells().size();
const int nseg_total = nseg_total_;
if ( !wellOps().has_multisegment_wells ){
@ -347,7 +425,7 @@ namespace Opm
std::vector<int> segment_cells;
segment_cells.reserve(nseg_total);
for (int w = 0; w < nw; ++w) {
const std::vector<int>& segment_cells_well = wells()[w]->segmentCells();
const std::vector<int>& segment_cells_well = msWells()[w]->segmentCells();
segment_cells.insert(segment_cells.end(), segment_cells_well.begin(), segment_cells_well.end());
}
assert(int(segment_cells.size()) == nseg_total);
@ -360,37 +438,37 @@ namespace Opm
// Compute PVT properties for segments.
std::vector<PhasePresence> segment_cond(nseg_total);
for (int s = 0; s < nseg_total; ++s) {
segment_cond[s] = pc[segment_cells[s]];
segment_cond[s] = (*phase_condition_)[segment_cells[s]];
}
std::vector<ADB> b_seg(np, ADB::null());
// Viscosities for different phases
std::vector<ADB> mu_seg(np, ADB::null());
ADB rsmax_seg = ADB::null();
ADB rvmax_seg = ADB::null();
const PhaseUsage& pu = fluid.phaseUsage();
const PhaseUsage& pu = fluid_->phaseUsage();
if (pu.phase_used[Water]) {
b_seg[pu.phase_pos[Water]] = fluid.bWat(segment_press, segment_temp, segment_cells);
mu_seg[pu.phase_pos[Water]] = fluid.muWat(segment_press, segment_temp, segment_cells);
b_seg[pu.phase_pos[Water]] = fluid_->bWat(segment_press, segment_temp, segment_cells);
mu_seg[pu.phase_pos[Water]] = fluid_->muWat(segment_press, segment_temp, segment_cells);
}
assert(active[Oil]);
assert((*active_)[Oil]);
const ADB segment_so = subset(state.saturation[pu.phase_pos[Oil]], segment_cells);
if (pu.phase_used[Oil]) {
const ADB segment_rs = subset(state.rs, segment_cells);
b_seg[pu.phase_pos[Oil]] = fluid.bOil(segment_press, segment_temp, segment_rs,
b_seg[pu.phase_pos[Oil]] = fluid_->bOil(segment_press, segment_temp, segment_rs,
segment_cond, segment_cells);
// rsmax_seg = fluidRsSat(segment_press, segment_so, segment_cells);
rsmax_seg = fluid.rsSat(segment_press, segment_so, segment_cells);
mu_seg[pu.phase_pos[Oil]] = fluid.muOil(segment_press, segment_temp, segment_rs,
rsmax_seg = fluid_->rsSat(segment_press, segment_so, segment_cells);
mu_seg[pu.phase_pos[Oil]] = fluid_->muOil(segment_press, segment_temp, segment_rs,
segment_cond, segment_cells);
}
assert(active[Gas]);
assert((*active_)[Gas]);
if (pu.phase_used[Gas]) {
const ADB segment_rv = subset(state.rv, segment_cells);
b_seg[pu.phase_pos[Gas]] = fluid.bGas(segment_press, segment_temp, segment_rv,
b_seg[pu.phase_pos[Gas]] = fluid_->bGas(segment_press, segment_temp, segment_rv,
segment_cond, segment_cells);
// rvmax_seg = fluidRvSat(segment_press, segment_so, segment_cells);
rvmax_seg = fluid.rvSat(segment_press, segment_so, segment_cells);
mu_seg[pu.phase_pos[Gas]] = fluid.muGas(segment_press, segment_temp, segment_rv,
rvmax_seg = fluid_->rvSat(segment_press, segment_so, segment_cells);
mu_seg[pu.phase_pos[Gas]] = fluid_->muGas(segment_press, segment_temp, segment_rv,
segment_cond, segment_cells);
}
@ -406,7 +484,7 @@ namespace Opm
std::vector<std::vector<double>> comp_frac(np, std::vector<double>(nseg_total, 0.0));
int start_segment = 0;
for (int w = 0; w < nw; ++w) {
WellMultiSegmentConstPtr well = wells()[w];
WellMultiSegmentConstPtr well = msWells()[w];
const int nseg = well->numberOfSegments();
const std::vector<double>& comp_frac_well = well->compFrac();
for (int phase = 0; phase < np; ++phase) {
@ -444,7 +522,7 @@ namespace Opm
ADB big_values = ADB::constant(Vector::Constant(nseg_total, 1.e100));
ADB mix_gas_oil = non_zero_mix_oilpos.select(mix[gaspos] / mix[oilpos], big_values);
ADB mix_oil_gas = non_zero_mix_gaspos.select(mix[oilpos] / mix[gaspos], big_values);
if (active[Oil]) {
if ((*active_)[Oil]) {
Vector selectorUnderRsmax = Vector::Zero(nseg_total);
Vector selectorAboveRsmax = Vector::Zero(nseg_total);
for (int s = 0; s < nseg_total; ++s) {
@ -456,7 +534,7 @@ namespace Opm
}
rs = non_zero_mix_oilpos.select(selectorAboveRsmax * rsmax_seg + selectorUnderRsmax * mix_gas_oil, rs);
}
if (active[Gas]) {
if ((*active_)[Gas]) {
Vector selectorUnderRvmax = Vector::Zero(nseg_total);
Vector selectorAboveRvmax = Vector::Zero(nseg_total);
for (int s = 0; s < nseg_total; ++s) {
@ -474,7 +552,7 @@ namespace Opm
for (int phase = 0; phase < np; ++phase) {
x[phase] = mix[phase];
}
if (active[Gas] && active[Oil]) {
if ((*active_)[Gas] && (*active_)[Oil]) {
x[gaspos] = (mix[gaspos] - mix[oilpos] * rs) / (Vector::Ones(nseg_total) - rs * rv);
x[oilpos] = (mix[oilpos] - mix[gaspos] * rv) / (Vector::Ones(nseg_total) - rs * rv);
}
@ -488,7 +566,7 @@ namespace Opm
// Compute segment densities.
ADB dens = ADB::constant(Vector::Zero(nseg_total));
for (int phase = 0; phase < np; ++phase) {
const Vector surface_density = fluid.surfaceDensity(phase, segment_cells);
const Vector surface_density = fluid_->surfaceDensity(phase, segment_cells);
dens += surface_density * mix[phase];
}
well_segment_densities_ = dens / volrat;
@ -504,7 +582,7 @@ namespace Opm
segment_mass_flow_rates_ = ADB::constant(Vector::Zero(nseg_total));
for (int phase = 0; phase < np; ++phase) {
// TODO: how to remove one repeated surfaceDensity()
const Vector surface_density = fluid.surfaceDensity(phase, segment_cells);
const Vector surface_density = fluid_->surfaceDensity(phase, segment_cells);
segment_mass_flow_rates_ += surface_density * segqs[phase];
}
@ -580,29 +658,28 @@ namespace Opm
addWellControlEq(const SolutionState& state,
const WellState& xw,
const Vector& aliveWells,
const std::vector<bool>& active,
LinearisedBlackoilResidual& residual)
{
// the name of the function is a a little misleading.
// Basically it is the function for the pressure equation.
// And also, it work as the control equation when it is the segment
if( wells().empty() ) return;
if( msWells().empty() ) return;
const int np = numPhases();
const int nw = wells().size();
const int nw = msWells().size();
const int nseg_total = nseg_total_;
ADB aqua = ADB::constant(Vector::Zero(nseg_total));
ADB liquid = ADB::constant(Vector::Zero(nseg_total));
ADB vapour = ADB::constant(Vector::Zero(nseg_total));
if (active[Water]) {
if ((*active_)[Water]) {
aqua += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Aqua * nseg_total));
}
if (active[Oil]) {
if ((*active_)[Oil]) {
liquid += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Liquid * nseg_total));
}
if (active[Gas]) {
if ((*active_)[Gas]) {
vapour += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Vapour * nseg_total));
}
@ -631,14 +708,14 @@ namespace Opm
//and gather info about current control
int start_segment = 0;
for (int w = 0; w < nw; ++w) {
const struct WellControls* wc = wells()[w]->wellControls();
const struct WellControls* wc = msWells()[w]->wellControls();
// 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];
const int nseg = wells()[w]->numberOfSegments();
const int nseg = msWells()[w]->numberOfSegments();
switch (well_controls_iget_type(wc, current)) {
case BHP:
@ -748,15 +825,15 @@ namespace Opm
updateWellControls(const bool terminal_output,
WellState& xw) const
{
if( wells().empty() ) return ;
if( msWells().empty() ) return ;
std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" };
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
const int np = numPhases();
const int nw = wells().size();
const int nw = msWells().size();
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells()[w]->wellControls();
const WellControls* wc = msWells()[w]->wellControls();
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
@ -775,7 +852,7 @@ namespace Opm
}
if (wellhelpers::constraintBroken(
xw.bhp(), xw.thp(), xw.wellRates(),
w, np, wells()[w]->wellType(), wc, ctrl_index)) {
w, np, msWells()[w]->wellType(), wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
}
@ -785,7 +862,7 @@ namespace Opm
// Constraint number ctrl_index was broken, switch to it.
if (terminal_output)
{
std::cout << "Switching control mode for well " << wells()[w]->name()
std::cout << "Switching control mode for well " << msWells()[w]->name()
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
}
@ -834,5 +911,250 @@ namespace Opm
}
// TODO: This is just a preliminary version, remains to be improved later when we decide a better way
// TODO: to intergrate the usual wells and multi-segment wells.
template <class SolutionState, class WellState>
void
MultisegmentWells::
computeWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const std::vector<ADB>& kr_adb,
const std::vector<ADB>& fluid_density)
{
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_total = nperf_total_;
const int nw = numWells();
const std::vector<int>& well_cells = wellOps().well_cells;
well_perforation_densities_ = Vector::Zero(nperf_total);
const Vector perf_press = Eigen::Map<const Vector>(xw.perfPress().data(), nperf_total);
Vector avg_press = perf_press * 0.0;
// for the non-segmented/regular wells, calculated the average pressures.
// If it is the top perforation, then average with the bhp().
// If it is not the top perforation, then average with the perforation above it().
int start_segment = 0;
for (int w = 0; w < nw; ++w) {
const int nseg = wells_multisegment_[w]->numberOfSegments();
if (wells_multisegment_[w]->isMultiSegmented()) {
// maybe we should give some reasonable values to prevent the following calculations fail
start_segment += nseg;
continue;
}
std::string well_name(wells_multisegment_[w]->name());
typedef typename WellState::SegmentedWellMapType::const_iterator const_iterator;
const_iterator it_well = xw.segmentedWellMap().find(well_name);
assert(it_well != xw.segmentedWellMap().end());
const int start_perforation = (*it_well).second.start_perforation;
const int end_perforation = start_perforation + (*it_well).second.number_of_perforations;
for (int perf = start_perforation; perf < end_perforation; ++perf) {
const double p_above = perf == start_perforation ? state.segp.value()[start_segment] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
start_segment += nseg;
}
assert(start_segment == xw.numSegments());
// 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_total);
for (int perf = 0; perf < nperf_total; ++perf) {
perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
}
const PhaseUsage& pu = fluid_->phaseUsage();
DataBlock b(nperf_total, pu.num_phases);
std::vector<double> rsmax_perf(nperf_total, 0.0);
std::vector<double> rvmax_perf(nperf_total, 0.0);
if (pu.phase_used[BlackoilPhases::Aqua]) {
const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
}
assert((*active_)[Oil]);
const Vector 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 Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rsmax_perf.assign(rssat.data(), rssat.data() + nperf_total);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = subset(state.rv, well_cells);
const Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf_total);
}
// b is row major, so can just copy data.
std::vector<double> b_perf(b.data(), b.data() + nperf_total * pu.num_phases);
const Vector& perfcelldepth = perf_cell_depth_;
std::vector<double> perf_cell_depth(perfcelldepth.data(), perfcelldepth.data() + nperf_total);
// Surface density.
// The compute density segment wants the surface densities as
// an np * number of wells cells array
Vector rho = superset(fluid_->surfaceDensity(0 , well_cells), Span(nperf_total, pu.num_phases, 0), nperf_total * pu.num_phases);
for (int phase = 1; phase < pu.num_phases; ++phase) {
rho += superset(fluid_->surfaceDensity(phase , well_cells), Span(nperf_total, pu.num_phases, phase), nperf_total * pu.num_phases);
}
std::vector<double> surf_dens_perf(rho.data(), rho.data() + nperf_total * pu.num_phases);
// 2. Compute densities
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_->phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// 3. Compute pressure deltas
std::vector<double> cdp =
WellDensitySegmented::computeConnectionPressureDelta(
wells(), perf_cell_depth, cd, gravity_);
// 4. Store the results
well_perforation_densities_ = Eigen::Map<const Vector>(cd.data(), nperf_total); // This one is not useful for segmented wells at all
well_perforation_pressure_diffs_ = Eigen::Map<const Vector>(cdp.data(), nperf_total);
if ( !wellOps().has_multisegment_wells ) {
well_perforation_cell_densities_ = Vector::Zero(nperf_total);
well_perforation_cell_pressure_diffs_ = Vector::Zero(nperf_total);
return;
}
// compute the average of the fluid densites in the well blocks.
// the average is weighted according to the fluid relative permeabilities.
// const std::vector<ADB> kr_adb = Base::computeRelPerm(state);
size_t temp_size = kr_adb.size();
std::vector<Vector> perf_kr;
for(size_t i = 0; i < temp_size; ++i) {
// const ADB kr_phase_adb = subset(kr_adb[i], well_cells);
const Vector kr_phase = (subset(kr_adb[i], well_cells)).value();
perf_kr.push_back(kr_phase);
}
// compute the averaged density for the well block
// TODO: for the non-segmented wells, they should be set to zero
// TODO: for the moment, they are still calculated, while not used later.
for (int i = 0; i < nperf_total; ++i) {
double sum_kr = 0.;
int np = perf_kr.size(); // make sure it is 3
for (int p = 0; p < np; ++p) {
sum_kr += perf_kr[p][i];
}
for (int p = 0; p < np; ++p) {
perf_kr[p][i] /= sum_kr;
}
}
Vector rho_avg_perf = Vector::Constant(nperf_total, 0.0);
// TODO: make sure the order of the density and the order of the kr are the same.
for (int phaseIdx = 0; phaseIdx < fluid_->numPhases(); ++phaseIdx) {
// const int canonicalPhaseIdx = canph_[phaseIdx];
// const ADB fluid_density = fluidDensity(canonicalPhaseIdx, rq_[phaseIdx].b, state.rs, state.rv);
const Vector rho_perf = subset(fluid_density[phaseIdx], well_cells).value();
// TODO: phaseIdx or canonicalPhaseIdx ?
rho_avg_perf += rho_perf * perf_kr[phaseIdx];
}
well_perforation_cell_densities_ = Eigen::Map<const Vector>(rho_avg_perf.data(), nperf_total);
well_perforation_cell_pressure_diffs_ = gravity_ * well_perforation_cell_densities_ * perf_cell_depth_diffs_;
}
template <class SolutionState>
void
MultisegmentWells::
variableStateExtractWellsVars(const std::vector<int>& indices,
std::vector<ADB>& vars,
SolutionState& state) const
{
// TODO: using the original Qs for the segment rates for now, to be fixed eventually.
// TODO: using the original Bhp for the segment pressures for now, to be fixed eventually.
// segment phase rates in surface volume
state.segqs = std::move(vars[indices[Qs]]);
// segment pressures
state.segp = std::move(vars[indices[Bhp]]);
// The qs and bhp are no longer primary variables, but could
// still be used in computations. They are identical to the
// pressures and flows of the top segments.
const int np = num_phases_;
const int ns = nseg_total_;
const int nw = numWells();
state.qs = ADB::constant(Vector::Zero(np * nw));
for (int phase = 0; phase < np; ++phase) {
// Extract segment fluxes for this phase (ns consecutive elements).
ADB segqs_phase = subset(state.segqs, Span(ns, 1, ns*phase));
// Extract top segment fluxes (= well fluxes)
ADB wellqs_phase = subset(segqs_phase, topWellSegments());
// Expand to full size of qs (which contains all phases) and add.
state.qs += superset(wellqs_phase, Span(nw, 1, nw * phase), nw * np);
}
state.bhp = subset(state.segp, topWellSegments());
}
template <class WellState>
void
MultisegmentWells::
variableWellStateInitials(const WellState& xw,
std::vector<Vector>& vars0) const
{
// Initial well rates
if ( wells_multisegment_.size() > 0 )
{
// Need to reshuffle well segment rates, from phase running fastest
const int nseg = xw.numSegments();
const int np = xw.numPhases();
// The transpose() below switches the ordering of the segment rates
const DataBlock segrates = Eigen::Map<const DataBlock>(& xw.segPhaseRates()[0], nseg, np).transpose();
// segment phase rates in surface volume
const Vector segqs = Eigen::Map<const Vector>(segrates.data(), nseg * np);
vars0.push_back(segqs);
// for the pressure of the segments
const Vector segp = Eigen::Map<const Vector>(& xw.segPress()[0], xw.segPress().size());
vars0.push_back(segp);
}
else
{
// push null sates for segqs and segp
vars0.push_back(Vector());
vars0.push_back(Vector());
}
}
}
#endif // OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED

3
opm/autodiff/SimulatorBase.hpp
View File

@ -90,6 +90,7 @@ namespace Opm
typedef typename Traits::OutputWriter OutputWriter;
typedef typename Traits::Grid Grid;
typedef typename Traits::Solver Solver;
typedef typename Traits::WellModel WellModel;
/// Initialise from parameters and objects to observe.
/// \param[in] param parameters, this class accepts the following:
@ -149,7 +150,7 @@ namespace Opm
WellState& well_state,
const Wells* wells);
std::unique_ptr<Solver> createSolver(const Wells* wells);
std::unique_ptr<Solver> createSolver(const WellModel& well_model);
void
computeRESV(const std::size_t step,

8
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -169,7 +169,9 @@ namespace Opm
// Run a multiple steps of the solver depending on the time step control.
solver_timer.start();
auto solver = asImpl().createSolver(wells);
const WellModel well_model(wells);
auto solver = asImpl().createSolver(well_model);
// If sub stepping is enabled allow the solver to sub cycle
// in case the report steps are too large for the solver to converge
@ -373,7 +375,7 @@ namespace Opm
{ }
template <class Implementation>
auto SimulatorBase<Implementation>::createSolver(const Wells* wells)
auto SimulatorBase<Implementation>::createSolver(const WellModel& well_model)
-> std::unique_ptr<Solver>
{
auto model = std::unique_ptr<Model>(new Model(model_param_,
@ -381,7 +383,7 @@ namespace Opm
props_,
geo_,
rock_comp_props_,
wells,
well_model,
solver_,
eclipse_state_,
has_disgas_,

2
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
View File

@ -28,6 +28,7 @@ namespace Opm {
template <class GridT>
class SimulatorFullyImplicitBlackoil;
class StandardWells;
template <class GridT>
struct SimulatorTraits<SimulatorFullyImplicitBlackoil<GridT> >
@ -38,6 +39,7 @@ struct SimulatorTraits<SimulatorFullyImplicitBlackoil<GridT> >
typedef GridT Grid;
typedef BlackoilModel<Grid> Model;
typedef NonlinearSolver<Model> Solver;
typedef StandardWells WellModel;
};
/// a simulator for the blackoil model

6
opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp
View File

@ -33,6 +33,8 @@ namespace Opm {
template <class GridT>
class SimulatorFullyImplicitBlackoilMultiSegment;
class MultisegmentWells;
template <class GridT>
struct SimulatorTraits<SimulatorFullyImplicitBlackoilMultiSegment<GridT> >
{
@ -42,6 +44,7 @@ struct SimulatorTraits<SimulatorFullyImplicitBlackoilMultiSegment<GridT> >
typedef GridT Grid;
typedef BlackoilMultiSegmentModel<Grid> Model;
typedef NonlinearSolver<Model> Solver;
typedef MultisegmentWells WellModel;
};
/// a simulator for the blackoil model
@ -56,6 +59,7 @@ public:
typedef typename Traits::ReservoirState ReservoirState;
typedef typename Traits::WellState WellState;
typedef typename Traits::Solver Solver;
typedef typename Traits::WellModel WellModel;
// forward the constructor to the base class
SimulatorFullyImplicitBlackoilMultiSegment(const parameter::ParameterGroup& param,
@ -80,7 +84,7 @@ public:
protected:
std::unique_ptr<Solver> createSolver(const Wells* wells, std::vector<WellMultiSegmentConstPtr>& wells_multisegment);
std::unique_ptr<Solver> createSolver(const WellModel& well_model);
using Base::output_writer_;
using Base::param_;

34
opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp
View File

@ -25,7 +25,7 @@ namespace Opm
template <class GridT>
auto SimulatorFullyImplicitBlackoilMultiSegment<GridT>::
createSolver(const Wells* wells, std::vector<WellMultiSegmentConstPtr>& wells_multisegment)
createSolver(const WellModel& well_model)
-> std::unique_ptr<Solver>
{
typedef typename Traits::Model Model;
@ -35,13 +35,12 @@ namespace Opm
props_,
geo_,
rock_comp_props_,
wells,
well_model,
solver_,
eclipse_state_,
has_disgas_,
has_vapoil_,
terminal_output_,
wells_multisegment));
terminal_output_));
if (!Base::threshold_pressures_by_face_.empty()) {
model->setThresholdPressures(Base::threshold_pressures_by_face_);
@ -112,30 +111,11 @@ namespace Opm
// well_state.init(wells, state, prev_well_state);
const std::vector<WellConstPtr>& wells_ecl = eclipse_state_->getSchedule()->getWells(timer.currentStepNum());
std::vector<WellMultiSegmentConstPtr> wells_multisegment;
wells_multisegment.reserve(wells_ecl.size());
for (size_t i = 0; i < wells_ecl.size(); ++i) {
// not processing SHUT wells.
if (wells_ecl[i]->getStatus(timer.currentStepNum()) == WellCommon::SHUT) {
continue;
}
// checking if the well can be found in the wells
const std::string& well_name = wells_ecl[i]->name();
// number of wells in wells
const int nw_wells = wells->number_of_wells;
int index_well;
for (index_well = 0; index_well < nw_wells; ++index_well) {
if (well_name == std::string(wells->name[index_well])) {
break;
}
}
const int current_time_step = timer.currentStepNum();
if (index_well != nw_wells) { // found in the wells
wells_multisegment.push_back(std::make_shared<WellMultiSegment>(wells_ecl[i], timer.currentStepNum(), wells));
}
}
const WellModel well_model(wells, wells_ecl, current_time_step);
well_state.init(wells_multisegment, state, prev_well_state);
well_state.init(well_model, state, prev_well_state);
// give the polymer and surfactant simulators the chance to do their stuff
Base::asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
@ -153,7 +133,7 @@ namespace Opm
// Run a multiple steps of the solver depending on the time step control.
solver_timer.start();
auto solver = createSolver(wells, wells_multisegment);
auto solver = createSolver(well_model);
// If sub stepping is enabled allow the solver to sub cycle
// in case the report steps are too large for the solver to converge

7
opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp
View File

@ -80,6 +80,8 @@ namespace Opm
template <class GridT>
class SimulatorFullyImplicitBlackoilSolvent;
class StandardWellsSolvent;
template<class GridT>
struct SimulatorTraits<SimulatorFullyImplicitBlackoilSolvent<GridT> >
{
@ -89,6 +91,8 @@ namespace Opm
typedef GridT Grid;
typedef BlackoilSolventModel<Grid> Model;
typedef NonlinearSolver<Model> Solver;
typedef StandardWellsSolvent WellModel;
};
/// Class collecting all necessary components for a blackoil simulation with polymer
@ -102,6 +106,7 @@ namespace Opm
typedef SimulatorTraits<ThisType> Traits;
typedef typename Traits::Solver Solver;
typedef typename Traits::WellModel WellModel;
public:
SimulatorFullyImplicitBlackoilSolvent(const parameter::ParameterGroup& param,
@ -120,7 +125,7 @@ namespace Opm
const std::vector<double>& threshold_pressures_by_face,
const bool solvent);
std::unique_ptr<Solver> createSolver(const Wells* wells);
std::unique_ptr<Solver> createSolver(const WellModel& well_model);
void handleAdditionalWellInflow(SimulatorTimer& timer,
WellsManager& wells_manager,

4
opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp
View File

@ -63,7 +63,7 @@ namespace Opm
template <class GridT>
auto SimulatorFullyImplicitBlackoilSolvent<GridT>::
createSolver(const Wells* wells)
createSolver(const WellModel& well_model)
-> std::unique_ptr<Solver>
{
typedef typename Traits::Model Model;
@ -75,7 +75,7 @@ namespace Opm
BaseType::geo_,
BaseType::rock_comp_props_,
solvent_props_,
wells,
well_model,
BaseType::solver_,
BaseType::eclipse_state_,
BaseType::has_disgas_,

106
opm/autodiff/StandardWells.hpp
View File

@ -60,7 +60,18 @@ namespace Opm {
Eigen::Dynamic,
Eigen::RowMajor>;
// --------- Public methods ---------
explicit StandardWells(const Wells* wells);
StandardWells(const Wells* wells_arg);
void init(const BlackoilPropsAdInterface* fluid_arg,
const std::vector<bool>* active_arg,
const std::vector<PhasePresence>* pc_arg,
const VFPProperties* vfp_properties_arg,
const double gravity_arg,
const Vector& depth_arg);
const WellOps& wellOps() const;
int numPhases() const { return num_phases_; };
const Wells& wells() const;
@ -70,7 +81,7 @@ namespace Opm {
/// return true if wells are available on this process
bool localWellsActive() const;
const WellOps& wellOps() const;
int numWellVars() const;
/// Density of each well perforation
Vector& wellPerforationDensities(); // mutable version kept for BlackoilMultisegmentModel
@ -80,40 +91,14 @@ namespace Opm {
Vector& wellPerforationPressureDiffs(); // mutable version kept for BlackoilMultisegmentModel
const Vector& wellPerforationPressureDiffs() const;
template <class SolutionState, class WellState>
void computePropertiesForWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
const std::vector<PhasePresence>& pc,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf);
template <class WellState>
void computeWellConnectionDensitesPressures(const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf,
const std::vector<double>& depth_perf,
const double grav);
template <class ReservoirResidualQuant, class SolutionState>
void extractWellPerfProperties(const SolutionState& state,
const std::vector<ReservoirResidualQuant>& rq,
const int np,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const;
template <class SolutionState>
void computeWellFlux(const SolutionState& state,
const Opm::PhaseUsage& phase_usage,
const std::vector<bool>& active,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
Vector& aliveWells,
@ -126,19 +111,11 @@ namespace Opm {
template <class WellState>
void updateWellState(const Vector& dwells,
const double gravity,
const double dpmaxrel,
const Opm::PhaseUsage& pu,
const std::vector<bool>& active,
const VFPProperties& vfp_properties,
WellState& well_state);
template <class WellState>
void updateWellControls(const Opm::PhaseUsage& pu,
const double gravity,
const VFPProperties& vfp_properties,
const bool terminal_output,
const std::vector<bool>& active,
void updateWellControls(const bool terminal_output,
WellState& xw) const;
// TODO: should LinearisedBlackoilResidual also be a template class?
@ -152,39 +129,30 @@ namespace Opm {
void addWellControlEq(const SolutionState& state,
const WellState& xw,
const Vector& aliveWells,
const std::vector<bool> active,
const VFPProperties& vfp_properties,
const double gravity,
LinearisedBlackoilResidual& residual);
template <class SolutionState, class WellState>
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
const std::vector<PhasePresence>& phaseCondition,
const Vector& depth,
const double gravity);
const WellState& xw);
// state0 is non-constant, while it will not be used outside of the function
template <class SolutionState, class WellState>
void
computeWellPotentials(SolutionState& state0,
const std::vector<ADB>& mob_perfcells,
computeWellPotentials(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
const Opm::PhaseUsage& pu,
const std::vector<bool> active,
const VFPProperties& vfp_properties,
const bool compute_well_potentials,
const bool gravity,
SolutionState& state0,
WellState& well_state);
template <class SolutionState>
void
variableStateExtractWellsVars(const std::vector<int>& indices,
std::vector<ADB>& vars,
SolutionState& state) const;
void
variableStateWellIndices(std::vector<int>& indices,
int& next) const;
std::vector<int>
variableWellStateIndices() const;
@ -197,8 +165,38 @@ namespace Opm {
bool wells_active_;
const Wells* wells_;
const WellOps wops_;
const int num_phases_;
const BlackoilPropsAdInterface* fluid_;
const std::vector<bool>* active_;
const std::vector<PhasePresence>* phase_condition_;
const VFPProperties* vfp_properties_;
double gravity_;
// the depth of the all the cell centers
// for standard Wells, it the same with the perforation depth
Vector perf_cell_depth_;
Vector well_perforation_densities_;
Vector well_perforation_pressure_diffs_;
// protected methods
template <class SolutionState, class WellState>
void computePropertiesForWellConnectionPressures(const SolutionState& state,
const WellState& xw,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf);
template <class WellState>
void computeWellConnectionDensitesPressures(const WellState& xw,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf,
const std::vector<double>& depth_perf,
const double grav);
};

17
opm/autodiff/StandardWellsSolvent.hpp
View File

@ -37,7 +37,7 @@ namespace Opm {
using Base::computeWellConnectionDensitesPressures;
// --------- Public methods ---------
explicit StandardWellsSolvent(const Wells* wells);
StandardWellsSolvent(const Wells* wells_arg);
// added the Solvent related
void initSolvent(const SolventPropsAdFromDeck* solvent_props,
@ -47,9 +47,6 @@ namespace Opm {
template <class SolutionState, class WellState>
void computePropertiesForWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
const std::vector<PhasePresence>& pc,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
@ -59,26 +56,20 @@ namespace Opm {
template <class ReservoirResidualQuant, class SolutionState>
void extractWellPerfProperties(const SolutionState& state,
const std::vector<ReservoirResidualQuant>& rq,
const int np,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const;
template <class SolutionState, class WellState>
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
const std::vector<PhasePresence>& phaseCondition,
const Vector& depth,
const double gravity);
const WellState& xw);
protected:
const SolventPropsAdFromDeck* solvent_props_;
int solvent_pos_;
bool has_solvent_;
using Base::phase_condition_;
};

56
opm/autodiff/StandardWellsSolvent_impl.hpp
View File

@ -63,9 +63,6 @@ namespace Opm
StandardWellsSolvent::
computePropertiesForWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
const std::vector<PhasePresence>& pc,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
@ -99,46 +96,46 @@ namespace Opm
const ADB avg_press_ad = ADB::constant(avg_press);
std::vector<PhasePresence> perf_cond(nperf);
for (int perf = 0; perf < nperf; ++perf) {
perf_cond[perf] = pc[well_cells[perf]];
perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
}
const PhaseUsage& pu = fluid.phaseUsage();
const PhaseUsage& pu = fluid_->phaseUsage();
DataBlock b(nperf, pu.num_phases);
const Vector bw = fluid.bWat(avg_press_ad, perf_temp, well_cells).value();
const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
if (pu.phase_used[BlackoilPhases::Aqua]) {
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
}
assert(active[Oil]);
assert(active[Gas]);
assert((*active_)[Oil]);
assert((*active_)[Gas]);
const ADB perf_rv = subset(state.rv, well_cells);
const ADB perf_rs = subset(state.rs, well_cells);
const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const Vector bo = fluid.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
//const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
// const Vector rssat = fluidRsSat(avg_press, perf_so, well_cells);
const Vector rssat = fluid.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
} else {
rsmax_perf.assign(0.0, nperf);
}
V surf_dens_copy = superset(fluid.surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
V surf_dens_copy = superset(fluid_->surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
for (int phase = 1; phase < pu.num_phases; ++phase) {
if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) {
continue; // the gas surface density is added after the solvent is accounted for.
}
surf_dens_copy += superset(fluid.surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
surf_dens_copy += superset(fluid_->surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
// Unclear wether the effective or the pure values should be used for the wells
// the current usage of unmodified properties values gives best match.
//V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value();
Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
Vector rhog = fluid.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
Vector rhog = fluid_->surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
// to handle solvent related
if (has_solvent_) {
@ -150,7 +147,7 @@ namespace Opm
const ADB zero = ADB::constant(Vector::Zero(nc));
const ADB& ss = state.solvent_saturation;
const ADB& sg = (active[ Gas ]
const ADB& sg = ((*active_)[ Gas ]
? state.saturation[ pu.phase_pos[ Gas ] ]
: zero);
@ -181,7 +178,7 @@ namespace Opm
surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases);
// const Vector rvsat = fluidRvSat(avg_press, perf_so, well_cells);
const Vector rvsat = fluid.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
} else {
rvmax_perf.assign(0.0, nperf);
@ -201,12 +198,7 @@ namespace Opm
void
StandardWellsSolvent::
computeWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
const std::vector<PhasePresence>& phaseCondition,
const Vector& depth,
const double gravity)
const WellState& xw)
{
if( ! localWellsActive() ) return ;
// 1. Compute properties required by computeConnectionPressureDelta().
@ -216,16 +208,13 @@ namespace Opm
std::vector<double> rsmax_perf;
std::vector<double> rvmax_perf;
std::vector<double> surf_dens_perf;
computePropertiesForWellConnectionPressures(state, xw, fluid, active, phaseCondition,
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// Extract well connection depths
// TODO: depth_perf should be a member of the StandardWells class
const Vector pdepth = subset(depth, wellOps().well_cells);
const Vector pdepth = perf_cell_depth_;
const int nperf = wells().well_connpos[wells().number_of_wells];
const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
computeWellConnectionDensitesPressures(xw, fluid, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity);
computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_);
}
@ -238,16 +227,14 @@ namespace Opm
StandardWellsSolvent::
extractWellPerfProperties(const SolutionState& state,
const std::vector<ReservoirResidualQuant>& rq,
const int np,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const
{
Base::extractWellPerfProperties(state, rq, np, fluid, active, mob_perfcells, b_perfcells);
Base::extractWellPerfProperties(state, rq, mob_perfcells, b_perfcells);
// handle the solvent related
if (has_solvent_) {
int gas_pos = fluid.phaseUsage().phase_pos[Gas];
const Opm::PhaseUsage& pu = fluid_->phaseUsage();
int gas_pos = pu.phase_pos[Gas];
const std::vector<int>& well_cells = wellOps().well_cells;
const int nperf = well_cells.size();
// Gas and solvent is combinded and solved together
@ -260,10 +247,9 @@ namespace Opm
// A weighted sum of the b-factors of gas and solvent are used.
const int nc = rq[solvent_pos_].mob.size();
const Opm::PhaseUsage& pu = fluid.phaseUsage();
const ADB zero = ADB::constant(Vector::Zero(nc));
const ADB& ss = state.solvent_saturation;
const ADB& sg = (active[ Gas ]
const ADB& sg = ((*active_)[ Gas ]
? state.saturation[ pu.phase_pos[ Gas ] ]
: zero);

369
opm/autodiff/StandardWells_impl.hpp
View File

@ -74,6 +74,11 @@ namespace Opm
StandardWells::StandardWells(const Wells* wells_arg)
: wells_(wells_arg)
, wops_(wells_arg)
, num_phases_(wells_arg->number_of_phases)
, fluid_(nullptr)
, active_(nullptr)
, phase_condition_(nullptr)
, vfp_properties_(nullptr)
, well_perforation_densities_(Vector())
, well_perforation_pressure_diffs_(Vector())
{
@ -83,6 +88,26 @@ namespace Opm
void
StandardWells::init(const BlackoilPropsAdInterface* fluid_arg,
const std::vector<bool>* active_arg,
const std::vector<PhasePresence>* pc_arg,
const VFPProperties* vfp_properties_arg,
const double gravity_arg,
const Vector& depth_arg)
{
fluid_ = fluid_arg;
active_ = active_arg;
phase_condition_ = pc_arg;
vfp_properties_ = vfp_properties_arg;
gravity_ = gravity_arg;
perf_cell_depth_ = subset(depth_arg, wellOps().well_cells);;
}
const Wells& StandardWells::wells() const
{
assert(wells_ != 0);
@ -120,6 +145,18 @@ namespace Opm
int
StandardWells::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;
}
const StandardWells::WellOps&
StandardWells::wellOps() const
{
@ -173,16 +210,11 @@ namespace Opm
StandardWells::
computePropertiesForWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
const std::vector<PhasePresence>& pc,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf)
{
static_cast<void>(active); // Silence unused argument warning in release mode.
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
@ -209,31 +241,31 @@ namespace Opm
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]];
perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
}
const PhaseUsage& pu = fluid.phaseUsage();
const PhaseUsage& pu = fluid_->phaseUsage();
DataBlock b(nperf, pu.num_phases);
if (pu.phase_used[BlackoilPhases::Aqua]) {
const Vector bw = fluid.bWat(avg_press_ad, perf_temp, well_cells).value();
const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
}
assert(active[Oil]);
assert((*active_)[Oil]);
const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const ADB perf_rs = (state.rs.size() > 0) ? subset(state.rs, well_cells) : ADB::null();
const Vector bo = fluid.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = (state.rv.size() > 0) ? subset(state.rv, well_cells) : ADB::null();
const Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
const Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
}
if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
const Vector rssat = fluid.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
const Vector rvsat = fluid.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
}
@ -243,9 +275,9 @@ namespace Opm
// Surface density.
// The compute density segment wants the surface densities as
// an np * number of wells cells array
Vector rho = superset(fluid.surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
Vector rho = superset(fluid_->surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
for (int phase = 1; phase < pu.num_phases; ++phase) {
rho += superset(fluid.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
rho += superset(fluid_->surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
}
surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases);
@ -259,7 +291,6 @@ namespace Opm
void
StandardWells::
computeWellConnectionDensitesPressures(const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
@ -270,7 +301,7 @@ namespace Opm
// Compute densities
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid.phaseUsage(),
wells(), xw, fluid_->phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
const int nperf = wells().well_connpos[wells().number_of_wells];
@ -293,12 +324,7 @@ namespace Opm
void
StandardWells::
computeWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
const std::vector<PhasePresence>& phaseCondition,
const Vector& depth,
const double gravity)
const WellState& xw)
{
if( ! localWellsActive() ) return ;
// 1. Compute properties required by computeConnectionPressureDelta().
@ -308,16 +334,13 @@ namespace Opm
std::vector<double> rsmax_perf;
std::vector<double> rvmax_perf;
std::vector<double> surf_dens_perf;
computePropertiesForWellConnectionPressures(state, xw, fluid, active, phaseCondition,
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// Extract well connection depths
// TODO: depth_perf should be a member of the StandardWells class
const Vector pdepth = subset(depth, wellOps().well_cells);
const Vector& pdepth = perf_cell_depth_;
const int nperf = wells().well_connpos[wells().number_of_wells];
const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
computeWellConnectionDensitesPressures(xw, fluid, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity);
computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_);
}
@ -330,9 +353,6 @@ namespace Opm
StandardWells::
extractWellPerfProperties(const SolutionState& /* state */,
const std::vector<ReservoirResidualQuant>& rq,
const int np,
const BlackoilPropsAdInterface& /* fluid */,
const std::vector<bool>& /* active */,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const
{
@ -344,9 +364,9 @@ namespace Opm
return;
} else {
const std::vector<int>& well_cells = wellOps().well_cells;
mob_perfcells.resize(np, ADB::null());
b_perfcells.resize(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
mob_perfcells.resize(num_phases_, ADB::null());
b_perfcells.resize(num_phases_, ADB::null());
for (int phase = 0; phase < num_phases_; ++phase) {
mob_perfcells[phase] = subset(rq[phase].mob, well_cells);
b_perfcells[phase] = subset(rq[phase].b, well_cells);
}
@ -362,8 +382,6 @@ namespace Opm
void
StandardWells::
computeWellFlux(const SolutionState& state,
const Opm::PhaseUsage& pu,
const std::vector<bool>& active,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
Vector& aliveWells,
@ -371,7 +389,7 @@ namespace Opm
{
if( ! localWellsActive() ) return ;
const int np = wells().number_of_phases;
const int np = num_phases_;
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
Vector Tw = Eigen::Map<const Vector>(wells().WI, nperf);
@ -428,7 +446,8 @@ namespace Opm
const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
cq_ps[phase] = b_perfcells[phase] * cq_p;
}
if (active[Oil] && active[Gas]) {
const Opm::PhaseUsage& pu = fluid_->phaseUsage();
if ((*active_)[Oil] && (*active_)[Gas]) {
const int oilpos = pu.phase_pos[Oil];
const int gaspos = pu.phase_pos[Gas];
const ADB cq_psOil = cq_ps[oilpos];
@ -475,12 +494,12 @@ namespace Opm
// compute volume ratio between connection at standard conditions
ADB volumeRatio = ADB::constant(Vector::Zero(nperf));
if (active[Water]) {
if ((*active_)[Water]) {
const int watpos = pu.phase_pos[Water];
volumeRatio += cmix_s[watpos] / b_perfcells[watpos];
}
if (active[Oil] && active[Gas]) {
if ((*active_)[Oil] && (*active_)[Gas]) {
// Incorporate RS/RV factors if both oil and gas active
const ADB& rv_perfcells = subset(state.rv, well_cells);
const ADB& rs_perfcells = subset(state.rs, well_cells);
@ -496,11 +515,11 @@ namespace Opm
volumeRatio += tmp_gas / b_perfcells[gaspos];
}
else {
if (active[Oil]) {
if ((*active_)[Oil]) {
const int oilpos = pu.phase_pos[Oil];
volumeRatio += cmix_s[oilpos] / b_perfcells[oilpos];
}
if (active[Gas]) {
if ((*active_)[Gas]) {
const int gaspos = pu.phase_pos[Gas];
volumeRatio += cmix_s[gaspos] / b_perfcells[gaspos];
}
@ -567,17 +586,12 @@ namespace Opm
void
StandardWells::
updateWellState(const Vector& dwells,
const double gravity,
const double dpmaxrel,
const Opm::PhaseUsage& pu,
const std::vector<bool>& active,
const VFPProperties& vfp_properties,
WellState& well_state)
{
if( localWellsActive() )
{
// TODO: these parameter should be stored in the StandardWells class
const int np = wells().number_of_phases;
const int np = num_phases_;
const int nw = wells().number_of_wells;
// Extract parts of dwells corresponding to each part.
@ -605,6 +619,8 @@ namespace Opm
const Vector bhp = bhp_old - dbhp_limited;
std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
const Opm::PhaseUsage& pu = fluid_->phaseUsage();
//Loop over all wells
#pragma omp parallel for schedule(static)
for (int w = 0; w < nw; ++w) {
@ -619,13 +635,13 @@ namespace Opm
double liquid = 0.0;
double vapour = 0.0;
if (active[ Water ]) {
if ((*active_)[ Water ]) {
aqua = wr[w*np + pu.phase_pos[ Water ] ];
}
if (active[ Oil ]) {
if ((*active_)[ Oil ]) {
liquid = wr[w*np + pu.phase_pos[ Oil ] ];
}
if (active[ Gas ]) {
if ((*active_)[ Gas ]) {
vapour = wr[w*np + pu.phase_pos[ Gas ] ];
}
@ -635,17 +651,17 @@ namespace Opm
const WellType& well_type = wells().type[w];
if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties.getInj()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity);
wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity_);
well_state.thp()[w] = vfp_properties.getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp);
well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp);
}
else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties.getProd()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity);
wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity_);
well_state.thp()[w] = vfp_properties.getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq);
well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq);
}
else {
OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
@ -666,11 +682,7 @@ namespace Opm
template <class WellState>
void
StandardWells::
updateWellControls(const Opm::PhaseUsage& pu,
const double gravity,
const VFPProperties& vfp_properties,
const bool terminal_output,
const std::vector<bool>& active,
updateWellControls(const bool terminal_output,
WellState& xw) const
{
if( !localWellsActive() ) return ;
@ -732,13 +744,15 @@ namespace Opm
double liquid = 0.0;
double vapour = 0.0;
if (active[ Water ]) {
const Opm::PhaseUsage& pu = fluid_->phaseUsage();
if ((*active_)[ Water ]) {
aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if (active[ Oil ]) {
if ((*active_)[ Oil ]) {
liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if (active[ Gas ]) {
if ((*active_)[ Gas ]) {
vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ];
}
@ -751,17 +765,17 @@ namespace Opm
if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties.getInj()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity);
wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity_);
xw.bhp()[w] = vfp_properties.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
}
else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties.getProd()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity);
wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity_);
xw.bhp()[w] = vfp_properties.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
}
else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
@ -850,12 +864,9 @@ namespace Opm
template <class SolutionState, class WellState>
void
StandardWells::addWellControlEq(const SolutionState& state,
const WellState& xw,
const Vector& aliveWells,
const std::vector<bool> active,
const VFPProperties& vfp_properties,
const double gravity,
LinearisedBlackoilResidual& residual)
const WellState& xw,
const Vector& aliveWells,
LinearisedBlackoilResidual& residual)
{
if( ! localWellsActive() ) return;
@ -866,13 +877,13 @@ namespace Opm
ADB liquid = ADB::constant(Vector::Zero(nw));
ADB vapour = ADB::constant(Vector::Zero(nw));
if (active[Water]) {
if ((*active_)[Water]) {
aqua += subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
}
if (active[Oil]) {
if ((*active_)[Oil]) {
liquid += subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
}
if (active[Gas]) {
if ((*active_)[Gas]) {
vapour += subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
}
@ -931,7 +942,7 @@ namespace Opm
thp_inj_target_v[w] = target;
alq_v[w] = -1e100;
vfp_ref_depth_v[w] = vfp_properties.getInj()->getTable(table_id)->getDatumDepth();
vfp_ref_depth_v[w] = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
thp_inj_elems.push_back(w);
}
@ -940,7 +951,7 @@ namespace Opm
thp_prod_target_v[w] = target;
alq_v[w] = well_controls_iget_alq(wc, current);
vfp_ref_depth_v[w] = vfp_properties.getProd()->getTable(table_id)->getDatumDepth();
vfp_ref_depth_v[w] = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
thp_prod_elems.push_back(w);
}
@ -978,12 +989,11 @@ namespace Opm
const ADB thp_inj_target = ADB::constant(thp_inj_target_v);
const ADB thp_prod_target = ADB::constant(thp_prod_target_v);
const ADB alq = ADB::constant(alq_v);
const ADB bhp_from_thp_inj = vfp_properties.getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target);
const ADB bhp_from_thp_prod = vfp_properties.getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq);
const ADB bhp_from_thp_inj = vfp_properties_->getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target);
const ADB bhp_from_thp_prod = vfp_properties_->getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq);
//Perform hydrostatic correction to computed targets
// double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const Vector dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, wellPerforationDensities(), gravity);
const Vector dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, wellPerforationDensities(), gravity_);
const ADB dp = ADB::constant(dp_v);
const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw);
const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw);
@ -1020,102 +1030,95 @@ namespace Opm
template <class SolutionState, class WellState>
void
StandardWells::computeWellPotentials(SolutionState& state0,
const std::vector<ADB>& mob_perfcells,
StandardWells::computeWellPotentials(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
const Opm::PhaseUsage& pu,
const std::vector<bool> active,
const VFPProperties& vfp_properties,
const bool compute_well_potentials,
const bool gravity,
SolutionState& state0,
WellState& well_state)
{
//only compute well potentials if they are needed
if (compute_well_potentials) {
const int nw = wells().number_of_wells;
const int np = wells().number_of_phases;
Vector bhps = Vector::Zero(nw);
for (int w = 0; w < nw; ++w) {
const WellControls* ctrl = wells().ctrls[w];
const int nwc = well_controls_get_num(ctrl);
//Loop over all controls until we find a BHP control
//or a THP control that specifies what we need.
//Pick the value that gives the most restrictive flow
for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
const int nw = wells().number_of_wells;
const int np = wells().number_of_phases;
const Opm::PhaseUsage& pu = fluid_->phaseUsage();
if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
bhps[w] = well_controls_iget_target(ctrl, ctrl_index);
}
Vector bhps = Vector::Zero(nw);
for (int w = 0; w < nw; ++w) {
const WellControls* ctrl = wells().ctrls[w];
const int nwc = well_controls_get_num(ctrl);
//Loop over all controls until we find a BHP control
//or a THP control that specifies what we need.
//Pick the value that gives the most restrictive flow
for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
if (active[ Water ]) {
aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if (active[ Oil ]) {
liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if (active[ Gas ]) {
vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
}
const int vfp = well_controls_iget_vfp(ctrl, ctrl_index);
const double& thp = well_controls_iget_target(ctrl, ctrl_index);
const double& alq = well_controls_iget_alq(ctrl, ctrl_index);
//Set *BHP* target by calculating bhp from THP
const WellType& well_type = wells().type[w];
if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties.getInj()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity);
const double bhp = vfp_properties.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
// apply the strictest of the bhp controlls i.e. smallest bhp for injectors
if ( bhp < bhps[w]) {
bhps[w] = bhp;
}
}
else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties.getProd()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity);
const double bhp = vfp_properties.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
// apply the strictest of the bhp controlls i.e. largest bhp for producers
if ( bhp > bhps[w]) {
bhps[w] = bhp;
}
}
else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
}
if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
bhps[w] = well_controls_iget_target(ctrl, ctrl_index);
}
if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
if ((*active_)[ Water ]) {
aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if ((*active_)[ Oil ]) {
liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if ((*active_)[ Gas ]) {
vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
}
const int vfp = well_controls_iget_vfp(ctrl, ctrl_index);
const double& thp = well_controls_iget_target(ctrl, ctrl_index);
const double& alq = well_controls_iget_alq(ctrl, ctrl_index);
//Set *BHP* target by calculating bhp from THP
const WellType& well_type = wells().type[w];
if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity_);
const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
// apply the strictest of the bhp controlls i.e. smallest bhp for injectors
if ( bhp < bhps[w]) {
bhps[w] = bhp;
}
}
else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity_);
const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
// apply the strictest of the bhp controlls i.e. largest bhp for producers
if ( bhp > bhps[w]) {
bhps[w] = bhp;
}
}
else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
}
}
// use bhp limit from control
state0.bhp = ADB::constant(bhps);
// compute well potentials
Vector aliveWells;
std::vector<ADB> well_potentials;
computeWellFlux(state0, pu, active, mob_perfcells, b_perfcells, aliveWells, well_potentials);
// store well potentials in the well state
// transform to a single vector instead of separate vectors pr phase
const int nperf = wells().well_connpos[nw];
V cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np);
}
well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np);
}
// use bhp limit from control
state0.bhp = ADB::constant(bhps);
// compute well potentials
Vector aliveWells;
std::vector<ADB> well_potentials;
computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, well_potentials);
// store well potentials in the well state
// transform to a single vector instead of separate vectors pr phase
const int nperf = wells().well_connpos[nw];
Vector cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np);
}
well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np);
}
@ -1134,6 +1137,24 @@ namespace Opm
template <class SolutionState>
void
StandardWells::
variableStateExtractWellsVars(const std::vector<int>& indices,
std::vector<ADB>& vars,
SolutionState& state) const
{
// Qs.
state.qs = std::move(vars[indices[Qs]]);
// Bhp.
state.bhp = std::move(vars[indices[Bhp]]);
}
std::vector<int>
StandardWells::variableWellStateIndices() const
{
@ -1167,12 +1188,12 @@ namespace Opm
// The transpose() below switches the ordering.
const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
const V qs = Eigen::Map<const V>(wrates.data(), nw*np);
const Vector qs = Eigen::Map<const V>(wrates.data(), nw*np);
vars0.push_back(qs);
// Initial well bottom-hole pressure.
assert (not xw.bhp().empty());
const V bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
const Vector bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
vars0.push_back(bhp);
}
else

6
opm/autodiff/WellStateMultiSegment.hpp
View File

@ -25,7 +25,8 @@
#include <opm/core/well_controls.h>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/WellMultiSegment.hpp>
// #include <opm/autodiff/WellMultiSegment.hpp>
#include <opm/autodiff/MultisegmentWells.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <vector>
#include <cassert>
@ -64,8 +65,9 @@ namespace Opm
/// to give useful initial values to the bhp(), wellRates()
/// and perfPhaseRates() fields, depending on controls
template <class ReservoirState, class PrevWellState>
void init(const std::vector<WellMultiSegmentConstPtr>& wells, const ReservoirState& state, const PrevWellState& prevState)
void init(const MultisegmentWells& ms_wells, const ReservoirState& state, const PrevWellState& prevState)
{
const std::vector<WellMultiSegmentConstPtr>& wells = ms_wells.msWells();
const int nw = wells.size();
nseg_ = 0;
nperf_ = 0;

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

@ -79,7 +79,7 @@ namespace Opm {
const DerivedGeology& geo,
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const Wells* wells,
const StandardWells& well_model,
const NewtonIterationBlackoilInterface& linsolver,
EclipseStateConstPtr eclipse_state,
const bool has_disgas,
@ -179,7 +179,7 @@ namespace Opm {
// --------- Protected methods ---------
// Need to declare Base members we want to use here.
using Base::stdWells;
using Base::wellModel;
using Base::wells;
using Base::wellsActive;
using Base::variableState;

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

@ -80,7 +80,7 @@ namespace Opm {
const DerivedGeology& geo,
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const Wells* wells_arg,
const StandardWells& well_model,
const NewtonIterationBlackoilInterface& linsolver,
EclipseStateConstPtr eclipse_state,
const bool has_disgas,
@ -92,7 +92,7 @@ namespace Opm {
const std::vector<double>& wells_perf_length,
const std::vector<double>& wells_bore_diameter,
const bool terminal_output)
: Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver, eclipse_state,
: Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver, eclipse_state,
has_disgas, has_vapoil, terminal_output),
polymer_props_ad_(polymer_props_ad),
has_polymer_(has_polymer),
@ -245,10 +245,10 @@ namespace Opm {
const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], AutoDiffGrid::numCells(grid_));
const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
// Compute polymer accumulation term.
rq_[poly_pos_].accum[aix] = pv_mult * rq_[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol)
rq_[poly_pos_].accum[aix] = pv_mult * rq_[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol)
+ pv_mult * rho_rock * (1. - phi) / phi * ads;
}
}
@ -498,10 +498,8 @@ namespace Opm {
// Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const V depth = cellCentroidsZToEigen(grid_);
// updateWellControls(well_state);
stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state);
wellModel().updateWellControls(terminal_output_, well_state);
// Create the primary variables.
SolutionState state = variableState(reservoir_state, well_state);
@ -514,7 +512,7 @@ namespace Opm {
// and well connection pressures.
computeAccum(state0, 0);
// computeWellConnectionPressures(state0, well_state);
stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity);
wellModel().computeWellConnectionPressures(state0, well_state);
}
// OPM_AD_DISKVAL(state.pressure);
@ -534,27 +532,18 @@ namespace Opm {
return;
}
V aliveWells;
const int np = wells().number_of_phases;
std::vector<ADB> cq_s(np, ADB::null());
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
std::vector<ADB> mob_perfcells(np, ADB::null());
std::vector<ADB> b_perfcells(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
mob_perfcells[phase] = subset(rq_[phase].mob, well_cells);
b_perfcells[phase] = subset(rq_[phase].b, well_cells);
}
std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells;
wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
}
stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s);
V aliveWells;
std::vector<ADB> cq_s;
wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
if (has_plyshlog_) {
std::vector<double> water_vel_wells;
@ -573,11 +562,11 @@ namespace Opm {
mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb;
}
stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s);
stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
stdWells().addWellFluxEq(cq_s, state, residual_);
wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
wellModel().addWellFluxEq(cq_s, state, residual_);
addWellContributionToMassBalanceEq(cq_s, state, well_state);
stdWells().addWellControlEq(state, well_state, aliveWells, active_, vfp_properties_, gravity, residual_);
wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
}
@ -736,8 +725,8 @@ namespace Opm {
ADB b_perfcells = subset(rq_[water_pos].b, well_cells);
const ADB& p_perfcells = subset(state.pressure, well_cells);
const V& cdp = stdWells().wellPerforationPressureDiffs();
const ADB perfpressure = (stdWells().wellOps().w2p * state.bhp) + cdp;
const V& cdp = wellModel().wellPerforationPressureDiffs();
const ADB perfpressure = (wellModel().wellOps().w2p * state.bhp) + cdp;
// Pressure drawdown (also used to determine direction of flow)
const ADB drawdown = p_perfcells - perfpressure;

6
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp
View File

@ -81,6 +81,8 @@ namespace Opm
template <class GridT>
class SimulatorFullyImplicitBlackoilPolymer;
class StandardWells;
template<class GridT>
struct SimulatorTraits<SimulatorFullyImplicitBlackoilPolymer<GridT> >
{
@ -90,6 +92,7 @@ namespace Opm
typedef GridT Grid;
typedef BlackoilPolymerModel<Grid> Model;
typedef NonlinearSolver<Model> Solver;
typedef StandardWells WellModel;
};
/// Class collecting all necessary components for a blackoil simulation with polymer
@ -103,6 +106,7 @@ namespace Opm
typedef SimulatorTraits<ThisType> Traits;
typedef typename Traits::Solver Solver;
typedef typename Traits::WellModel WellModel;
public:
SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param,
@ -123,7 +127,7 @@ namespace Opm
Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face);
std::unique_ptr<Solver> createSolver(const Wells* wells);
std::unique_ptr<Solver> createSolver(const WellModel& well_model);
void handleAdditionalWellInflow(SimulatorTimer& timer,

4
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp
View File

@ -62,7 +62,7 @@ namespace Opm
template <class GridT>
auto SimulatorFullyImplicitBlackoilPolymer<GridT>::
createSolver(const Wells* wells)
createSolver(const WellModel& well_model)
-> std::unique_ptr<Solver>
{
typedef typename Traits::Model Model;
@ -74,7 +74,7 @@ namespace Opm
BaseType::geo_,
BaseType::rock_comp_props_,
polymer_props_,
wells,
well_model,
BaseType::solver_,
BaseType::eclipse_state_,
BaseType::has_disgas_,

6
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp
View File

@ -71,6 +71,8 @@ namespace Opm
template <class GridT>
class SimulatorFullyImplicitCompressiblePolymer;
class StandardWells;
template <class GridT>
struct SimulatorTraits<SimulatorFullyImplicitCompressiblePolymer<GridT> >
{
@ -79,6 +81,7 @@ namespace Opm
typedef BlackoilOutputWriter OutputWriter;
typedef GridT Grid;
typedef FullyImplicitCompressiblePolymerSolver Solver;
typedef StandardWells WellModel;
/// Dummy class, this Solver does not use a Model.
struct Model
{
@ -94,6 +97,7 @@ namespace Opm
typedef SimulatorFullyImplicitCompressiblePolymer ThisType;
typedef SimulatorBase<ThisType> BaseType;
typedef typename BaseType::Solver Solver;
typedef typename BaseType::WellModel WellModel;
public:
/// Initialise from parameters and objects to observe.
@ -109,7 +113,7 @@ namespace Opm
NewtonIterationBlackoilInterface& linsolver,
const double* gravity);
std::unique_ptr<Solver> createSolver(const Wells* wells);
std::unique_ptr<Solver> createSolver(const WellModel& well_model);
void handleAdditionalWellInflow(SimulatorTimer& timer,
WellsManager& wells_manager,

6
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer_impl.hpp
View File

@ -57,7 +57,7 @@ SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param
template <class GridT>
auto SimulatorFullyImplicitCompressiblePolymer<GridT>::
createSolver(const Wells* wells)
createSolver(const WellModel& well_model)
-> std::unique_ptr<Solver>
{
return std::unique_ptr<Solver>(new Solver(BaseType::grid_,
@ -65,7 +65,9 @@ createSolver(const Wells* wells)
BaseType::geo_,
BaseType::rock_comp_props_,
polymer_props_,
*wells,
// *wells,
// TODO: it is resulted from refactoring of other simulators.
well_model.wells(),
BaseType::solver_));
}

28
tests/test_multisegmentwells.cpp
View File

@ -106,31 +106,7 @@ struct SetupMSW {
const Wells* wells = wells_manager.c_wells();
const std::vector<Opm::WellConstPtr>& wells_ecl = ecl_state->getSchedule()->getWells(current_timestep);
std::vector<Opm::WellMultiSegmentConstPtr> wells_multisegment;
wells_multisegment.reserve(wells_ecl.size());
for (size_t i = 0; i < wells_ecl.size(); ++i) {
// not processing SHUT wells.
if (wells_ecl[i]->getStatus(current_timestep) == Opm::WellCommon::SHUT) {
continue;
}
// checking if the well can be found in the wells
const std::string& well_name = wells_ecl[i]->name();
// number of wells in wells
const int nw_wells = wells->number_of_wells;
int index_well;
for (index_well = 0; index_well < nw_wells; ++index_well) {
if (well_name == std::string(wells->name[index_well])) {
break;
}
}
// have to be able to be found
assert(index_well != nw_wells);
wells_multisegment.push_back(std::make_shared<Opm::WellMultiSegment>(wells_ecl[i], current_timestep, wells));
}
ms_wells.reset(new Opm::MultisegmentWells(wells_multisegment, wells->number_of_phases));
ms_wells.reset(new Opm::MultisegmentWells(wells, wells_ecl, current_timestep));
};
std::shared_ptr<const Opm::MultisegmentWells> ms_wells;
@ -175,7 +151,7 @@ BOOST_AUTO_TEST_CASE(testStructure)
BOOST_CHECK_EQUAL(nperf, ms_wells->numPerf());
BOOST_CHECK_EQUAL(nw, ms_wells->topWellSegments().size());
BOOST_CHECK_EQUAL(nw, ms_wells->wells().size());
BOOST_CHECK_EQUAL(nw, ms_wells->msWells().size());
BOOST_CHECK_EQUAL(0, ms_wells->topWellSegments()[0]);
BOOST_CHECK_EQUAL(1, ms_wells->topWellSegments()[1]);
}