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 #647 from GitPaean/wells_refactoring_onlytemplate_function

Browse Source 
Moving more well related functions to StandardWells
This commit is contained in:
dr-robertk 2016-04-16 12:37:42 +02:00
commit 93f48d6ad9
15 changed files with 1502 additions and 1003 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists_files.cmake
View File

@ -47,7 +47,6 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/VFPProdProperties.cpp opm/autodiff/VFPProdProperties.cpp
opm/autodiff/VFPInjProperties.cpp opm/autodiff/VFPInjProperties.cpp
opm/autodiff/WellMultiSegment.cpp opm/autodiff/WellMultiSegment.cpp
opm/autodiff/StandardWells.cpp
opm/autodiff/BlackoilSolventState.cpp opm/autodiff/BlackoilSolventState.cpp
opm/autodiff/ThreadHandle.hpp opm/autodiff/ThreadHandle.hpp
opm/polymer/PolymerState.cpp opm/polymer/PolymerState.cpp
@ -194,7 +193,9 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/VFPInjProperties.hpp opm/autodiff/VFPInjProperties.hpp
opm/autodiff/WellStateMultiSegment.hpp opm/autodiff/WellStateMultiSegment.hpp
opm/autodiff/WellMultiSegment.hpp opm/autodiff/WellMultiSegment.hpp
opm/autodiff/WellHelpers.hpp
opm/autodiff/StandardWells.hpp opm/autodiff/StandardWells.hpp
opm/autodiff/StandardWellsSolvent.hpp
opm/polymer/CompressibleTpfaPolymer.hpp opm/polymer/CompressibleTpfaPolymer.hpp
opm/polymer/GravityColumnSolverPolymer.hpp opm/polymer/GravityColumnSolverPolymer.hpp
opm/polymer/GravityColumnSolverPolymer_impl.hpp opm/polymer/GravityColumnSolverPolymer_impl.hpp

32
opm/autodiff/BlackoilModelBase.hpp
View File

@ -383,39 +383,28 @@ namespace Opm {
void computeWellConnectionPressures(const SolutionState& state, void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw); const WellState& xw);
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);
void void
assembleMassBalanceEq(const SolutionState& state); assembleMassBalanceEq(const SolutionState& state);
// TODO: only kept for now due to flow_multisegment
// will be removed soon
void void
extractWellPerfProperties(const SolutionState& state, extractWellPerfProperties(const SolutionState& state,
std::vector<ADB>& mob_perfcells, std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const; 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 bool
solveWellEq(const std::vector<ADB>& mob_perfcells, solveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells, const std::vector<ADB>& b_perfcells,
SolutionState& state, SolutionState& state,
WellState& well_state); WellState& well_state);
void
computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& aliveWells,
std::vector<ADB>& cq_s) const;
void
updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw) const;
void void
addWellFluxEq(const std::vector<ADB>& cq_s, addWellFluxEq(const std::vector<ADB>& cq_s,
const SolutionState& state); const SolutionState& state);
@ -430,11 +419,6 @@ namespace Opm {
const WellState& xw, const WellState& xw,
const V& aliveWells); const V& aliveWells);
void updateWellControls(WellState& xw) const;
void updateWellState(const V& dwells,
WellState& well_state);
bool getWellConvergence(const int iteration); bool getWellConvergence(const int iteration);
bool isVFPActive() const; bool isVFPActive() const;

