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

Remove unused code and remove Eigen vectors

Browse Source 
-- isRS and phaseCondition is removed and hydroCarbonState in the state
is used instead
-- input of pressurediffs to computeHydrostaticCorrection() is changed
to double from Vector in WellHelpers.hpp
This commit is contained in:
Tor Harald Sandve
2016-09-07 11:21:51 +02:00
parent 83ff3271af
commit a4dcc4b13d
5 changed files with 84 additions and 780 deletions
Download Patch File Download Diff File Expand all files Collapse all files

416
opm/autodiff/BlackoilModelEbos.hpp
View File

@@ -153,21 +153,15 @@ namespace Opm {
, grid_(ebosSimulator_.gridManager().grid()) , grid_(ebosSimulator_.gridManager().grid())
, fluid_ (fluid) , fluid_ (fluid)
, geo_ (geo) , geo_ (geo)
, rock_comp_props_(rock_comp_props)
, vfp_properties_( , vfp_properties_(
eclState().getTableManager().getVFPInjTables(), eclState().getTableManager().getVFPInjTables(),
eclState().getTableManager().getVFPProdTables()) eclState().getTableManager().getVFPProdTables())
, linsolver_ (linsolver) , linsolver_ (linsolver)
, active_(detail::activePhases(fluid.phaseUsage())) , active_(detail::activePhases(fluid.phaseUsage()))
, cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid_)))
, has_disgas_(FluidSystem::enableDissolvedGas()) , has_disgas_(FluidSystem::enableDissolvedGas())
, has_vapoil_(FluidSystem::enableVaporizedOil()) , has_vapoil_(FluidSystem::enableVaporizedOil())
, param_( param ) , param_( param )
, phaseCondition_(AutoDiffGrid::numCells(grid_))
, well_model_ (well_model) , well_model_ (well_model)
, isRs_(V::Zero(AutoDiffGrid::numCells(grid_)))
, isRv_(V::Zero(AutoDiffGrid::numCells(grid_)))
, isSg_(V::Zero(AutoDiffGrid::numCells(grid_)))
, terminal_output_ (terminal_output) , terminal_output_ (terminal_output)
, current_relaxation_(1.0) , current_relaxation_(1.0)
, isBeginReportStep_(false) , isBeginReportStep_(false)
@@ -176,7 +170,7 @@ namespace Opm {
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const std::vector<double> pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); const std::vector<double> pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
const std::vector<double> depth(geo_.z().data(), geo_.z().data() + geo_.z().size()); const std::vector<double> depth(geo_.z().data(), geo_.z().data() + geo_.z().size());
well_model_.init(&fluid_, &active_, &phaseCondition_, &vfp_properties_, gravity, depth, pv); well_model_.init(&fluid_, &active_, &vfp_properties_, gravity, depth, pv);
wellModel().setWellsActive( localWellsActive() ); wellModel().setWellsActive( localWellsActive() );
global_nc_ = Opm::AutoDiffGrid::numCells(grid_); global_nc_ = Opm::AutoDiffGrid::numCells(grid_);
} }
@@ -188,13 +182,10 @@ namespace Opm {
/// \param[in] timer simulation timer /// \param[in] timer simulation timer
/// \param[in, out] reservoir_state reservoir state variables /// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables /// \param[in, out] well_state well state variables
void prepareStep(const SimulatorTimerInterface& timer, void prepareStep(const SimulatorTimerInterface& /*timer*/,
const ReservoirState& reservoir_state, const ReservoirState& /*reservoir_state*/,
const WellState& /* well_state */) const WellState& /* well_state */)
{ {
if (active_[Gas]) {
updatePrimalVariableFromState(reservoir_state);
}
} }
@@ -219,7 +210,7 @@ namespace Opm {
// the mass balance for each active phase, the well flux and the well equations. // the mass balance for each active phase, the well flux and the well equations.
residual_norms_history_.clear(); residual_norms_history_.clear();
current_relaxation_ = 1.0; current_relaxation_ = 1.0;
dx_old_ = 0.0; //V::Zero(sizeNonLinear()); dx_old_ = 0.0;
} }
IterationReport iter_report = assemble(timer, iteration, reservoir_state, well_state); IterationReport iter_report = assemble(timer, iteration, reservoir_state, well_state);
std::vector<double> residual_norms; std::vector<double> residual_norms;
@@ -236,7 +227,7 @@ namespace Opm {
const int nw = wellModel().wells().number_of_wells; const int nw = wellModel().wells().number_of_wells;
BVector x(nc); BVector x(nc);
BVector xw(nw); BVector xw(nw);
V dx = solveJacobianSystem(result, x, xw); solveJacobianSystem(result, x, xw);
// Stabilize the nonlinear update. // Stabilize the nonlinear update.
bool isOscillate = false; bool isOscillate = false;
@@ -255,24 +246,9 @@ namespace Opm {
// Apply the update, applying model-dependent // Apply the update, applying model-dependent
// limitations and chopping of the update. // limitations and chopping of the update.
//auto well_state2 = well_state;
//auto reservoir_state2 = reservoir_state;
//updateState(dx, reservoir_state2, well_state);
updateState(x,reservoir_state); updateState(x,reservoir_state);
wellModel().updateWellState(xw, well_state); wellModel().updateWellState(xw, well_state);
// for (int c = 0; c < nc; ++c) {
// printIf(c, reservoir_state.pressure()[c], reservoir_state2.pressure()[c], 1e5, "pressure");
// printIf(c, reservoir_state.rv()[c], reservoir_state2.rv()[c] , 1e-3, "rv");
// printIf(c, reservoir_state.gasoilratio()[c], reservoir_state2.gasoilratio()[c], 1, "rs");
// printIf(c, reservoir_state.saturation()[3*c+0], reservoir_state2.saturation()[3*c+0], 1e-3, "sw");
// printIf(c, reservoir_state.saturation()[3*c+1], reservoir_state2.saturation()[3*c+1], 1e-3, "so");
// printIf(c, reservoir_state.saturation()[3*c+2], reservoir_state2.saturation()[3*c+2], 1e-3, "sg");
// printIf(c, reservoir_state.hydroCarbonState()[c], reservoir_state2.hydroCarbonState()[c], 1e-6,"state");
// }
} }
const bool failed = false; // Not needed in this model. const bool failed = false; // Not needed in this model.
const int linear_iters = must_solve ? result.iterations : 0; const int linear_iters = must_solve ? result.iterations : 0;
@@ -322,13 +298,6 @@ namespace Opm {
OPM_THROW(Opm::NumericalProblem,"no convergence"); OPM_THROW(Opm::NumericalProblem,"no convergence");
} }
typedef double Scalar;
typedef Dune::FieldVector<Scalar, 3 > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, 3, 3 > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
typedef Dune::MatrixAdapter<Mat,BVector,BVector> Operator;
auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
auto& ebosResid = ebosSimulator_.model().linearizer().residual(); auto& ebosResid = ebosSimulator_.model().linearizer().residual();
@@ -396,7 +365,7 @@ namespace Opm {
/// Solve the Jacobian system Jx = r where J is the Jacobian and /// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual. /// r is the residual.
V solveJacobianSystem(Dune::InverseOperatorResult& result, BVector& x, BVector& xw) const void solveJacobianSystem(Dune::InverseOperatorResult& result, BVector& x, BVector& xw) const
{ {
typedef double Scalar; typedef double Scalar;
@@ -415,49 +384,17 @@ namespace Opm {
SeqPreconditioner precond(opA.getmat(), relax); SeqPreconditioner precond(opA.getmat(), relax);
Dune::SeqScalarProduct<BVector> sp; Dune::SeqScalarProduct<BVector> sp;
wellModel().apply(ebosJac, ebosResid); wellModel().apply(ebosJac, ebosResid);
Dune::BiCGSTABSolver<BVector> linsolve(opA, sp, precond, Dune::BiCGSTABSolver<BVector> linsolve(opA, sp, precond,
0.01, 0.01,
100, 100,
false); false);
const int np = numPhases();
const int nc = AutoDiffGrid::numCells(grid_);
// Solve system. // Solve system.
//BVector x(ebosJac.M());
x = 0.0; x = 0.0;
linsolve.apply(x, ebosResid, result); linsolve.apply(x, ebosResid, result);
const int nw = wellModel().wells().number_of_wells;
//BVector xw(nw);
xw = 0.0; xw = 0.0;
wellModel().recoverVariable(x, xw); wellModel().recoverVariable(x, xw);
// convert to V;
V dx( sizeNonLinear() );
for( int p=0; p<np; ++p) {
for (int i = 0; i < nc; ++i) {
int idx = i + nc*ebosCompToFlowPhaseIdx(p);
dx(idx) = x[i][p];
}
for (int w = 0; w < nw; ++w) {
int idx = w + nw*p + nc*np;
dx(idx) = xw[w][flowPhaseToEbosCompIdx(p)];
}
}
//V dx2 = linsolver_.computeNewtonIncrement(residual_);
//std::cout << "------dx------- " << std::endl;
//std::cout << dx << std::endl;
//std::cout << "------dx2------- " << std::endl;
//std::cout << dx2 << std::endl;
//return dx;
return dx;
} }
/// Apply an update to the primary variables, chopped if appropriate. /// Apply an update to the primary variables, chopped if appropriate.
@@ -559,6 +496,7 @@ namespace Opm {
so = 0; so = 0;
sg = 1.0 - sw - so; sg = 1.0 - sw - so;
double& rv = reservoir_state.rv()[cell_idx]; double& rv = reservoir_state.rv()[cell_idx];
// use gas pressure?
double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(0, reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(0, reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
rv = rvSat*(1-epsilon); rv = rvSat*(1-epsilon);
} }
@@ -595,251 +533,8 @@ namespace Opm {
} }
} }
//wellModel().updateWellState(dwells, well_state);
// Update phase conditions used for property calculations.
updatePhaseCondFromPrimalVariable(reservoir_state);
} }
/// Apply an update to the primary variables, chopped if appropriate.
/// \param[in] dx updates to apply to primary variables
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
void updateState(const V& dx,
ReservoirState& reservoir_state,
WellState& well_state)
{
using namespace Opm::AutoDiffGrid;
const int np = fluid_.numPhases();
const int nc = numCells(grid_);
const V null;
assert(null.size() == 0);
const V zero = V::Zero(nc);
// Extract parts of dx corresponding to each part.
const V dp = subset(dx, Span(nc));
int varstart = nc;
const V dsw = active_[Water] ? subset(dx, Span(nc, 1, varstart)) : null;
varstart += dsw.size();
const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null;
varstart += dxvar.size();
// Extract well parts np phase rates + bhp
const V dwells = subset(dx, Span(wellModel().numWellVars(), 1, varstart));
varstart += dwells.size();
assert(varstart == dx.size());
// Pressure update.
const double dpmaxrel = dpMaxRel();
const V p_old = Eigen::Map<const V>(&reservoir_state.pressure()[0], nc, 1);
const V absdpmax = dpmaxrel*p_old.abs();
const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
const V p = (p_old - dp_limited).max(zero);
std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
// Saturation updates.
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
const double dsmax = dsMax();
V so;
V sw;
V sg;
{
V maxVal = zero;
V dso = zero;
if (active_[Water]){
maxVal = dsw.abs().max(maxVal);
dso = dso - dsw;
}
V dsg;
if (active_[Gas]){
dsg = isSg_ * dxvar - isRv_ * dsw;
maxVal = dsg.abs().max(maxVal);
dso = dso - dsg;
}
maxVal = dso.abs().max(maxVal);
V step = dsmax/maxVal;
step = step.min(1.);
if (active_[Water]) {
const int pos = pu.phase_pos[ Water ];
const V sw_old = s_old.col(pos);
sw = sw_old - step * dsw;
}
if (active_[Gas]) {
const int pos = pu.phase_pos[ Gas ];
const V sg_old = s_old.col(pos);
sg = sg_old - step * dsg;
}
assert(active_[Oil]);
const int pos = pu.phase_pos[ Oil ];
const V so_old = s_old.col(pos);
so = so_old - step * dso;
}
// Appleyard chop process.
if (active_[Gas]) {
auto ixg = sg < 0;
for (int c = 0; c < nc; ++c) {
if (ixg[c]) {
if (active_[Water]) {
sw[c] = sw[c] / (1-sg[c]);
}
so[c] = so[c] / (1-sg[c]);
sg[c] = 0;
}
}
}
if (active_[Oil]) {
auto ixo = so < 0;
for (int c = 0; c < nc; ++c) {
if (ixo[c]) {
if (active_[Water]) {
sw[c] = sw[c] / (1-so[c]);
}
if (active_[Gas]) {
sg[c] = sg[c] / (1-so[c]);
}
so[c] = 0;
}
}
}
if (active_[Water]) {
auto ixw = sw < 0;
for (int c = 0; c < nc; ++c) {
if (ixw[c]) {
so[c] = so[c] / (1-sw[c]);
if (active_[Gas]) {
sg[c] = sg[c] / (1-sw[c]);
}
sw[c] = 0;
}
}
}
//const V sumSat = sw + so + sg;
//sw = sw / sumSat;
//so = so / sumSat;
//sg = sg / sumSat;
// Update the reservoir_state
if (active_[Water]) {
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
}
}
if (active_[Gas]) {
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
}
}
if (active_[ Oil ]) {
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pu.phase_pos[ Oil ]] = so[c];
}
}
// Update rs and rv
const double drmaxrel = drMaxRel();
V rs;
if (has_disgas_) {
const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
const V drs = isRs_ * dxvar;
const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drmaxrel);
rs = rs_old - drs_limited;
}
V rv;
if (has_vapoil_) {
const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
const V drv = isRv_ * dxvar;
const V drv_limited = sign(drv) * drv.abs().min(rv_old.abs()*drmaxrel);
rv = rv_old - drv_limited;
}
// Sg is used as primal variable for water only cells.
const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
auto watOnly = sw > (1 - epsilon);
// phase translation sg <-> rs
std::vector<HydroCarbonState>& hydroCarbonState = reservoir_state.hydroCarbonState();
std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::GasAndOil);
if (has_disgas_) {
const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
const V rsSat = fluidRsSat(p, so, cells_);
// The obvious case
auto hasGas = (sg > 0 && isRs_ == 0);
// Set oil saturated if previous rs is sufficiently large
const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs_ == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
auto useSg = watOnly || hasGas || gasVaporized;
for (int c = 0; c < nc; ++c) {
if (useSg[c]) {
rs[c] = rsSat[c];
} else {
hydroCarbonState[c] = HydroCarbonState::OilOnly;
}
}
}
// phase transitions so <-> rv
if (has_vapoil_) {
// The gas pressure is needed for the rvSat calculations
const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas));
const V gaspress = computeGasPressure(p, sw, so, sg);
const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_);
const V rvSat = fluidRvSat(gaspress, so, cells_);
// The obvious case
auto hasOil = (so > 0 && isRv_ == 0);
// Set oil saturated if previous rv is sufficiently large
const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv_ == 1) && (rv_old > rvSat0 * (1-epsilon)) );
auto useSg = watOnly || hasOil || oilCondensed;
for (int c = 0; c < nc; ++c) {
if (useSg[c]) {
rv[c] = rvSat[c];
} else {
hydroCarbonState[c] = HydroCarbonState::GasOnly;
}
}
}
// Update the reservoir_state
if (has_disgas_) {
std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
}
if (has_vapoil_) {
std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
}
wellModel().updateWellState(dwells, well_state);
// Update phase conditions used for property calculations.
updatePhaseCondFromPrimalVariable(reservoir_state);
}
/// Return true if output to cout is wanted. /// Return true if output to cout is wanted.
bool terminalOutputEnabled() const bool terminalOutputEnabled() const
@@ -1008,7 +703,6 @@ namespace Opm {
const Grid& grid_; const Grid& grid_;
const BlackoilPropsAdInterface& fluid_; const BlackoilPropsAdInterface& fluid_;
const DerivedGeology& geo_; const DerivedGeology& geo_;
const RockCompressibility* rock_comp_props_;
VFPProperties vfp_properties_; VFPProperties vfp_properties_;
const NewtonIterationBlackoilInterface& linsolver_; const NewtonIterationBlackoilInterface& linsolver_;
// For each canonical phase -> true if active // For each canonical phase -> true if active
@@ -1019,15 +713,10 @@ namespace Opm {
const bool has_vapoil_; const bool has_vapoil_;
ModelParameters param_; ModelParameters param_;
std::vector<PhasePresence> phaseCondition_;
// Well Model // Well Model
StandardWellsDense<FluidSystem, BlackoilIndices> well_model_; StandardWellsDense<FluidSystem, BlackoilIndices> well_model_;
V isRs_;
V isRv_;
V isSg_;
/// \brief Whether we print something to std::cout /// \brief Whether we print something to std::cout
bool terminal_output_; bool terminal_output_;
/// \brief The number of cells of the global grid. /// \brief The number of cells of the global grid.
@@ -1058,22 +747,6 @@ namespace Opm {
bool localWellsActive() const { return well_model_.localWellsActive(); } bool localWellsActive() const { return well_model_.localWellsActive(); }
V
fluidRsSat(const V& p,
const V& satOil,
const std::vector<int>& cells) const
{
return fluid_.rsSat(ADB::constant(p), ADB::constant(satOil), cells).value();
}
V
fluidRvSat(const V& p,
const V& satOil,
const std::vector<int>& cells) const
{
return fluid_.rvSat(ADB::constant(p), ADB::constant(satOil), cells).value();
}
void convertInput( const int iterationIdx, void convertInput( const int iterationIdx,
const ReservoirState& reservoirState, const ReservoirState& reservoirState,
Simulator& simulator ) const Simulator& simulator ) const
@@ -1095,13 +768,13 @@ namespace Opm {
cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]]; cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
// set switching variable and interpretation // set switching variable and interpretation
if( isRs_[ cellIdx ] && has_disgas_ ) if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ )
{ {
cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx]; cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx];
cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx]; cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs ); cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs );
} }
else if( isRv_[ cellIdx ] && has_vapoil_ ) else if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasOnly && has_vapoil_ )
{ {
// this case (-> gas only with vaporized oil in the gas) is // this case (-> gas only with vaporized oil in the gas) is
// relatively expensive as it requires to compute the capillary // relatively expensive as it requires to compute the capillary
@@ -1141,7 +814,7 @@ namespace Opm {
} }
else else
{ {
assert(isSg_[cellIdx]); assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil);
cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]]; cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]];
cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ]; cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ];
cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg ); cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg );
@@ -1186,8 +859,7 @@ namespace Opm {
{ {
const int numPhases = wells().number_of_phases; const int numPhases = wells().number_of_phases;
const int numCells = ebosJac.N(); const int numCells = ebosJac.N();
const int cols = ebosJac.M(); assert( numCells == ebosJac.M());
assert( numCells == cols );
// write the right-hand-side values from the ebosJac into the objects // write the right-hand-side values from the ebosJac into the objects
// allocated above. // allocated above.
@@ -1299,72 +971,6 @@ namespace Opm {
} }
V computeGasPressure(const V& po,
const V& sw,
const V& so,
const V& sg) const
{
assert (active_[Gas]);
std::vector<ADB> cp = fluid_.capPress(ADB::constant(sw),
ADB::constant(so),
ADB::constant(sg),
cells_);
return cp[Gas].value() + po;
}
const std::vector<PhasePresence>
phaseCondition() const {return phaseCondition_;}
/// update the primal variable for Sg, Rv or Rs. The Gas phase must
/// be active to call this method.
void
updatePrimalVariableFromState(const ReservoirState& state)
{
updatePhaseCondFromPrimalVariable(state);
}
/// Update the phaseCondition_ member based on the primalVariable_ member.
/// Also updates isRs_, isRv_ and isSg_;
void
updatePhaseCondFromPrimalVariable(const ReservoirState& state)
{
const int nc = Opm::AutoDiffGrid::numCells(grid_);
isRs_ = V::Zero(nc);
isRv_ = V::Zero(nc);
isSg_ = V::Zero(nc);
if (! (active_[Gas] && active_[Oil])) {
// updatePhaseCondFromPrimarVariable() logic requires active gas and oil phase.
phaseCondition_.assign(nc, PhasePresence());
return;
}
for (int c = 0; c < nc; ++c) {
phaseCondition_[c] = PhasePresence(); // No free phases.
phaseCondition_[c].setFreeWater(); // Not necessary for property calculation usage.
switch (state.hydroCarbonState()[c]) {
case HydroCarbonState::GasAndOil:
phaseCondition_[c].setFreeOil();
phaseCondition_[c].setFreeGas();
isSg_[c] = 1;
break;
case HydroCarbonState::OilOnly:
phaseCondition_[c].setFreeOil();
isRs_[c] = 1;
break;
case HydroCarbonState::GasOnly:
phaseCondition_[c].setFreeGas();
isRv_[c] = 1;
break;
default:
OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << c << ": " << state.hydroCarbonState()[c]);
}
}
}
double dpMaxRel() const { return param_.dp_max_rel_; } double dpMaxRel() const { return param_.dp_max_rel_; }
double dsMax() const { return param_.ds_max_; } double dsMax() const { return param_.ds_max_; }
double drMaxRel() const { return param_.dr_max_rel_; } double drMaxRel() const { return param_.dr_max_rel_; }

12
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@@ -174,13 +174,13 @@ public:
if( ! restorefilename.empty() ) if( ! restorefilename.empty() )
{ {
// -1 means that we'll take the last report step that was written // -1 means that we'll take the last report step that was written
const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) ); //const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
// output_writer_.restore( timer, // output_writer_.restore( timer,
// state, // state,
// prev_well_state, // prev_well_state,
// restorefilename, // restorefilename,
// desiredRestoreStep ); // desiredRestoreStep );
} }
unsigned int totalLinearizations = 0; unsigned int totalLinearizations = 0;