714
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -29,6 +29,7 @@
#include <opm/autodiff/AutoDiffBlock.hpp> #include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp> #include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/GridHelpers.hpp> #include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/WellHelpers.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp> #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/GeoProps.hpp> #include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp> #include <opm/autodiff/WellDensitySegmented.hpp>
@ -758,81 +759,8 @@ namespace detail {
} }
} }
template <class Grid, class Implementation>
void BlackoilModelBase<Grid, Implementation>::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)
{
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
// Compute the average pressure in each well block
const V perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
V avg_press = perf_press*0;
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
}
const std::vector<int>& well_cells = stdWells().wellOps().well_cells;
// Use cell values for the temperature as the wells don't knows its temperature yet.
const ADB perf_temp = subset(state.temperature, well_cells);
// Compute b, rsmax, rvmax values for perforations.
// Evaluate the properties using average well block pressures
// and cell values for rs, rv, phase condition and temperature.
const ADB avg_press_ad = ADB::constant(avg_press);
std::vector<PhasePresence> perf_cond(nperf);
const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf; ++perf) {
perf_cond[perf] = pc[well_cells[perf]];
}
const PhaseUsage& pu = fluid_.phaseUsage();
DataBlock b(nperf, pu.num_phases);
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);
} else {
rsmax_perf.assign(nperf, 0.0);
}
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);
} else {
rvmax_perf.assign(nperf, 0.0);
}
// b is row major, so can just copy data.
b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
// 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, 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);
}
surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases);
}
template <class Grid, class Implementation> template <class Grid, class Implementation>
@ -849,33 +777,19 @@ namespace detail {
std::vector<double> rsmax_perf; std::vector<double> rsmax_perf;
std::vector<double> rvmax_perf; std::vector<double> rvmax_perf;
std::vector<double> surf_dens_perf; std::vector<double> surf_dens_perf;
asImpl().computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); asImpl().stdWells().computePropertiesForWellConnectionPressures(state, xw, fluid_, active_, phaseCondition_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// 2. Compute densities
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_.phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
const int nperf = wells().well_connpos[wells().number_of_wells];
const std::vector<int>& well_cells = stdWells().wellOps().well_cells;
// Extract well connection depths. // Extract well connection depths.
const V depth = cellCentroidsZToEigen(grid_); const StandardWells::Vector depth = cellCentroidsZToEigen(grid_);
const V pdepth = subset(depth, well_cells); const StandardWells::Vector pdepth = subset(depth, asImpl().stdWells().wellOps().well_cells);
std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf); const int nperf = wells().well_connpos[wells().number_of_wells];
const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
// Gravity // Gravity
double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); const double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
// 3. Compute pressure deltas asImpl().stdWells().computeWellConnectionDensitesPressures(xw, fluid_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, grav);
std::vector<double> cdp =
WellDensitySegmented::computeConnectionPressureDelta(
wells(), perf_depth, cd, grav);
// 4. Store the results
stdWells().wellPerforationDensities() = Eigen::Map<const V>(cd.data(), nperf);
stdWells().wellPerforationPressureDiffs() = Eigen::Map<const V>(cdp.data(), nperf);
} }
@ -903,7 +817,10 @@ namespace detail {
// Possibly switch well controls and updating well state to // Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells // get reasonable initial conditions for the wells
asImpl().updateWellControls(well_state); // asImpl().updateWellControls(well_state);
// asImpl().stdWells().updateWellControls(well_state);
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
asImpl().stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state);
// Create the primary variables. // Create the primary variables.
SolutionState state = asImpl().variableState(reservoir_state, well_state); SolutionState state = asImpl().variableState(reservoir_state, well_state);
@ -938,15 +855,15 @@ namespace detail {
std::vector<ADB> mob_perfcells; std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells; std::vector<ADB> b_perfcells;
asImpl().extractWellPerfProperties(state, mob_perfcells, b_perfcells); asImpl().stdWells().extractWellPerfProperties(state, rq_, fluid_.numPhases(), fluid_, active_, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) { if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step // solve the well equations as a pre-processing step
asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state); asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
} }
V aliveWells; V aliveWells;
std::vector<ADB> cq_s; std::vector<ADB> cq_s;
asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); asImpl().stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s);
asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
asImpl().addWellFluxEq(cq_s, state); asImpl().addWellFluxEq(cq_s, state);
asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
asImpl().addWellControlEq(state, well_state, aliveWells); asImpl().addWellControlEq(state, well_state, aliveWells);
@ -1111,190 +1028,6 @@ namespace detail {
template <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& aliveWells,
std::vector<ADB>& cq_s) const
{
if( ! localWellsActive() ) return ;
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
V Tw = Eigen::Map<const V>(wells().WI, nperf);
const std::vector<int>& well_cells = stdWells().wellOps().well_cells;
// pressure diffs computed already (once per step, not changing per iteration)
const V& cdp = stdWells().wellPerforationPressureDiffs();
// Extract needed quantities for the perforation cells
const ADB& p_perfcells = subset(state.pressure, well_cells);
const ADB& rv_perfcells = subset(state.rv, well_cells);
const ADB& rs_perfcells = subset(state.rs, well_cells);
// Perforation pressure
const ADB perfpressure = (stdWells().wellOps().w2p * state.bhp) + cdp;
// Pressure drawdown (also used to determine direction of flow)
const ADB drawdown = p_perfcells - perfpressure;
// Compute vectors with zero and ones that
// selects the wanted quantities.
// selects injection perforations
V selectInjectingPerforations = V::Zero(nperf);
// selects producing perforations
V selectProducingPerforations = V::Zero(nperf);
for (int c = 0; c < nperf; ++c){
if (drawdown.value()[c] < 0)
selectInjectingPerforations[c] = 1;
else
selectProducingPerforations[c] = 1;
}
// Handle cross flow
const V numInjectingPerforations = (stdWells().wellOps().p2w * ADB::constant(selectInjectingPerforations)).value();
const V numProducingPerforations = (stdWells().wellOps().p2w * ADB::constant(selectProducingPerforations)).value();
for (int w = 0; w < nw; ++w) {
if (!wells().allow_cf[w]) {
for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) {
// Crossflow is not allowed; reverse flow is prevented.
// At least one of the perforation must be open in order to have a meeningful
// equation to solve. For the special case where all perforations have reverse flow,
// and the target rate is non-zero all of the perforations are keept open.
if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
selectProducingPerforations[perf] = 0.0;
} else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
selectInjectingPerforations[perf] = 0.0;
}
}
}
}
// HANDLE FLOW INTO WELLBORE
// compute phase volumetric rates at standard conditions
std::vector<ADB> cq_ps(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
cq_ps[phase] = b_perfcells[phase] * cq_p;
}
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];
const ADB cq_psGas = cq_ps[gaspos];
cq_ps[gaspos] += rs_perfcells * cq_psOil;
cq_ps[oilpos] += rv_perfcells * cq_psGas;
}
// HANDLE FLOW OUT FROM WELLBORE
// Using total mobilities
ADB total_mob = mob_perfcells[0];
for (int phase = 1; phase < np; ++phase) {
total_mob += mob_perfcells[phase];
}
// injection perforations total volume rates
const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
// compute wellbore mixture for injecting perforations
// The wellbore mixture depends on the inflow from the reservoar
// and the well injection rates.
// compute avg. and total wellbore phase volumetric rates at standard conds
const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
std::vector<ADB> wbq(np, ADB::null());
ADB wbqt = ADB::constant(V::Zero(nw));
for (int phase = 0; phase < np; ++phase) {
const ADB& q_ps = stdWells().wellOps().p2w * cq_ps[phase];
const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw));
Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
const int pos = pu.phase_pos[phase];
wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(V::Zero(nw)))) - q_ps;
wbqt += wbq[phase];
}
// compute wellbore mixture at standard conditions.
Selector<double> notDeadWells_selector(wbqt.value(), Selector<double>::Zero);
std::vector<ADB> cmix_s(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
const int pos = pu.phase_pos[phase];
cmix_s[phase] = stdWells().wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
}
// compute volume ratio between connection at standard conditions
ADB volumeRatio = ADB::constant(V::Zero(nperf));
const ADB d = V::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
for (int phase = 0; phase < np; ++phase) {
ADB tmp = cmix_s[phase];
if (phase == Oil && active_[Gas]) {
const int gaspos = pu.phase_pos[Gas];
tmp -= rv_perfcells * cmix_s[gaspos] / d;
}
if (phase == Gas && active_[Oil]) {
const int oilpos = pu.phase_pos[Oil];
tmp -= rs_perfcells * cmix_s[oilpos] / d;
}
volumeRatio += tmp / b_perfcells[phase];
}
// injecting connections total volumerates at standard conditions
ADB cqt_is = cqt_i/volumeRatio;
// connection phase volumerates at standard conditions
cq_s.resize(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
}
// check for dead wells (used in the well controll equations)
aliveWells = V::Constant(nw, 1.0);
for (int w = 0; w < nw; ++w) {
if (wbqt.value()[w] == 0) {
aliveWells[w] = 0.0;
}
}
}
template <class Grid, class Implementation>
void BlackoilModelBase<Grid, Implementation>::updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw) const
{
if ( !asImpl().localWellsActive() )
{
// If there are no wells in the subdomain of the proces then
// cq_s has zero size and will cause a segmentation fault below.
return;
}
// Update the perforation phase rates (used to calculate the pressure drop in the wellbore).
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np);
}
xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np);
// Update the perforation pressures.
const V& cdp = stdWells().wellPerforationPressureDiffs();
const V perfpressure = (stdWells().wellOps().w2p * state.bhp.value().matrix()).array() + cdp;
xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf);
}
template <class Grid, class Implementation> template <class Grid, class Implementation>
void BlackoilModelBase<Grid, Implementation>::addWellFluxEq(const std::vector<ADB>& cq_s, void BlackoilModelBase<Grid, Implementation>::addWellFluxEq(const std::vector<ADB>& cq_s,
const SolutionState& state) const SolutionState& state)
@ -1321,138 +1054,6 @@ namespace detail {
namespace detail
{
inline
double rateToCompare(const std::vector<double>& well_phase_flow_rate,
const int well,
const int num_phases,
const double* distr)
{
double rate = 0.0;
for (int phase = 0; phase < num_phases; ++phase) {
// Important: well_phase_flow_rate is ordered with all phase rates for first
// well first, then all phase rates for second well etc.
rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase];
}
return rate;
}
inline
bool constraintBroken(const std::vector<double>& bhp,
const std::vector<double>& thp,
const std::vector<double>& well_phase_flow_rate,
const int well,
const int num_phases,
const WellType& well_type,
const WellControls* wc,
const int ctrl_index)
{
const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
const double target = well_controls_iget_target(wc, ctrl_index);
const double* distr = well_controls_iget_distr(wc, ctrl_index);
bool broken = false;
switch (well_type) {
case INJECTOR:
{
switch (ctrl_type) {
case BHP:
broken = bhp[well] > target;
break;
case THP:
broken = thp[well] > target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
broken = rateToCompare(well_phase_flow_rate,
well, num_phases, distr) > target;
break;
}
}
break;
case PRODUCER:
{
switch (ctrl_type) {
case BHP:
broken = bhp[well] < target;
break;
case THP:
broken = thp[well] < target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
// Note that the rates compared below are negative,
// so breaking the constraints means: too high flow rate
// (as for injection).
broken = rateToCompare(well_phase_flow_rate,
well, num_phases, distr) < target;
break;
}
}
break;
default:
OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells.");
}
return broken;
}
} // namespace detail
namespace detail {
/**
* Simple hydrostatic correction for VFP table
* @param wells - wells struct
* @param w Well number
* @param vfp_table VFP table
* @param well_perforation_densities Densities at well perforations
* @param gravity Gravitational constant (e.g., 9.81...)
*/
inline
double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth,
const ADB::V& well_perforation_densities, const double gravity) {
if ( wells.well_connpos[w] == wells.well_connpos[w+1] )
{
// This is a well with no perforations.
// If this is the last well we would subscript over the
// bounds below.
// we assume well_perforation_densities to be 0
return 0;
}
const double well_ref_depth = wells.depth_ref[w];
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;
return dp;
}
inline
ADB::V computeHydrostaticCorrection(const Wells& wells, const ADB::V vfp_ref_depth,
const ADB::V& well_perforation_densities, const double gravity) {
const int nw = wells.number_of_wells;
ADB::V retval = ADB::V::Zero(nw);
#pragma omp parallel for schedule(static)
for (int i=0; i<nw; ++i) {
retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities, gravity);
}
return retval;
}
} //Namespace
template <class Grid, class Implementation> template <class Grid, class Implementation>
bool BlackoilModelBase<Grid, Implementation>::isVFPActive() const bool BlackoilModelBase<Grid, Implementation>::isVFPActive() const
{ {
@ -1485,156 +1086,6 @@ namespace detail {
} }
template <class Grid, class Implementation>
void BlackoilModelBase<Grid, Implementation>::updateWellControls(WellState& xw) const
{
if( ! localWellsActive() ) 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 = wells().number_of_phases;
const int nw = wells().number_of_wells;
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
#pragma omp parallel for schedule(dynamic)
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
int current = xw.currentControls()[w];
// Loop over all controls except the current one, and also
// skip any RESERVOIR_RATE controls, since we cannot
// handle those.
const int nwc = well_controls_get_num(wc);
int ctrl_index = 0;
for (; ctrl_index < nwc; ++ctrl_index) {
if (ctrl_index == current) {
// This is the currently used control, so it is
// used as an equation. So this is not used as an
// inequality constraint, and therefore skipped.
continue;
}
if (detail::constraintBroken(
xw.bhp(), xw.thp(), xw.wellRates(),
w, np, wells().type[w], wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
}
}
if (ctrl_index != nwc) {
// Constraint number ctrl_index was broken, switch to it.
if (terminal_output_)
{
std::cout << "Switching control mode for well " << wells().name[w]
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
}
xw.currentControls()[w] = ctrl_index;
current = xw.currentControls()[w];
}
//Get gravity for THP hydrostatic corrrection
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
// Updating well state and primary variables.
// Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
const double target = well_controls_iget_target(wc, current);
const double* distr = well_controls_iget_distr(wc, current);
switch (well_controls_iget_type(wc, current)) {
case BHP:
xw.bhp()[w] = target;
break;
case THP: {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
if (active_[ Water ]) {
aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if (active_[ Oil ]) {
liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if (active_[ Gas ]) {
vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ];
}
const int vfp = well_controls_iget_vfp(wc, current);
const double& thp = well_controls_iget_target(wc, current);
const double& alq = well_controls_iget_alq(wc, current);
//Set *BHP* target by calculating bhp from THP
const WellType& well_type = wells().type[w];
if (well_type == INJECTOR) {
double dp = detail::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(),
stdWells().wellPerforationDensities(), gravity);
xw.bhp()[w] = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
}
else if (well_type == PRODUCER) {
double dp = detail::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(),
stdWells().wellPerforationDensities(), gravity);
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");
}
break;
}
case RESERVOIR_RATE:
// No direct change to any observable quantity at
// surface condition. In this case, use existing
// flow rates as initial conditions as reservoir
// rate acts only in aggregate.
break;
case SURFACE_RATE:
// assign target value as initial guess for injectors and
// single phase producers (orat, grat, wrat)
const WellType& well_type = wells().type[w];
if (well_type == INJECTOR) {
for (int phase = 0; phase < np; ++phase) {
const double& compi = wells().comp_frac[np * w + phase];
if (compi > 0.0) {
xw.wellRates()[np*w + phase] = target * compi;
}
}
} else if (well_type == PRODUCER) {
// only set target as initial rates for single phase
// producers. (orat, grat and wrat, and not lrat)
// lrat will result in numPhasesWithTargetsUnderThisControl == 2
int numPhasesWithTargetsUnderThisControl = 0;
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0) {
numPhasesWithTargetsUnderThisControl += 1;
}
}
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) {
xw.wellRates()[np*w + phase] = target * distr[phase];
}
}
} else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
break;
}
}
}
template <class Grid, class Implementation> template <class Grid, class Implementation>
@ -1674,8 +1125,8 @@ namespace detail {
SolutionState wellSolutionState = state0; SolutionState wellSolutionState = state0;
asImpl().variableStateExtractWellsVars(indices, vars, wellSolutionState); asImpl().variableStateExtractWellsVars(indices, vars, wellSolutionState);
asImpl().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); asImpl().stdWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
asImpl().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
asImpl().addWellFluxEq(cq_s, wellSolutionState); asImpl().addWellFluxEq(cq_s, wellSolutionState);
asImpl().addWellControlEq(wellSolutionState, well_state, aliveWells); asImpl().addWellControlEq(wellSolutionState, well_state, aliveWells);
converged = getWellConvergence(it); converged = getWellConvergence(it);
@ -1700,8 +1151,10 @@ namespace detail {
ADB::V total_residual_v = total_residual.value(); ADB::V total_residual_v = total_residual.value();
const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
assert(dx.size() == total_residual_v.size()); assert(dx.size() == total_residual_v.size());
asImpl().updateWellState(dx.array(), well_state); // asImpl().updateWellState(dx.array(), well_state);
asImpl().updateWellControls(well_state); const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
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);
} }
} while (it < 15); } while (it < 15);
@ -1810,7 +1263,7 @@ namespace detail {
case THP: case THP:
{ {
const int perf = wells().well_connpos[w]; const int perf = wells().well_connpos[w];
rho_v[w] = stdWells().wellPerforationDensities()[perf]; rho_v[w] = asImpl().stdWells().wellPerforationDensities()[perf];
const int table_id = well_controls_iget_vfp(wc, current); const int table_id = well_controls_iget_vfp(wc, current);
const double target = well_controls_iget_target(wc, current); const double target = well_controls_iget_target(wc, current);
@ -1873,7 +1326,7 @@ namespace detail {
//Perform hydrostatic correction to computed targets //Perform hydrostatic correction to computed targets
double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const ADB::V dp_v = detail::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, stdWells().wellPerforationDensities(), gravity); const ADB::V dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, asImpl().stdWells().wellPerforationDensities(), gravity);
const ADB dp = ADB::constant(dp_v); const ADB dp = ADB::constant(dp_v);
const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw); 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); const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw);
@ -2231,6 +1684,9 @@ namespace detail {
} }
// 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().updateWellState(dwells,well_state);
// Update phase conditions used for property calculations. // Update phase conditions used for property calculations.
@ -2241,106 +1697,6 @@ namespace detail {
template <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::updateWellState(const V& 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 V dqs = subset(dwells, Span(np*nw, 1, varstart));
varstart += dqs.size();
const V dbhp = subset(dwells, Span(nw, 1, varstart));
varstart += dbhp.size();
assert(varstart == dwells.size());
const double dpmaxrel = dpMaxRel();
// 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 DataBlock wwr = Eigen::Map<const DataBlock>(dqs.data(), np, nw).transpose();
const V dwr = Eigen::Map<const V>(wwr.data(), nw*np);
const V wr_old = Eigen::Map<const V>(&well_state.wellRates()[0], nw*np);
const V wr = wr_old - dwr;
std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
// Bhp update.
const V bhp_old = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel);
const V bhp = bhp_old - dbhp_limited;
std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
//Get gravity for THP hydrostatic correction
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
// Thp update
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 = wr[w*np + pu.phase_pos[ Water ] ];
}
if (active_[ Oil ]) {
liquid = wr[w*np + pu.phase_pos[ Oil ] ];
}
if (active_[ Gas ]) {
vapour = wr[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 = detail::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getInj()->getTable(table_id)->getDatumDepth(),
stdWells().wellPerforationDensities(), gravity);
well_state.thp()[w] = vfp_properties_.getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp);
}
else if (well_type == PRODUCER) {
double dp = detail::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getProd()->getTable(table_id)->getDatumDepth(),
stdWells().wellPerforationDensities(), gravity);
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");
}
//Assume only one THP control specified for each well
break;
}
}
}
}
}
template <class Grid, class Implementation> template <class Grid, class Implementation>
std::vector<ADB> std::vector<ADB>
BlackoilModelBase<Grid, Implementation>::computeRelPerm(const SolutionState& state) const BlackoilModelBase<Grid, Implementation>::computeRelPerm(const SolutionState& state) const
@ -3171,6 +2527,24 @@ namespace detail {
// TODO: only kept for now due to flow_multisegment
// will be removed soon
template <class Grid, class Implementation>
void
BlackoilModelBase<Grid, 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);
}
} // namespace Opm } // namespace Opm
#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED

8
opm/autodiff/BlackoilMultiSegmentModel.hpp
View File

@ -222,6 +222,7 @@ namespace Opm {
using Base::convergenceReduction; using Base::convergenceReduction;
using Base::maxResidualAllowed; using Base::maxResidualAllowed;
using Base::variableState; using Base::variableState;
using Base::variableWellStateIndices;
using Base::asImpl; using Base::asImpl;
const std::vector<WellMultiSegmentConstPtr>& wellsMultiSegment() const { return wells_multisegment_; } const std::vector<WellMultiSegmentConstPtr>& wellsMultiSegment() const { return wells_multisegment_; }
@ -239,6 +240,13 @@ namespace Opm {
void computeWellConnectionPressures(const SolutionState& state, void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw); 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);
bool bool
solveWellEq(const std::vector<ADB>& mob_perfcells, solveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells, const std::vector<ADB>& b_perfcells,

110
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
View File

@ -964,7 +964,7 @@ namespace Opm {
// inequality constraint, and therefore skipped. // inequality constraint, and therefore skipped.
continue; continue;
} }
if (detail::constraintBroken( if (wellhelpers::constraintBroken(
xw.bhp(), xw.thp(), xw.wellRates(), xw.bhp(), xw.thp(), xw.wellRates(),
w, np, wellsMultiSegment()[w]->wellType(), wc, ctrl_index)) { w, np, wellsMultiSegment()[w]->wellType(), wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop. // ctrl_index will be the index of the broken constraint after the loop.
@ -1033,7 +1033,7 @@ namespace Opm {
SolutionState& state, SolutionState& state,
WellState& well_state) WellState& well_state)
{ {
const bool converged = Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); const bool converged = baseSolveWellEq(mob_perfcells, b_perfcells, state, well_state);
if (converged) { if (converged) {
// We must now update the state.segp and state.segqs members, // We must now update the state.segp and state.segqs members,
@ -1547,6 +1547,112 @@ namespace Opm {
} }
/// 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 = 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);
computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
addWellFluxEq(cq_s, wellSolutionState);
addWellControlEq(wellSolutionState, well_state, aliveWells);
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);
updateWellState(dx.array(), well_state);
updateWellControls(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;
}
} // namespace Opm } // namespace Opm
#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED

17
opm/autodiff/BlackoilSolventModel.hpp
View File

@ -25,6 +25,7 @@
#include <opm/autodiff/BlackoilSolventState.hpp> #include <opm/autodiff/BlackoilSolventState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp> #include <opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp> #include <opm/autodiff/SolventPropsAdFromDeck.hpp>
#include <opm/autodiff/StandardWellsSolvent.hpp>
namespace Opm { namespace Opm {
@ -103,6 +104,9 @@ namespace Opm {
const bool is_miscible_; const bool is_miscible_;
std::vector<ADB> mu_eff_; std::vector<ADB> mu_eff_;
std::vector<ADB> b_eff_; std::vector<ADB> b_eff_;
StandardWellsSolvent std_wells_;
const StandardWellsSolvent& stdWells() const { return std_wells_; }
StandardWellsSolvent& stdWells() { return std_wells_; }
// Need to declare Base members we want to use here. // Need to declare Base members we want to use here.
@ -130,7 +134,7 @@ namespace Opm {
// --------- Protected methods --------- // --------- Protected methods ---------
// Need to declare Base members we want to use here. // Need to declare Base members we want to use here.
using Base::stdWells; // using Base::stdWells;
using Base::wells; using Base::wells;
using Base::variableState; using Base::variableState;
using Base::computeGasPressure; using Base::computeGasPressure;
@ -145,10 +149,10 @@ namespace Opm {
using Base::dsMax; using Base::dsMax;
using Base::drMaxRel; using Base::drMaxRel;
using Base::maxResidualAllowed; using Base::maxResidualAllowed;
using Base::updateWellControls; // using Base::updateWellControls;
using Base::computeWellConnectionPressures; using Base::computeWellConnectionPressures;
using Base::addWellControlEq; using Base::addWellControlEq;
using Base::computePropertiesForWellConnectionPressures; // using Base::computePropertiesForWellConnectionPressures;
std::vector<ADB> std::vector<ADB>
computeRelPerm(const SolutionState& state) const; computeRelPerm(const SolutionState& state) const;
@ -202,13 +206,6 @@ namespace Opm {
const SolutionState& state, const SolutionState& state,
WellState& xw); WellState& xw);
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);
void updateEquationsScaling(); void updateEquationsScaling();
void void

129
opm/autodiff/BlackoilSolventModel_impl.hpp
View File

@ -89,7 +89,8 @@ namespace Opm {
has_solvent_(has_solvent), has_solvent_(has_solvent),
solvent_pos_(detail::solventPos(fluid.phaseUsage())), solvent_pos_(detail::solventPos(fluid.phaseUsage())),
solvent_props_(solvent_props), solvent_props_(solvent_props),
is_miscible_(is_miscible) is_miscible_(is_miscible),
std_wells_(wells_arg, solvent_props, solvent_pos_)
{ {
if (has_solvent_) { if (has_solvent_) {
@ -381,132 +382,6 @@ namespace Opm {
} }
} }
template <class Grid>
void BlackoilSolventModel<Grid>::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)
{
using namespace Opm::AutoDiffGrid;
// 1. Compute properties required by computeConnectionPressureDelta().
// Note that some of the complexity of this part is due to the function
// taking std::vector<double> arguments, and not Eigen objects.
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
// Compute the average pressure in each well block
const V perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
V avg_press = perf_press*0;
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
}
// Use cell values for the temperature as the wells don't knows its temperature yet.
const ADB perf_temp = subset(state.temperature, well_cells);
// Compute b, rsmax, rvmax values for perforations.
// Evaluate the properties using average well block pressures
// and cell values for rs, rv, phase condition and temperature.
const ADB avg_press_ad = ADB::constant(avg_press);
std::vector<PhasePresence> perf_cond(nperf);
const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf; ++perf) {
perf_cond[perf] = pc[well_cells[perf]];
}
const PhaseUsage& pu = fluid_.phaseUsage();
DataBlock b(nperf, pu.num_phases);
const V 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]);
const ADB perf_rv = subset(state.rv, well_cells);
const ADB perf_rs = subset(state.rs, well_cells);
const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const V 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 V rssat = fluidRsSat(avg_press, perf_so, well_cells);
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);
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);
}
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();
V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
V rhog = fluid_.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
if (has_solvent_) {
const V bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value();
//const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value();
// A weighted sum of the b-factors of gas and solvent are used.
const int nc = Opm::AutoDiffGrid::numCells(grid_);
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);
V F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value();
V injectedSolventFraction = Eigen::Map<const V>(&xw.solventFraction()[0], nperf);
V isProducer = V::Zero(nperf);
V ones = V::Constant(nperf,1.0);
for (int w = 0; w < nw; ++w) {
if(wells().type[w] == PRODUCER) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
isProducer[perf] = 1;
}
}
}
F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction;
bg = bg * (ones - F_solvent);
bg = bg + F_solvent * bs;
const V& rhos = solvent_props_.solventSurfaceDensity(well_cells);
rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
}
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases);
const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
} else {
rvmax_perf.assign(0.0, nperf);
}
// b and surf_dens_perf is row major, so can just copy data.
b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases);
}