398
opm/autodiff/StandardWellsDense.hpp
View File

@@ -63,18 +63,6 @@ namespace Opm {
template<typename FluidSystem, typename BlackoilIndices> template<typename FluidSystem, typename BlackoilIndices>
class StandardWellsDense { class StandardWellsDense {
public: public:
struct WellOps {
WellOps(const Wells* wells)
: well_cells()
{
if( wells )
{
well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]);
}
}
std::vector<int> well_cells; // the set of perforated cells
};
// --------- Types --------- // --------- Types ---------
using ADB = AutoDiffBlock<double>; using ADB = AutoDiffBlock<double>;
@@ -109,11 +97,9 @@ namespace Opm {
, wells_(wells_arg) , wells_(wells_arg)
, fluid_(nullptr) , fluid_(nullptr)
, active_(nullptr) , active_(nullptr)
, phase_condition_(nullptr)
, vfp_properties_(nullptr) , vfp_properties_(nullptr)
, well_perforation_densities_(Vector()) , well_perforation_densities_(wells_arg->well_connpos[wells_arg->number_of_wells])
, well_perforation_pressure_diffs_(Vector()) , well_perforation_pressure_diffs_(wells_arg->well_connpos[wells_arg->number_of_wells])
, store_well_perforation_fluxes_(false)
, wellVariables_(wells_arg->number_of_wells * wells_arg->number_of_phases) , wellVariables_(wells_arg->number_of_wells * wells_arg->number_of_phases)
, F0_(wells_arg->number_of_wells * wells_arg->number_of_phases) , F0_(wells_arg->number_of_wells * wells_arg->number_of_phases)
, param_(param) , param_(param)
@@ -178,7 +164,6 @@ namespace Opm {
BVector rhs(nc); BVector rhs(nc);
BVector resWell(nw); BVector resWell(nw);
const V& cdp = wellPerforationPressureDiffs();
const double volume = 0.002831684659200; // 0.1 cu ft; const double volume = 0.002831684659200; // 0.1 cu ft;
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
@@ -187,7 +172,7 @@ namespace Opm {
const int cell_idx = wells().well_cells[perf]; const int cell_idx = wells().well_cells[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
std::vector<EvalWell> cq_s(np); std::vector<EvalWell> cq_s(np);
computeWellFlux(w, wells().WI[perf], intQuants, cdp[perf], cq_s); computeWellFlux(w, wells().WI[perf], intQuants, wellPerforationPressureDiffs()[perf], cq_s);
for (int p1 = 0; p1 < np; ++p1) { for (int p1 = 0; p1 < np; ++p1) {
// subtract sum of phase fluxes in the reservoir equation. // subtract sum of phase fluxes in the reservoir equation.
@@ -208,7 +193,7 @@ namespace Opm {
well_state.perfPhaseRates()[perf*np + p1] = cq_s[p1].value; well_state.perfPhaseRates()[perf*np + p1] = cq_s[p1].value;
} }
// Store the perforation pressure for later usage. // Store the perforation pressure for later usage.
well_state.perfPress()[perf] = well_state.bhp()[w] + cdp[perf]; well_state.perfPress()[perf] = well_state.bhp()[w] + wellPerforationPressureDiffs()[perf];
} }
// add vol * dF/dt + Q to the well equations; // add vol * dF/dt + Q to the well equations;
@@ -369,7 +354,6 @@ namespace Opm {
void init(const BlackoilPropsAdInterface* fluid_arg, void init(const BlackoilPropsAdInterface* fluid_arg,
const std::vector<bool>* active_arg, const std::vector<bool>* active_arg,
const std::vector<PhasePresence>* pc_arg,
const VFPProperties* vfp_properties_arg, const VFPProperties* vfp_properties_arg,
const double gravity_arg, const double gravity_arg,
const std::vector<double>& depth_arg, const std::vector<double>& depth_arg,
@@ -377,7 +361,6 @@ namespace Opm {
{ {
fluid_ = fluid_arg; fluid_ = fluid_arg;
active_ = active_arg; active_ = active_arg;
phase_condition_ = pc_arg;
vfp_properties_ = vfp_properties_arg; vfp_properties_ = vfp_properties_arg;
gravity_ = gravity_arg; gravity_ = gravity_arg;
cell_depths_ = extractPerfData(depth_arg); cell_depths_ = extractPerfData(depth_arg);
@@ -450,20 +433,12 @@ namespace Opm {
} }
/// Density of each well perforation /// Density of each well perforation
Vector& wellPerforationDensities() // mutable version kept for BlackoilMultisegmentModel std::vector<double> wellPerforationDensities() const {
{
return well_perforation_densities_;
}
const Vector& wellPerforationDensities() const {
return well_perforation_densities_; return well_perforation_densities_;
} }
/// Diff to bhp for each well perforation. /// Diff to bhp for each well perforation.
Vector& wellPerforationPressureDiffs() { // mutable version kept for BlackoilMultisegmentModel std::vector<double> wellPerforationPressureDiffs() const
return well_perforation_pressure_diffs_;
}
const Vector& wellPerforationPressureDiffs() const
{ {
return well_perforation_pressure_diffs_; return well_perforation_pressure_diffs_;
} }
@@ -648,7 +623,6 @@ namespace Opm {
const double dt, const double dt,
WellState& well_state) WellState& well_state)
{ {
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells; const int nw = wells().number_of_wells;
WellState well_state0 = well_state; WellState well_state0 = well_state;
@@ -665,28 +639,10 @@ namespace Opm {
++it; ++it;
if( localWellsActive() ) if( localWellsActive() )
{ {
BVector dx_new (nw); BVector dx_well (nw);
invDuneD_.mv(resWell_, dx_new); invDuneD_.mv(resWell_, dx_well);
V dx_new_eigen(np*nw);
for( int p=0; p<np; ++p) {
for (int w = 0; w < nw; ++w) {
int idx = w + nw*p;
dx_new_eigen(idx) = dx_new[w][flowPhaseToEbosCompIdx(p)];
}
}
//auto well_state2 = well_state;
//updateWellState(dx_new_eigen.array(), well_state);
updateWellState(dx_new, well_state);
// for (int w = 0; w < nw; ++w) {
// printIf(w, well_state.wellSolutions()[w], well_state2.wellSolutions()[w], 1e-3 ,"well solution 1");
// printIf(w, well_state.wellSolutions()[nw + w], well_state2.wellSolutions()[nw + w], 1e-3 ,"well solution 2");
// printIf(w, well_state.wellSolutions()[2*nw + w], well_state2.wellSolutions()[2*nw + w], 1e-3 ,"well solution 3");
// }
updateWellState(dx_well, well_state);
updateWellControls(well_state); updateWellControls(well_state);
setWellVariables(well_state); setWellVariables(well_state);
} }
@@ -858,28 +814,28 @@ namespace Opm {
FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
} }
const PhasePresence& perf_cond = (*phase_condition_)[wells().well_cells[perf]]; //const PhasePresence& perf_cond = (*phase_condition_)[wells().well_cells[perf]];
if (pu.phase_used[BlackoilPhases::Vapour]) { if (pu.phase_used[BlackoilPhases::Vapour]) {
int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases; int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases;
if (perf_cond.hasFreeOil()) { //if (perf_cond.hasFreeOil()) {
b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
} //}
else { // else {
double rv = fs.Rv().value; // double rv = fs.Rv().value;
b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); // b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv);
} // }
} }
if (pu.phase_used[BlackoilPhases::Liquid]) { if (pu.phase_used[BlackoilPhases::Liquid]) {
int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases;
if (perf_cond.hasFreeGas()) { //if (perf_cond.hasFreeGas()) {
b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
} //}
else { //else {
double rs = fs.Rs().value; // double rs = fs.Rs().value;
b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); // b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs);
} //}
} }
if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
@@ -906,7 +862,7 @@ namespace Opm {
double dFLimit = 0.2; double dFLimit = 0.2;
double dBHPLimit = 2.0; double dBHPLimit = 2.0;
const Vector xvar_well_old = Eigen::Map<const Vector>(&well_state.wellSolutions()[0], nw*np); std::vector<double> xvar_well_old = well_state.wellSolutions();
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
@@ -971,11 +927,6 @@ namespace Opm {
F[p] /= g[p]; F[p] /= g[p];
} }
} }
//std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl;
// const double dFw = dxvar_well[nw + w]; // const double dFw = dxvar_well[nw + w];
@@ -1099,17 +1050,18 @@ namespace Opm {
int table_id = well_controls_iget_vfp(wc, ctrl_index); int table_id = well_controls_iget_vfp(wc, ctrl_index);
const WellType& well_type = wells().type[w]; const WellType& well_type = wells().type[w];
const int perf = wells().well_connpos[w];
if (well_type == INJECTOR) { if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(), wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity_); wellPerforationDensities()[perf], gravity_);
well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp); well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp);
} }
else if (well_type == PRODUCER) { else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(), wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity_); wellPerforationDensities()[perf], gravity_);
well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp, alq); well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp, alq);
} }
@@ -1127,253 +1079,6 @@ namespace Opm {
template <class WellState>
void updateWellState(const Vector& dwells,
WellState& well_state)
{
if( localWellsActive() )
{
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
// Extract parts of dwells corresponding to each part.
int varstart = 0;
const Vector dxvar_well = subset(dwells, Span(np*nw, 1, varstart));
//const Vector dqs = subset(dwells, Span(np*nw, 1, varstart));
varstart += dxvar_well.size();
//const Vector dbhp = subset(dwells, Span(nw, 1, varstart));
//varstart += dbhp.size();
assert(varstart == dwells.size());
// Qs update.
// Since we need to update the wellrates, that are ordered by wells,
// from dqs which are ordered by phase, the simplest is to compute
// dwr, which is the data from dqs but ordered by wells.
const Vector xvar_well_old = Eigen::Map<const Vector>(&well_state.wellSolutions()[0], nw*np);
double dFLimit = 0.2;
double dBHPLimit = 2;
double dTotalRateLimit = 0.5;
//std::cout << "dxvar_well "<<dxvar_well << std::endl;
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
const int current = well_state.currentControls()[w];
const double target_rate = well_controls_iget_target(wc, current);
const double* distr = well_controls_iget_distr(wc, current);
std::vector<double> F(np,0.0);
const int sign2 = dxvar_well[nw + w] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dxvar_well[nw + w]),dFLimit);
well_state.wellSolutions()[nw + w] = xvar_well_old[nw + w] - dx2_limited;
const int sign3 = dxvar_well[2*nw + w] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dxvar_well[2*nw + w]),dFLimit);
well_state.wellSolutions()[2*nw + w] = xvar_well_old[2*nw + w] - dx3_limited;
F[Water] = well_state.wellSolutions()[nw + w];
F[Gas] = well_state.wellSolutions()[2*nw + w];
F[Oil] = 1.0 - F[Water] - F[Gas];
// const double dFw = dxvar_well[nw + w];
// const double dFg = dxvar_well[nw*2 + w];
// const double dFo = - dFw - dFg;
// //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl;
// double step = dFLimit / std::max(std::abs(dFw),std::max(std::abs(dFg),std::abs(dFo))); //)) / dFLimit;
// step = std::min(step, 1.0);
// //std::cout << step << std::endl;
// F[Water] = xvar_well_old[nw + w] - step*dFw;
// F[Gas] = xvar_well_old[2*nw + w] - step*dFg;
// F[Oil] = (1.0 - xvar_well_old[2*nw + w] - xvar_well_old[nw + w]) - step * dFo;
if (F[Water] < 0.0) {
F[Gas] /= (1.0 - F[Water]);
F[Oil] /= (1.0 - F[Water]);
F[Water] = 0.0;
}
if (F[Gas] < 0.0) {
F[Water] /= (1.0 - F[Gas]);
F[Oil] /= (1.0 - F[Gas]);
F[Gas] = 0.0;
}
if (F[Oil] < 0.0) {
F[Water] /= (1.0 - F[Oil]);
F[Gas] /= (1.0 - F[Oil]);
F[Oil] = 0.0;
}
well_state.wellSolutions()[nw + w] = F[Water];
well_state.wellSolutions()[2*nw + w] = F[Gas];
//std::cout << wells().name[w] << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl;
std::vector<double> g = {1,1,0.01};
if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) {
for (int p = 0; p < np; ++p) {
F[p] /= distr[p];
}
} else {
for (int p = 0; p < np; ++p) {
F[p] /= g[p];
}
}
//std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl;
// const double dFw = dxvar_well[nw + w];
// const double dFg = dxvar_well[nw*2 + w];
// const double dFo = - dFw - dFg;
//std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl;
// double step = dFLimit / std::max(std::abs(dFw),std::max(std::abs(dFg),std::abs(dFo))); //)) / dFLimit;
// step = std::min(step, 1.0);
// std::cout << step << std::endl;
// F[Water] = xvar_well_old[nw + w] - step*dFw;
// F[Gas] = xvar_well_old[2*nw + w] - step*dFg;
// F[Oil] = (1.0 - xvar_well_old[2*nw + w] - xvar_well_old[nw + w]) - step * dFo;
// double sumF = F[Water]+F[Gas]+F[Oil];
// F[Water] /= sumF;
// F[Gas] /= sumF;
// F[Oil] /= sumF;
// well_state.wellSolutions()[nw + w] = F[Water];
// well_state.wellSolutions()[2 * nw + w] = F[Gas];
switch (well_controls_iget_type(wc, current)) {
case BHP:
{
//const int sign1 = dxvar_well[w] > 0 ? 1: -1;
//const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dTotalRateLimit);
well_state.wellSolutions()[w] = xvar_well_old[w] - dxvar_well[w];
switch (wells().type[w]) {
case INJECTOR:
for (int p = 0; p < np; ++p) {
const double comp_frac = wells().comp_frac[np*w + p];
//if (comp_frac > 0) {
well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[w];
//}
}
break;
case PRODUCER:
for (int p = 0; p < np; ++p) {
well_state.wellRates()[w*np + p] = well_state.wellSolutions()[w] * F[p];
}
break;
}
}
break;
case SURFACE_RATE:
{
const int sign1 = dxvar_well[w] > 0 ? 1: -1;
const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dBHPLimit);
well_state.wellSolutions()[w] = xvar_well_old[w] - dx1_limited;
//const int sign = (dxvar_well1[w] < 0) ? -1 : 1;
//well_state.bhp()[w] -= sign * std::min( std::abs(dxvar_well1[w]), std::abs(well_state.bhp()[w])*dpmaxrel) ;
well_state.bhp()[w] = well_state.wellSolutions()[w];
if (wells().type[w]==PRODUCER) {
double F_target = 0.0;
for (int p = 0; p < np; ++p) {
F_target += wells().comp_frac[np*w + p] * F[p];
}
for (int p = 0; p < np; ++p) {
//std::cout << F[p] << std::endl;
well_state.wellRates()[np*w + p] = F[p] * target_rate /F_target;
}
} else {
for (int p = 0; p < np; ++p) {
//std::cout << wells().comp_frac[np*w + p] << " " <<distr[p] << std::endl;
well_state.wellRates()[w*np + p] = wells().comp_frac[np*w + p] * target_rate;
}
}
}
break;
case RESERVOIR_RATE: {
const int sign1 = dxvar_well[w] > 0 ? 1: -1;
const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dBHPLimit);
well_state.wellSolutions()[w] = xvar_well_old[w] - dx1_limited;
//const int sign = (dxvar_well1[w] < 0) ? -1 : 1;
//well_state.bhp()[w] -= sign * std::min( std::abs(dxvar_well1[w]), std::abs(well_state.bhp()[w])*dpmaxrel) ;
well_state.bhp()[w] = well_state.wellSolutions()[w];
for (int p = 0; p < np; ++p) {
well_state.wellRates()[np*w + p] = F[p] * target_rate;
}
}
break;
}
}
const Opm::PhaseUsage& pu = fluid_->phaseUsage();
//Loop over all wells
#pragma omp parallel for schedule(static)
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
const int nwc = well_controls_get_num(wc);
//Loop over all controls until we find a THP control
//that specifies what we need...
//Will only update THP for wells with THP control
for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(wc, 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 ] ];
}
double alq = well_controls_iget_alq(wc, ctrl_index);
int table_id = well_controls_iget_vfp(wc, ctrl_index);
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_);
well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp);
}
else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity_);
well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp, alq);
}
else {
OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
}
//Assume only one THP control specified for each well
break;
}
}
}
}
}
template <class WellState> template <class WellState>
void updateWellControls(WellState& xw) void updateWellControls(WellState& xw)
{ {
@@ -1458,18 +1163,19 @@ namespace Opm {
//Set *BHP* target by calculating bhp from THP //Set *BHP* target by calculating bhp from THP
const WellType& well_type = wells().type[w]; const WellType& well_type = wells().type[w];
const int perf = wells().well_connpos[w];
if (well_type == INJECTOR) { if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity_); wellPerforationDensities()[perf], 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) { else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity_); wellPerforationDensities()[perf], 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;
} }
@@ -1569,24 +1275,6 @@ namespace Opm {
} }
/// If set, computeWellFlux() will additionally store the
/// total reservoir volume perforation fluxes.
void setStoreWellPerforationFluxesFlag(const bool store_fluxes)
{
store_well_perforation_fluxes_ = store_fluxes;
}
/// Retrieves the stored fluxes. It is an error to call this
/// unless setStoreWellPerforationFluxesFlag(true) has been
/// called.
const Vector& getStoredWellPerforationFluxes() const
{
assert(store_well_perforation_fluxes_);
return well_perforation_fluxes_;
}
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
{ {
const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx };
@@ -1705,21 +1393,17 @@ namespace Opm {
const std::vector<double>& depth_perf, const std::vector<double>& depth_perf,
const double grav) { const double grav) {
// Compute densities // Compute densities
std::vector<double> cd = well_perforation_densities_ =
WellDensitySegmented::computeConnectionDensities( WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_->phaseUsage(), wells(), xw, fluid_->phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
const int nperf = wells().well_connpos[wells().number_of_wells];
// Compute pressure deltas // Compute pressure deltas
std::vector<double> cdp = well_perforation_pressure_diffs_ =
WellDensitySegmented::computeConnectionPressureDelta( WellDensitySegmented::computeConnectionPressureDelta(
wells(), depth_perf, cd, grav); wells(), depth_perf, well_perforation_densities_, grav);
// Store the results
well_perforation_densities_ = Eigen::Map<const Vector>(cd.data(), nperf);
well_perforation_pressure_diffs_ = Eigen::Map<const Vector>(cdp.data(), nperf);
} }
@@ -1731,7 +1415,6 @@ namespace Opm {
const BlackoilPropsAdInterface* fluid_; const BlackoilPropsAdInterface* fluid_;
const std::vector<bool>* active_; const std::vector<bool>* active_;
const std::vector<PhasePresence>* phase_condition_;
const VFPProperties* vfp_properties_; const VFPProperties* vfp_properties_;
double gravity_; double gravity_;
// the depth of the all the cell centers // the depth of the all the cell centers
@@ -1739,11 +1422,8 @@ namespace Opm {
std::vector<double> cell_depths_; std::vector<double> cell_depths_;
std::vector<double> pv_; std::vector<double> pv_;
Vector well_perforation_densities_; std::vector<double> well_perforation_densities_;
Vector well_perforation_pressure_diffs_; std::vector<double> well_perforation_pressure_diffs_;
bool store_well_perforation_fluxes_;
Vector well_perforation_fluxes_;
std::vector<EvalWell> wellVariables_; std::vector<EvalWell> wellVariables_;
std::vector<double> F0_; std::vector<double> F0_;
@@ -1770,7 +1450,7 @@ namespace Opm {
EvalWell getQs(const int wellIdx, const int phaseIdx) const { EvalWell getQs(const int wellIdx, const int phaseIdx) const {
EvalWell qs = 0.0; EvalWell qs = 0.0;
const WellControls* wc = wells().ctrls[wellIdx]; const WellControls* wc = wells().ctrls[wellIdx];
const int np = fluid_->numPhases(); const int np = wells().number_of_phases;
const double target_rate = well_controls_get_current_target(wc); const double target_rate = well_controls_get_current_target(wc);
if (wells().type[wellIdx] == INJECTOR) { if (wells().type[wellIdx] == INJECTOR) {
@@ -1811,7 +1491,7 @@ namespace Opm {
} }
EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const { EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const {
assert(fluid_.numPhases() == 3); assert(wells().number_of_phases == 3);
const int nw = wells().number_of_wells; const int nw = wells().number_of_wells;
if (phaseIdx == Water) { if (phaseIdx == Water) {
return wellVariables_[nw + wellIdx]; return wellVariables_[nw + wellIdx];

15
opm/autodiff/StandardWells_impl.hpp
View File

@@ -672,17 +672,18 @@ namespace Opm
int table_id = well_controls_iget_vfp(wc, ctrl_index); int table_id = well_controls_iget_vfp(wc, ctrl_index);
const WellType& well_type = wells().type[w]; const WellType& well_type = wells().type[w];
const int perf = wells().well_connpos[w]; //first perforation.
if (well_type == INJECTOR) { if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(), wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity_); wellPerforationDensities()[perf], 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) { else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(), wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity_); wellPerforationDensities()[perf], 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);
} }
@@ -786,17 +787,18 @@ namespace Opm
//Set *BHP* target by calculating bhp from THP //Set *BHP* target by calculating bhp from THP
const WellType& well_type = wells().type[w]; const WellType& well_type = wells().type[w];
const int perf = wells().well_connpos[w]; // first perforation.
if (well_type == INJECTOR) { if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity_); wellPerforationDensities()[perf], 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) { else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity_); wellPerforationDensities()[perf], 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;
} }
@@ -1096,11 +1098,12 @@ namespace Opm
//Set *BHP* target by calculating bhp from THP //Set *BHP* target by calculating bhp from THP
const WellType& well_type = wells().type[w]; const WellType& well_type = wells().type[w];
const int perf = wells().well_connpos[w]; //first perforation
if (well_type == INJECTOR) { if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity_); wellPerforationDensities()[perf], gravity_);
const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; 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 // apply the strictest of the bhp controlls i.e. smallest bhp for injectors
if ( bhp < bhps[w]) { if ( bhp < bhps[w]) {
@@ -1110,7 +1113,7 @@ namespace Opm
else if (well_type == PRODUCER) { else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity_); wellPerforationDensities()[perf], gravity_);
const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; 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 // apply the strictest of the bhp controlls i.e. largest bhp for producers

23
opm/autodiff/WellHelpers.hpp
View File

@@ -131,7 +131,7 @@ namespace Opm {
*/ */
inline inline
double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth, double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth,
const Vector& well_perforation_densities, const double gravity) { const double& rho, const double gravity) {
if ( wells.well_connpos[w] == wells.well_connpos[w+1] ) if ( wells.well_connpos[w] == wells.well_connpos[w+1] )
{ {
// This is a well with no perforations. // This is a well with no perforations.
@@ -142,13 +142,12 @@ namespace Opm {
} }
const double well_ref_depth = wells.depth_ref[w]; const double well_ref_depth = wells.depth_ref[w];
const double dh = vfp_ref_depth - well_ref_depth; const double dh = vfp_ref_depth - well_ref_depth;
const int perf = wells.well_connpos[w];
const double rho = well_perforation_densities[perf];
const double dp = rho*gravity*dh; const double dp = rho*gravity*dh;
return dp; return dp;
} }
inline inline
Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth, Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth,
const Vector& well_perforation_densities, const double gravity) { const Vector& well_perforation_densities, const double gravity) {
@@ -157,7 +156,23 @@ namespace Opm {
#pragma omp parallel for schedule(static) #pragma omp parallel for schedule(static)
for (int i=0; i<nw; ++i) { for (int i=0; i<nw; ++i) {
retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities, gravity); const int perf = wells.well_connpos[i];
retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities[perf], gravity);
}
return retval;
}
inline
std::vector<double> computeHydrostaticCorrection(const Wells& wells, const std::vector<double>& vfp_ref_depth,
const std::vector<double>& well_perforation_densities, const double gravity) {
const int nw = wells.number_of_wells;
std::vector<double> retval(nw,0.0);
#pragma omp parallel for schedule(static)
for (int i=0; i<nw; ++i) {
const int perf = wells.well_connpos[i];
retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities[perf], gravity);
} }
return retval; return retval;