159
opm/autodiff/StandardWells.cpp
View File

@ -1,159 +0,0 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/autodiff/StandardWells.hpp>
namespace Opm
{
StandardWells::
WellOps::WellOps(const Wells* wells)
: w2p(),
p2w(),
well_cells()
{
if( wells )
{
w2p = Eigen::SparseMatrix<double>(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells);
p2w = Eigen::SparseMatrix<double>(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]);
const int nw = wells->number_of_wells;
const int* const wpos = wells->well_connpos;
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> scatter, gather;
scatter.reserve(wpos[nw]);
gather .reserve(wpos[nw]);
for (int w = 0, i = 0; w < nw; ++w) {
for (; i < wpos[ w + 1 ]; ++i) {
scatter.push_back(Tri(i, w, 1.0));
gather .push_back(Tri(w, i, 1.0));
}
}
w2p.setFromTriplets(scatter.begin(), scatter.end());
p2w.setFromTriplets(gather .begin(), gather .end());
well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]);
}
}
StandardWells::StandardWells(const Wells* wells_arg)
: wells_(wells_arg)
, wops_(wells_arg)
, well_perforation_densities_(Vector())
, well_perforation_pressure_diffs_(Vector())
{
}
const Wells& StandardWells::wells() const
{
assert(wells_ != 0);
return *(wells_);
}
bool StandardWells::wellsActive() const
{
return wells_active_;
}
void StandardWells::setWellsActive(const bool wells_active)
{
wells_active_ = wells_active;
}
bool StandardWells::localWellsActive() const
{
return wells_ ? (wells_->number_of_wells > 0 ) : false;
}
const StandardWells::WellOps&
StandardWells::wellOps() const
{
return wops_;
}
Vector& StandardWells::wellPerforationDensities()
{
return well_perforation_densities_;
}
const Vector& StandardWells::wellPerforationDensities() const
{
return well_perforation_densities_;
}
Vector& StandardWells::wellPerforationPressureDiffs()
{
return well_perforation_pressure_diffs_;
}
const Vector& StandardWells::wellPerforationPressureDiffs() const
{
return well_perforation_pressure_diffs_;
}
}

84
opm/autodiff/StandardWells.hpp
View File

@ -31,13 +31,12 @@
#include <opm/core/wells.h> #include <opm/core/wells.h>
#include <opm/autodiff/AutoDiffBlock.hpp> #include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
namespace Opm { namespace Opm {
// --------- Types ---------
typedef AutoDiffBlock<double> ADB;
typedef ADB::V Vector;
/// Class for handling the standard well model. /// Class for handling the standard well model.
class StandardWells { class StandardWells {
@ -50,6 +49,16 @@ namespace Opm {
}; };
public: public:
// --------- Types ---------
using ADB = AutoDiffBlock<double>;
using Vector = ADB::V;
// copied from BlackoilModelBase
// should put to somewhere better
using DataBlock = Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor>;
// --------- Public methods --------- // --------- Public methods ---------
explicit StandardWells(const Wells* wells); explicit StandardWells(const Wells* wells);
@ -64,13 +73,75 @@ namespace Opm {
const WellOps& wellOps() const; const WellOps& wellOps() const;
/// Density of each well perforation /// Density of each well perforation
Vector& wellPerforationDensities(); Vector& wellPerforationDensities(); // mutable version kept for BlackoilMultisegmentModel
const Vector& wellPerforationDensities() const; const Vector& wellPerforationDensities() const;
/// Diff to bhp for each well perforation. /// Diff to bhp for each well perforation.
Vector& wellPerforationPressureDiffs(); Vector& wellPerforationPressureDiffs(); // mutable version kept for BlackoilMultisegmentModel
const Vector& wellPerforationPressureDiffs() const; 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,
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;
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,
WellState& xw) const;
protected: protected:
bool wells_active_; bool wells_active_;
const Wells* wells_; const Wells* wells_;
@ -81,4 +152,7 @@ namespace Opm {
} // namespace Opm } // namespace Opm
#include "StandardWells_impl.hpp"
#endif #endif

72
opm/autodiff/StandardWellsSolvent.hpp Normal file
View File

@ -0,0 +1,72 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_STANDARDWELLSSOLVENT_HEADER_INCLUDED
#define OPM_STANDARDWELLSSOLVENT_HEADER_INCLUDED
#include <opm/autodiff/StandardWells.hpp>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
namespace Opm {
/// Class for handling the standard well model for solvent model
class StandardWellsSolvent : public StandardWells
{
public:
using Base = StandardWells;
// --------- Public methods ---------
explicit StandardWellsSolvent(const Wells* wells, const SolventPropsAdFromDeck& solvent_props, const int solvent_pos);
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);
// TODO: fluid and active may be can put in the member list
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;
protected:
const SolventPropsAdFromDeck& solvent_props_;
const int solvent_pos_;
};
} // namespace Opm
#include "StandardWellsSolvent_impl.hpp"
#endif

231
opm/autodiff/StandardWellsSolvent_impl.hpp Normal file
View File

@ -0,0 +1,231 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/autodiff/StandardWellsSolvent.hpp>
namespace Opm
{
StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg,
const SolventPropsAdFromDeck& solvent_props,
const int solvent_pos)
: Base(wells_arg)
, solvent_props_(solvent_props)
, solvent_pos_(solvent_pos)
{
}
template<class SolutionState, class WellState>
void
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,
std::vector<double>& surf_dens_perf)
{
// 1. Compute properties required by computeConnectionPressureDelta().
// Note that some of the complexity of this part is due to the function
// taking std::vector<double> arguments, and not Eigen objects.
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
// Compute the average pressure in each well block
const Vector perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
Vector avg_press = perf_press*0;
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
}
const std::vector<int>& well_cells = wellOps().well_cells;
// 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);
for (int perf = 0; perf < nperf; ++perf) {
perf_cond[perf] = pc[well_cells[perf]];
}
const PhaseUsage& pu = fluid.phaseUsage();
DataBlock b(nperf, pu.num_phases);
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]);
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 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();
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);
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);
}
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);
// to handle solvent related
{
const Vector bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value();
//const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value();
// number of cells
const int nc = state.pressure.size();
const ADB zero = ADB::constant(Vector::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);
Vector F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value();
Vector injectedSolventFraction = Eigen::Map<const Vector>(&xw.solventFraction()[0], nperf);
Vector isProducer = Vector::Zero(nperf);
Vector ones = Vector::Constant(nperf,1.0);
for (int w = 0; w < nw; ++w) {
if(wells().type[w] == PRODUCER) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
isProducer[perf] = 1;
}
}
}
F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction;
bg = bg * (ones - F_solvent);
bg = bg + F_solvent * bs;
const Vector& rhos = solvent_props_.solventSurfaceDensity(well_cells);
rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
}
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
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();
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
} else {
rvmax_perf.assign(0.0, nperf);
}
// b and surf_dens_perf is row major, so can just copy data.
b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases);
}
template <class ReservoirResidualQuant, class SolutionState>
void
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);
// handle the solvent related
{
int gas_pos = fluid.phaseUsage().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
// 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 = 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 ]
? 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);
Vector ones = Vector::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));
}
}
}

763
opm/autodiff/StandardWells_impl.hpp Normal file
View File

@ -0,0 +1,763 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/autodiff/StandardWells.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/WellHelpers.hpp>
namespace Opm
{
StandardWells::
WellOps::WellOps(const Wells* wells)
: w2p(),
p2w(),
well_cells()
{
if( wells )
{
w2p = Eigen::SparseMatrix<double>(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells);
p2w = Eigen::SparseMatrix<double>(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]);
const int nw = wells->number_of_wells;
const int* const wpos = wells->well_connpos;
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> scatter, gather;
scatter.reserve(wpos[nw]);
gather .reserve(wpos[nw]);
for (int w = 0, i = 0; w < nw; ++w) {
for (; i < wpos[ w + 1 ]; ++i) {
scatter.push_back(Tri(i, w, 1.0));
gather .push_back(Tri(w, i, 1.0));
}
}
w2p.setFromTriplets(scatter.begin(), scatter.end());
p2w.setFromTriplets(gather .begin(), gather .end());
well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]);
}
}
StandardWells::StandardWells(const Wells* wells_arg)
: wells_(wells_arg)
, wops_(wells_arg)
, well_perforation_densities_(Vector())
, well_perforation_pressure_diffs_(Vector())
{
}
const Wells& StandardWells::wells() const
{
assert(wells_ != 0);
return *(wells_);
}
bool StandardWells::wellsActive() const
{
return wells_active_;
}
void StandardWells::setWellsActive(const bool wells_active)
{
wells_active_ = wells_active;
}
bool StandardWells::localWellsActive() const
{
return wells_ ? (wells_->number_of_wells > 0 ) : false;
}
const StandardWells::WellOps&
StandardWells::wellOps() const
{
return wops_;
}
StandardWells::Vector& StandardWells::wellPerforationDensities()
{
return well_perforation_densities_;
}
const StandardWells::Vector&
StandardWells::wellPerforationDensities() const
{
return well_perforation_densities_;
}
StandardWells::Vector&
StandardWells::wellPerforationPressureDiffs()
{
return well_perforation_pressure_diffs_;
}
const StandardWells::Vector&
StandardWells::wellPerforationPressureDiffs() const
{
return well_perforation_pressure_diffs_;
}
template<class SolutionState, class WellState>
void
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)
{
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
// Compute the average pressure in each well block
const Vector perf_press = Eigen::Map<const Vector>(xw.perfPress().data(), nperf);
Vector avg_press = perf_press*0;
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
}
const std::vector<int>& well_cells = wellOps().well_cells;
// Use cell values for the temperature as the wells don't knows its temperature yet.
const ADB perf_temp = subset(state.temperature, well_cells);
// Compute b, rsmax, rvmax values for perforations.
// Evaluate the properties using average well block pressures
// and cell values for rs, rv, phase condition and temperature.
const ADB avg_press_ad = ADB::constant(avg_press);
std::vector<PhasePresence> perf_cond(nperf);
// const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf; ++perf) {
perf_cond[perf] = pc[well_cells[perf]];
}
const PhaseUsage& pu = fluid.phaseUsage();
DataBlock b(nperf, pu.num_phases);
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 V rssat = fluidRsSat(avg_press, perf_so, well_cells);
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(nperf, 0.0);
}
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 V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
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(nperf, 0.0);
}
// b is row major, so can just copy data.
b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
// 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);
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);
}
surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases);
}
template <class WellState>
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,
const std::vector<double>& surf_dens_perf,
const std::vector<double>& depth_perf,
const double grav)
{
// Compute densities
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid.phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
const int nperf = wells().well_connpos[wells().number_of_wells];
// Compute pressure deltas
std::vector<double> cdp =
WellDensitySegmented::computeConnectionPressureDelta(
wells(), depth_perf, cd, 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);
}
template <class ReservoirResidualQuant, class SolutionState>
void
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
{
// 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(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);
}
}
}
template <class SolutionState>
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,
std::vector<ADB>& cq_s) const
{
if( ! localWellsActive() ) return ;
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
Vector Tw = Eigen::Map<const Vector>(wells().WI, nperf);
const std::vector<int>& well_cells = wellOps().well_cells;
// pressure diffs computed already (once per step, not changing per iteration)
const Vector& cdp = wellPerforationPressureDiffs();
// Extract needed quantities for the perforation cells
const ADB& p_perfcells = subset(state.pressure, well_cells);
const ADB& rv_perfcells = subset(state.rv, well_cells);
const ADB& rs_perfcells = subset(state.rs, well_cells);
// Perforation pressure
const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp;
// Pressure drawdown (also used to determine direction of flow)
const ADB drawdown = p_perfcells - perfpressure;
// Compute vectors with zero and ones that
// selects the wanted quantities.
// selects injection perforations
Vector selectInjectingPerforations = Vector::Zero(nperf);
// selects producing perforations
Vector selectProducingPerforations = Vector::Zero(nperf);
for (int c = 0; c < nperf; ++c){
if (drawdown.value()[c] < 0)
selectInjectingPerforations[c] = 1;
else
selectProducingPerforations[c] = 1;
}
// Handle cross flow
const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value();
const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value();
for (int w = 0; w < nw; ++w) {
if (!wells().allow_cf[w]) {
for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) {
// Crossflow is not allowed; reverse flow is prevented.
// At least one of the perforation must be open in order to have a meeningful
// equation to solve. For the special case where all perforations have reverse flow,
// and the target rate is non-zero all of the perforations are keept open.
if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
selectProducingPerforations[perf] = 0.0;
} else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
selectInjectingPerforations[perf] = 0.0;
}
}
}
}
// HANDLE FLOW INTO WELLBORE
// compute phase volumetric rates at standard conditions
std::vector<ADB> cq_ps(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
cq_ps[phase] = b_perfcells[phase] * cq_p;
}
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];
const ADB cq_psGas = cq_ps[gaspos];
cq_ps[gaspos] += rs_perfcells * cq_psOil;
cq_ps[oilpos] += rv_perfcells * cq_psGas;
}
// HANDLE FLOW OUT FROM WELLBORE
// Using total mobilities
ADB total_mob = mob_perfcells[0];
for (int phase = 1; phase < np; ++phase) {
total_mob += mob_perfcells[phase];
}
// injection perforations total volume rates
const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
// compute wellbore mixture for injecting perforations
// The wellbore mixture depends on the inflow from the reservoar
// and the well injection rates.
// compute avg. and total wellbore phase volumetric rates at standard conds
const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
std::vector<ADB> wbq(np, ADB::null());
ADB wbqt = ADB::constant(Vector::Zero(nw));
for (int phase = 0; phase < np; ++phase) {
const ADB& q_ps = wellOps().p2w * cq_ps[phase];
const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw));
Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
const int pos = pu.phase_pos[phase];
wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps;
wbqt += wbq[phase];
}
// compute wellbore mixture at standard conditions.
Selector<double> notDeadWells_selector(wbqt.value(), Selector<double>::Zero);
std::vector<ADB> cmix_s(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
const int pos = pu.phase_pos[phase];
cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
}
// compute volume ratio between connection at standard conditions
ADB volumeRatio = ADB::constant(Vector::Zero(nperf));
const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
for (int phase = 0; phase < np; ++phase) {
ADB tmp = cmix_s[phase];
if (phase == Oil && active[Gas]) {
const int gaspos = pu.phase_pos[Gas];
tmp -= rv_perfcells * cmix_s[gaspos] / d;
}
if (phase == Gas && active[Oil]) {
const int oilpos = pu.phase_pos[Oil];
tmp -= rs_perfcells * cmix_s[oilpos] / d;
}
volumeRatio += tmp / b_perfcells[phase];
}
// injecting connections total volumerates at standard conditions
ADB cqt_is = cqt_i/volumeRatio;
// connection phase volumerates at standard conditions
cq_s.resize(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
}
// check for dead wells (used in the well controll equations)
aliveWells = Vector::Constant(nw, 1.0);
for (int w = 0; w < nw; ++w) {
if (wbqt.value()[w] == 0) {
aliveWells[w] = 0.0;
}
}
}
template <class SolutionState, class WellState>
void
StandardWells::
updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw) const
{
if ( !localWellsActive() )
{
// If there are no wells in the subdomain of the proces then
// cq_s has zero size and will cause a segmentation fault below.
return;
}
// Update the perforation phase rates (used to calculate the pressure drop in the wellbore).
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
Vector cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np);
}
xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np);
// Update the perforation pressures.
const Vector& cdp = wellPerforationPressureDiffs();
const Vector perfpressure = (wellOps().w2p * state.bhp.value().matrix()).array() + cdp;
xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf);
}
template <class WellState>
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 nw = wells().number_of_wells;
// Extract parts of dwells corresponding to each part.
int varstart = 0;
const Vector dqs = subset(dwells, Span(np*nw, 1, varstart));
varstart += dqs.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 DataBlock wwr = Eigen::Map<const DataBlock>(dqs.data(), np, nw).transpose();
const Vector dwr = Eigen::Map<const Vector>(wwr.data(), nw*np);
const Vector wr_old = Eigen::Map<const Vector>(&well_state.wellRates()[0], nw*np);
const Vector wr = wr_old - dwr;
std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
// Bhp update.
const Vector bhp_old = Eigen::Map<const Vector>(&well_state.bhp()[0], nw, 1);
const Vector dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel);
const Vector bhp = bhp_old - dbhp_limited;
std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
//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 = wr[w*np + pu.phase_pos[ Water ] ];
}
if (active[ Oil ]) {
liquid = wr[w*np + pu.phase_pos[ Oil ] ];
}
if (active[ Gas ]) {
vapour = wr[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, 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, 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>
void
StandardWells::
updateWellControls(const Opm::PhaseUsage& pu,
const double gravity,
const VFPProperties& vfp_properties,
const bool terminal_output,
const std::vector<bool>& active,
WellState& xw) const
{
if( !localWellsActive() ) 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 = wells().number_of_phases;
const int nw = wells().number_of_wells;
#pragma omp parallel for schedule(dynamic)
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
int current = xw.currentControls()[w];
// Loop over all controls except the current one, and also
// skip any RESERVOIR_RATE controls, since we cannot
// handle those.
const int nwc = well_controls_get_num(wc);
int ctrl_index = 0;
for (; ctrl_index < nwc; ++ctrl_index) {
if (ctrl_index == current) {
// This is the currently used control, so it is
// used as an equation. So this is not used as an
// inequality constraint, and therefore skipped.
continue;
}
if (wellhelpers::constraintBroken(
xw.bhp(), xw.thp(), xw.wellRates(),
w, np, wells().type[w], wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
}
}
if (ctrl_index != nwc) {
// Constraint number ctrl_index was broken, switch to it.
if (terminal_output)
{
std::cout << "Switching control mode for well " << wells().name[w]
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
}
xw.currentControls()[w] = ctrl_index;
current = xw.currentControls()[w];
}
// Updating well state and primary variables.
// Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
const double target = well_controls_iget_target(wc, current);
const double* distr = well_controls_iget_distr(wc, current);
switch (well_controls_iget_type(wc, current)) {
case BHP:
xw.bhp()[w] = target;
break;
case THP: {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
if (active[ Water ]) {
aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if (active[ Oil ]) {
liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if (active[ Gas ]) {
vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ];
}
const int vfp = well_controls_iget_vfp(wc, current);
const double& thp = well_controls_iget_target(wc, current);
const double& alq = well_controls_iget_alq(wc, current);
//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);
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);
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");
}
break;
}
case RESERVOIR_RATE:
// No direct change to any observable quantity at
// surface condition. In this case, use existing
// flow rates as initial conditions as reservoir
// rate acts only in aggregate.
break;
case SURFACE_RATE:
// assign target value as initial guess for injectors and
// single phase producers (orat, grat, wrat)
const WellType& well_type = wells().type[w];
if (well_type == INJECTOR) {
for (int phase = 0; phase < np; ++phase) {
const double& compi = wells().comp_frac[np * w + phase];
if (compi > 0.0) {
xw.wellRates()[np*w + phase] = target * compi;
}
}
} else if (well_type == PRODUCER) {
// only set target as initial rates for single phase
// producers. (orat, grat and wrat, and not lrat)
// lrat will result in numPhasesWithTargetsUnderThisControl == 2
int numPhasesWithTargetsUnderThisControl = 0;
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0) {
numPhasesWithTargetsUnderThisControl += 1;
}
}
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) {
xw.wellRates()[np*w + phase] = target * distr[phase];
}
}
} else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
break;
}
}
}
}

170
opm/autodiff/WellHelpers.hpp Normal file
View File

@ -0,0 +1,170 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLHELPERS_HEADER_INCLUDED
#define OPM_WELLHELPERS_HEADER_INCLUDED
#include <opm/core/wells.h>
#include <opm/autodiff/AutoDiffBlock.hpp>
// #include <opm/autodiff/AutoDiffHelpers.hpp>
#include <vector>
namespace Opm {
namespace wellhelpers
{
inline
double rateToCompare(const std::vector<double>& well_phase_flow_rate,
const int well,
const int num_phases,
const double* distr)
{
double rate = 0.0;
for (int phase = 0; phase < num_phases; ++phase) {
// Important: well_phase_flow_rate is ordered with all phase rates for first
// well first, then all phase rates for second well etc.
rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase];
}
return rate;
}
inline
bool constraintBroken(const std::vector<double>& bhp,
const std::vector<double>& thp,
const std::vector<double>& well_phase_flow_rate,
const int well,
const int num_phases,
const WellType& well_type,
const WellControls* wc,
const int ctrl_index)
{
const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
const double target = well_controls_iget_target(wc, ctrl_index);
const double* distr = well_controls_iget_distr(wc, ctrl_index);
bool broken = false;
switch (well_type) {
case INJECTOR:
{
switch (ctrl_type) {
case BHP:
broken = bhp[well] > target;
break;
case THP:
broken = thp[well] > target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
broken = rateToCompare(well_phase_flow_rate,
well, num_phases, distr) > target;
break;
}
}
break;
case PRODUCER:
{
switch (ctrl_type) {
case BHP:
broken = bhp[well] < target;
break;
case THP:
broken = thp[well] < target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
// Note that the rates compared below are negative,
// so breaking the constraints means: too high flow rate
// (as for injection).
broken = rateToCompare(well_phase_flow_rate,
well, num_phases, distr) < target;
break;
}
}
break;
default:
OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells.");
}
return broken;
}
// --------- Types ---------
using Vector = AutoDiffBlock<double>::V;
/**
* Simple hydrostatic correction for VFP table
* @param wells - wells struct
* @param w Well number
* @param vfp_table VFP table
* @param well_perforation_densities Densities at well perforations
* @param gravity Gravitational constant (e.g., 9.81...)
*/
inline
double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth,
const Vector& well_perforation_densities, const double gravity) {
if ( wells.well_connpos[w] == wells.well_connpos[w+1] )
{
// This is a well with no perforations.
// If this is the last well we would subscript over the
// bounds below.
// we assume well_perforation_densities to be 0
return 0;
}
const double well_ref_depth = wells.depth_ref[w];
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;
return dp;
}
inline
Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth,
const Vector& well_perforation_densities, const double gravity) {
const int nw = wells.number_of_wells;
Vector retval = Vector::Zero(nw);
#pragma omp parallel for schedule(static)
for (int i=0; i<nw; ++i) {
retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities, gravity);
}
return retval;
}
} // namespace wellhelpers
}
#endif

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

@ -174,6 +174,7 @@ namespace Opm {
using Base::terminal_output_; using Base::terminal_output_;
using Base::primalVariable_; using Base::primalVariable_;
using Base::pvdt_; using Base::pvdt_;
using Base::vfp_properties_;
// --------- Protected methods --------- // --------- Protected methods ---------
@ -199,7 +200,7 @@ namespace Opm {
using Base::drMaxRel; using Base::drMaxRel;
using Base::maxResidualAllowed; using Base::maxResidualAllowed;
using Base::updateWellControls; // using Base::updateWellControls;
using Base::computeWellConnectionPressures; using Base::computeWellConnectionPressures;
using Base::addWellControlEq; using Base::addWellControlEq;
using Base::computeRelPerm; using Base::computeRelPerm;

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

@ -498,7 +498,9 @@ namespace Opm {
// Possibly switch well controls and updating well state to // Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells // get reasonable initial conditions for the wells
updateWellControls(well_state); const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
// updateWellControls(well_state);
stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state);
// Create the primary variables. // Create the primary variables.
SolutionState state = variableState(reservoir_state, well_state); SolutionState state = variableState(reservoir_state, well_state);
@ -550,7 +552,7 @@ namespace Opm {
Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
} }
Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s);
if (has_plyshlog_) { if (has_plyshlog_) {
std::vector<double> water_vel_wells; std::vector<double> water_vel_wells;
@ -569,8 +571,8 @@ namespace Opm {
mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb; mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb;
} }
Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s);
Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state); stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
Base::addWellFluxEq(cq_s, state); Base::addWellFluxEq(cq_s, state);
addWellContributionToMassBalanceEq(cq_s, state, well_state); addWellContributionToMassBalanceEq(cq_s, state, well_state);
addWellControlEq(state, well_state, aliveWells); addWellControlEq(state, well_state, aliveWells);