OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-30 11:06:55 -06:00
Code Issues Packages Projects Releases Wiki Activity
bc8da79779
opm-simulators/opm/simulators/wells/MultisegmentWell_impl.hpp

1954 lines
84 KiB
C++
Raw Blame History

/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 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/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/MSW/Valve.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <string>
#include <algorithm>
#if HAVE_CUDA || HAVE_OPENCL
#include <opm/simulators/linalg/bda/WellContributions.hpp>
#endif
namespace Opm
{
template <typename TypeTag>
MultisegmentWell<TypeTag>::
MultisegmentWell(const Well& well,
const ParallelWellInfo& pw_info,
const int time_step,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData>& perf_data)
: Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
, MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
, segment_fluid_initial_(this->numberOfSegments(), std::vector<double>(this->num_components_, 0.0))
{
// not handling solvent or polymer for now with multisegment well
if constexpr (has_solvent) {
OPM_THROW(std::runtime_error, "solvent is not supported by multisegment well yet");
}
if constexpr (has_polymer) {
OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet");
}
if constexpr (Base::has_energy) {
OPM_THROW(std::runtime_error, "energy is not supported by multisegment well yet");
}
if constexpr (Base::has_foam) {
OPM_THROW(std::runtime_error, "foam is not supported by multisegment well yet");
}
if constexpr (Base::has_brine) {
OPM_THROW(std::runtime_error, "brine is not supported by multisegment well yet");
}
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg,
const double gravity_arg,
const int num_cells,
const std::vector< Scalar >& B_avg)
{
Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg);
// TODO: for StandardWell, we need to update the perf depth here using depth_arg.
// for MultisegmentWell, it is much more complicated.
// It can be specified directly, it can be calculated from the segment depth,
// it can also use the cell center, which is the same for StandardWell.
// For the last case, should we update the depth with the depth_arg? For the
// future, it can be a source of wrong result with Multisegment well.
// An indicator from the opm-parser should indicate what kind of depth we should use here.
// \Note: we do not update the depth here. And it looks like for now, we only have the option to use
// specified perforation depth
this->initMatrixAndVectors(num_cells);
// calcuate the depth difference between the perforations and the perforated grid block
for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
const int cell_idx = this->well_cells_[perf];
this->cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - this->perf_depth_[perf];
}
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
initPrimaryVariablesEvaluation() const
{
this->MSWEval::initPrimaryVariablesEvaluation();
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
updatePrimaryVariables(const WellState& well_state, DeferredLogger& /* deferred_logger */) const
{
this->MSWEval::updatePrimaryVariables(well_state);
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
updateWellStateWithTarget(const Simulator& ebos_simulator,
const GroupState& group_state,
WellState& well_state,
DeferredLogger& deferred_logger) const
{
Base::updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
// scale segment rates based on the wellRates
// and segment pressure based on bhp
this->scaleSegmentRatesWithWellRates(well_state);
this->scaleSegmentPressuresWithBhp(well_state);
}
template <typename TypeTag>
ConvergenceReport
MultisegmentWell<TypeTag>::
getWellConvergence(const WellState& well_state,
const std::vector<double>& B_avg,
DeferredLogger& deferred_logger,
const bool relax_tolerance) const
{
return this->MSWEval::getWellConvergence(well_state,
B_avg,
deferred_logger,
this->param_.max_residual_allowed_,
this->param_.tolerance_wells_,
this->param_.relaxed_tolerance_flow_well_,
this->param_.tolerance_pressure_ms_wells_,
this->param_.relaxed_tolerance_pressure_ms_well_,
relax_tolerance);
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
apply(const BVector& x, BVector& Ax) const
{
if (!this->isOperable() && !this->wellIsStopped()) return;
if ( this->param_.matrix_add_well_contributions_ )
{
// Contributions are already in the matrix itself
return;
}
BVectorWell Bx(this->duneB_.N());
this->duneB_.mv(x, Bx);
// invDBx = duneD^-1 * Bx_
const BVectorWell invDBx = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, Bx);
// Ax = Ax - duneC_^T * invDBx
this->duneC_.mmtv(invDBx,Ax);
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
apply(BVector& r) const
{
if (!this->isOperable() && !this->wellIsStopped()) return;
// invDrw_ = duneD^-1 * resWell_
const BVectorWell invDrw = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
// r = r - duneC_^T * invDrw
this->duneC_.mmtv(invDrw, r);
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
recoverWellSolutionAndUpdateWellState(const BVector& x,
WellState& well_state,
DeferredLogger& deferred_logger) const
{
if (!this->isOperable() && !this->wellIsStopped()) return;
BVectorWell xw(1);
this->recoverSolutionWell(x, xw);
updateWellState(xw, well_state, deferred_logger);
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials,
DeferredLogger& deferred_logger)
{
const int np = this->number_of_phases_;
well_potentials.resize(np, 0.0);
// Stopped wells have zero potential.
if (this->wellIsStopped()) {
return;
}
// If the well is pressure controlled the potential equals the rate.
bool thp_controlled_well = false;
bool bhp_controlled_well = false;
const auto& ws = well_state.well(this->index_of_well_);
if (this->isInjector()) {
const Well::InjectorCMode& current = ws.injection_cmode;
if (current == Well::InjectorCMode::THP) {
thp_controlled_well = true;
}
if (current == Well::InjectorCMode::BHP) {
bhp_controlled_well = true;
}
} else {
const Well::ProducerCMode& current = ws.production_cmode;
if (current == Well::ProducerCMode::THP) {
thp_controlled_well = true;
}
if (current == Well::ProducerCMode::BHP) {
bhp_controlled_well = true;
}
}
if (thp_controlled_well || bhp_controlled_well) {
double total_rate = 0.0;
for (int phase = 0; phase < np; ++phase){
total_rate += ws.surface_rates[phase];
}
// for pressure controlled wells the well rates are the potentials
// if the rates are trivial we are most probably looking at the newly
// opened well and we therefore make the affort of computing the potentials anyway.
if (std::abs(total_rate) > 0) {
for (int phase = 0; phase < np; ++phase){
well_potentials[phase] = ws.surface_rates[phase];
}
return;
}
}
debug_cost_counter_ = 0;
// does the well have a THP related constraint?
const auto& summaryState = ebosSimulator.vanguard().summaryState();
if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger);
} else {
well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger);
}
deferred_logger.debug("Cost in iterations of finding well potential for well "
+ this->name() + ": " + std::to_string(debug_cost_counter_));
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeWellRatesAtBhpLimit(const Simulator& ebosSimulator,
std::vector<double>& well_flux,
DeferredLogger& deferred_logger) const
{
if (this->well_ecl_.isInjector()) {
const auto controls = this->well_ecl_.injectionControls(ebosSimulator.vanguard().summaryState());
computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
} else {
const auto controls = this->well_ecl_.productionControls(ebosSimulator.vanguard().summaryState());
computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
}
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeWellRatesWithBhp(const Simulator& ebosSimulator,
const Scalar& bhp,
std::vector<double>& well_flux,
DeferredLogger& deferred_logger) const
{
const int np = this->number_of_phases_;
well_flux.resize(np, 0.0);
const bool allow_cf = this->getAllowCrossFlow();
const int nseg = this->numberOfSegments();
const WellState& well_state = ebosSimulator.problem().wellModel().wellState();
const auto& ws = well_state.well(this->indexOfWell());
auto segments_copy = ws.segments;
segments_copy.scale_pressure(bhp);
const auto& segment_pressure = segments_copy.pressure;
for (int seg = 0; seg < nseg; ++seg) {
for (const int perf : this->segment_perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
// flux for each perforation
std::vector<Scalar> mob(this->num_components_, 0.);
getMobilityScalar(ebosSimulator, perf, mob);
double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
const double Tw = this->well_index_[perf] * trans_mult;
const Scalar seg_pressure = segment_pressure[seg];
std::vector<Scalar> cq_s(this->num_components_, 0.);
computePerfRateScalar(intQuants, mob, Tw, seg, perf, seg_pressure,
allow_cf, cq_s, deferred_logger);
for(int p = 0; p < np; ++p) {
well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
}
}
}
this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeWellRatesWithBhpIterations(const Simulator& ebosSimulator,
const Scalar& bhp,
std::vector<double>& well_flux,
DeferredLogger& deferred_logger) const
{
// creating a copy of the well itself, to avoid messing up the explicit informations
// during this copy, the only information not copied properly is the well controls
MultisegmentWell<TypeTag> well_copy(*this);
well_copy.debug_cost_counter_ = 0;
// store a copy of the well state, we don't want to update the real well state
WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
const auto& group_state = ebosSimulator.problem().wellModel().groupState();
auto& ws = well_state_copy.well(this->index_of_well_);
// Get the current controls.
const auto& summary_state = ebosSimulator.vanguard().summaryState();
auto inj_controls = well_copy.well_ecl_.isInjector()
? well_copy.well_ecl_.injectionControls(summary_state)
: Well::InjectionControls(0);
auto prod_controls = well_copy.well_ecl_.isProducer()
? well_copy.well_ecl_.productionControls(summary_state) :
Well::ProductionControls(0);
// Set current control to bhp, and bhp value in state, modify bhp limit in control object.
if (well_copy.well_ecl_.isInjector()) {
inj_controls.bhp_limit = bhp;
ws.injection_cmode = Well::InjectorCMode::BHP;
} else {
prod_controls.bhp_limit = bhp;
ws.production_cmode = Well::ProducerCMode::BHP;
}
ws.bhp = bhp;
well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
// initialized the well rates with the potentials i.e. the well rates based on bhp
const int np = this->number_of_phases_;
const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
for (int phase = 0; phase < np; ++phase){
ws.surface_rates[phase] = sign * ws.well_potentials[phase];
}
well_copy.scaleSegmentRatesWithWellRates(well_state_copy);
well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
const double dt = ebosSimulator.timeStepSize();
// iterate to get a solution at the given bhp.
well_copy.iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
deferred_logger);
// compute the potential and store in the flux vector.
well_flux.clear();
well_flux.resize(np, 0.0);
for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
const EvalWell rate = well_copy.getQs(compIdx);
well_flux[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
}
debug_cost_counter_ += well_copy.debug_cost_counter_;
}
template<typename TypeTag>
std::vector<double>
MultisegmentWell<TypeTag>::
computeWellPotentialWithTHP(const Simulator& ebos_simulator,
DeferredLogger& deferred_logger) const
{
std::vector<double> potentials(this->number_of_phases_, 0.0);
const auto& summary_state = ebos_simulator.vanguard().summaryState();
const auto& well = this->well_ecl_;
if (well.isInjector()){
auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
if (bhp_at_thp_limit) {
const auto& controls = well.injectionControls(summary_state);
const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
deferred_logger.debug("Converged thp based potential calculation for well "
+ this->name() + ", at bhp = " + std::to_string(bhp));
} else {
deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
"Failed in getting converged thp based potential calculation for well "
+ this->name() + ". Instead the bhp based value is used");
const auto& controls = well.injectionControls(summary_state);
const double bhp = controls.bhp_limit;
computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
}
} else {
auto bhp_at_thp_limit = computeBhpAtThpLimitProd(ebos_simulator, summary_state, deferred_logger);
if (bhp_at_thp_limit) {
const auto& controls = well.productionControls(summary_state);
const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
deferred_logger.debug("Converged thp based potential calculation for well "
+ this->name() + ", at bhp = " + std::to_string(bhp));
} else {
deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
"Failed in getting converged thp based potential calculation for well "
+ this->name() + ". Instead the bhp based value is used");
const auto& controls = well.productionControls(summary_state);
const double bhp = controls.bhp_limit;
computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
}
}
return potentials;
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
{
if (!this->isOperable() && !this->wellIsStopped()) return;
// We assemble the well equations, then we check the convergence,
// which is why we do not put the assembleWellEq here.
const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
updateWellState(dx_well, well_state, deferred_logger);
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
computePerfCellPressDiffs(const Simulator& ebosSimulator)
{
for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
std::vector<double> kr(this->number_of_phases_, 0.0);
std::vector<double> density(this->number_of_phases_, 0.0);
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = intQuants.fluidState();
double sum_kr = 0.;
const PhaseUsage& pu = this->phaseUsage();
if (pu.phase_used[Water]) {
const int water_pos = pu.phase_pos[Water];
kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
sum_kr += kr[water_pos];
density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
}
if (pu.phase_used[Oil]) {
const int oil_pos = pu.phase_pos[Oil];
kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
sum_kr += kr[oil_pos];
density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
}
if (pu.phase_used[Gas]) {
const int gas_pos = pu.phase_pos[Gas];
kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
sum_kr += kr[gas_pos];
density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
}
assert(sum_kr != 0.);
// calculate the average density
double average_density = 0.;
for (int p = 0; p < this->number_of_phases_; ++p) {
average_density += kr[p] * density[p];
}
average_density /= sum_kr;
this->cell_perforation_pressure_diffs_[perf] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[perf];
}
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeInitialSegmentFluids(const Simulator& ebos_simulator)
{
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
// TODO: trying to reduce the times for the surfaceVolumeFraction calculation
const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value();
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
segment_fluid_initial_[seg][comp_idx] = surface_volume * this->surfaceVolumeFraction(seg, comp_idx).value();
}
}
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
updateWellState(const BVectorWell& dwells,
WellState& well_state,
DeferredLogger& deferred_logger,
const double relaxation_factor) const
{
if (!this->isOperable() && !this->wellIsStopped()) return;
const double dFLimit = this->param_.dwell_fraction_max_;
const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
this->MSWEval::updateWellState(dwells,
relaxation_factor,
dFLimit,
max_pressure_change);
this->updateWellStateFromPrimaryVariables(well_state, getRefDensity(), deferred_logger);
Base::calculateReservoirRates(well_state.well(this->index_of_well_));
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state,
DeferredLogger& deferred_logger)
{
updatePrimaryVariables(well_state, deferred_logger);
initPrimaryVariablesEvaluation();
computePerfCellPressDiffs(ebosSimulator);
computeInitialSegmentFluids(ebosSimulator);
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
updateProductivityIndex(const Simulator& ebosSimulator,
const WellProdIndexCalculator& wellPICalc,
WellState& well_state,
DeferredLogger& deferred_logger) const
{
auto fluidState = [&ebosSimulator, this](const int perf)
{
const auto cell_idx = this->well_cells_[perf];
return ebosSimulator.model()
.cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
};
const int np = this->number_of_phases_;
auto setToZero = [np](double* x) -> void
{
std::fill_n(x, np, 0.0);
};
auto addVector = [np](const double* src, double* dest) -> void
{
std::transform(src, src + np, dest, dest, std::plus<>{});
};
auto& ws = well_state.well(this->index_of_well_);
auto& perf_data = ws.perf_data;
auto* connPI = perf_data.prod_index.data();
auto* wellPI = ws.productivity_index.data();
setToZero(wellPI);
const auto preferred_phase = this->well_ecl_.getPreferredPhase();
auto subsetPerfID = 0;
for ( const auto& perf : *this->perf_data_){
auto allPerfID = perf.ecl_index;
auto connPICalc = [&wellPICalc, allPerfID](const double mobility) -> double
{
return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
};
std::vector<Scalar> mob(this->num_components_, 0.0);
getMobilityScalar(ebosSimulator, static_cast<int>(subsetPerfID), mob);
const auto& fs = fluidState(subsetPerfID);
setToZero(connPI);
if (this->isInjector()) {
this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
mob, connPI, deferred_logger);
}
else { // Production or zero flow rate
this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
}
addVector(connPI, wellPI);
++subsetPerfID;
connPI += np;
}
assert (static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
"Internal logic error in processing connections for PI/II");
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
addWellContributions(SparseMatrixAdapter& jacobian) const
{
const auto invDuneD = mswellhelpers::invertWithUMFPack<DiagMatWell, BVectorWell>(this->duneD_, this->duneDSolver_);
// We need to change matrix A as follows
// A -= C^T D^-1 B
// D is a (nseg x nseg) block matrix with (4 x 4) blocks.
// B and C are (nseg x ncells) block matrices with (4 x 4 blocks).
// They have nonzeros at (i, j) only if this well has a
// perforation at cell j connected to segment i. The code
// assumes that no cell is connected to more than one segment,
// i.e. the columns of B/C have no more than one nonzero.
for (size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) {
for (auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) {
const auto row_index = colC.index();
for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
OffDiagMatrixBlockWellType tmp1;
Detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type());
typename SparseMatrixAdapter::MatrixBlock tmp2;
Detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type());
jacobian.addToBlock(row_index, col_index, tmp2);
}
}
}
}
}
template<typename TypeTag>
template<class Value>
void
MultisegmentWell<TypeTag>::
computePerfRate(const Value& pressure_cell,
const Value& rs,
const Value& rv,
const std::vector<Value>& b_perfcells,
const std::vector<Value>& mob_perfcells,
const double Tw,
const int perf,
const Value& segment_pressure,
const Value& segment_density,
const bool& allow_cf,
const std::vector<Value>& cmix_s,
std::vector<Value>& cq_s,
Value& perf_press,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate,
DeferredLogger& deferred_logger) const
{
// pressure difference between the segment and the perforation
const Value perf_seg_press_diff = this->gravity() * segment_density * this->perforation_segment_depth_diffs_[perf];
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
perf_press = pressure_cell - cell_perf_press_diff;
// Pressure drawdown (also used to determine direction of flow)
// TODO: not 100% sure about the sign of the seg_perf_press_diff
const Value drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
// producing perforations
if ( drawdown > 0.0) {
// Do nothing is crossflow is not allowed
if (!allow_cf && this->isInjector()) {
return;
}
// compute component volumetric rates at standard conditions
for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown);
cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const Value cq_s_oil = cq_s[oilCompIdx];
const Value cq_s_gas = cq_s[gasCompIdx];
cq_s[gasCompIdx] += rs * cq_s_oil;
cq_s[oilCompIdx] += rv * cq_s_gas;
}
} else { // injecting perforations
// Do nothing if crossflow is not allowed
if (!allow_cf && this->isProducer()) {
return;
}
// for injecting perforations, we use total mobility
Value total_mob = mob_perfcells[0];
for (int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
total_mob += mob_perfcells[comp_idx];
}
// injection perforations total volume rates
const Value cqt_i = - Tw * (total_mob * drawdown);
// compute volume ratio between connection and at standard conditions
Value volume_ratio = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// Incorporate RS/RV factors if both oil and gas active
// TODO: not sure we use rs rv from the perforation cells when handling injecting perforations
// basically, for injecting perforations, the wellbore is the upstreaming side.
const Value d = 1.0 - rv * rs;
if (getValue(d) == 0.0) {
OPM_DEFLOG_THROW(NumericalIssue, "Zero d value obtained for well " << this->name()
<< " during flux calculation"
<< " with rs " << rs << " and rv " << rv, deferred_logger);
}
const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
} else { // not having gas and oil at the same time
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
}
}
// injecting connections total volumerates at standard conditions
Value cqt_is = cqt_i / volume_ratio;
for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
}
} // end for injection perforations
// calculating the perforation solution gas rate and solution oil rates
if (this->isProducer()) {
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
// s means standard condition, r means reservoir condition
// q_os = q_or * b_o + rv * q_gr * b_g
// q_gs = q_gr * g_g + rs * q_or * b_o
// d = 1.0 - rs * rv
// q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
const double d = 1.0 - getValue(rv) * getValue(rs);
// vaporized oil into gas
// rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
// dissolved of gas in oil
// rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
}
}
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
computePerfRateEval(const IntensiveQuantities& int_quants,
const std::vector<EvalWell>& mob_perfcells,
const double Tw,
const int seg,
const int perf,
const EvalWell& segment_pressure,
const bool& allow_cf,
std::vector<EvalWell>& cq_s,
EvalWell& perf_press,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate,
DeferredLogger& deferred_logger) const
{
const auto& fs = int_quants.fluidState();
const EvalWell pressure_cell = this->extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
const EvalWell rs = this->extendEval(fs.Rs());
const EvalWell rv = this->extendEval(fs.Rv());
// not using number_of_phases_ because of solvent
std::vector<EvalWell> b_perfcells(this->num_components_, 0.0);
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
b_perfcells[compIdx] = this->extendEval(fs.invB(phaseIdx));
}
std::vector<EvalWell> cmix_s(this->numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
cmix_s[comp_idx] = this->surfaceVolumeFraction(seg, comp_idx);
}
this->computePerfRate(pressure_cell,
rs,
rv,
b_perfcells,
mob_perfcells,
Tw,
perf,
segment_pressure,
this->segment_densities_[seg],
allow_cf,
cmix_s,
cq_s,
perf_press,
perf_dis_gas_rate,
perf_vap_oil_rate,
deferred_logger);
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
computePerfRateScalar(const IntensiveQuantities& int_quants,
const std::vector<Scalar>& mob_perfcells,
const double Tw,
const int seg,
const int perf,
const Scalar& segment_pressure,
const bool& allow_cf,
std::vector<Scalar>& cq_s,
DeferredLogger& deferred_logger) const
{
const auto& fs = int_quants.fluidState();
const Scalar pressure_cell = getValue(fs.pressure(FluidSystem::oilPhaseIdx));
const Scalar rs = getValue(fs.Rs());
const Scalar rv = getValue(fs.Rv());
// not using number_of_phases_ because of solvent
std::vector<Scalar> b_perfcells(this->num_components_, 0.0);
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
b_perfcells[compIdx] = getValue(fs.invB(phaseIdx));
}
std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
cmix_s[comp_idx] = getValue(this->surfaceVolumeFraction(seg, comp_idx));
}
Scalar perf_dis_gas_rate = 0.0;
Scalar perf_vap_oil_rate = 0.0;
Scalar perf_press = 0.0;
this->computePerfRate(pressure_cell,
rs,
rv,
b_perfcells,
mob_perfcells,
Tw,
perf,
segment_pressure,
getValue(this->segment_densities_[seg]),
allow_cf,
cmix_s,
cq_s,
perf_press,
perf_dis_gas_rate,
perf_vap_oil_rate,
deferred_logger);
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeSegmentFluidProperties(const Simulator& ebosSimulator)
{
// TODO: the concept of phases and components are rather confusing in this function.
// needs to be addressed sooner or later.
// get the temperature for later use. It is only useful when we are not handling
// thermal related simulation
// basically, it is a single value for all the segments
EvalWell temperature;
EvalWell saltConcentration;
// not sure how to handle the pvt region related to segment
// for the current approach, we use the pvt region of the first perforated cell
// although there are some text indicating using the pvt region of the lowest
// perforated cell
// TODO: later to investigate how to handle the pvt region
int pvt_region_index;
{
// using the first perforated cell
const int cell_idx = this->well_cells_[0];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& fs = intQuants.fluidState();
temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
saltConcentration = this->extendEval(fs.saltConcentration());
pvt_region_index = fs.pvtRegionIndex();
}
this->MSWEval::computeSegmentFluidProperties(temperature,
saltConcentration,
pvt_region_index);
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
getMobilityEval(const Simulator& ebosSimulator,
const int perf,
std::vector<EvalWell>& mob) const
{
// TODO: most of this function, if not the whole function, can be moved to the base class
const int cell_idx = this->well_cells_[perf];
assert (int(mob.size()) == this->num_components_);
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
// either use mobility of the perforation cell or calcualte its own
// based on passing the saturation table index
const int satid = this->saturation_table_number_[perf] - 1;
const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
}
// if (has_solvent) {
// mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
// }
} else {
const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
// reset the satnumvalue back to original
materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
// compute the mobility
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
}
}
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
getMobilityScalar(const Simulator& ebosSimulator,
const int perf,
std::vector<Scalar>& mob) const
{
// TODO: most of this function, if not the whole function, can be moved to the base class
const int cell_idx = this->well_cells_[perf];
assert (int(mob.size()) == this->num_components_);
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
// either use mobility of the perforation cell or calcualte its own
// based on passing the saturation table index
const int satid = this->saturation_table_number_[perf] - 1;
const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
}
// if (has_solvent) {
// mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
// }
} else {
const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
Scalar relativePerms[3] = { 0.0, 0.0, 0.0 };
MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
// reset the satnumvalue back to original
materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
// compute the mobility
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
mob[activeCompIdx] = relativePerms[phaseIdx] / getValue(intQuants.fluidState().viscosity(phaseIdx));
}
}
}
template<typename TypeTag>
double
MultisegmentWell<TypeTag>::
getRefDensity() const
{
return this->segment_densities_[0].value();
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
checkOperabilityUnderBHPLimitProducer(const WellState& /*well_state*/, const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
{
const auto& summaryState = ebos_simulator.vanguard().summaryState();
const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
// Crude but works: default is one atmosphere.
// TODO: a better way to detect whether the BHP is defaulted or not
const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
// if the BHP limit is not defaulted or the well does not have a THP limit
// we need to check the BHP limit
double temp = 0;
for (int p = 0; p < this->number_of_phases_; ++p) {
temp += this->ipr_a_[p] - this->ipr_b_[p] * bhp_limit;
}
if (temp < 0.) {
this->operability_status_.operable_under_only_bhp_limit = false;
}
// checking whether running under BHP limit will violate THP limit
if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
// option 1: calculate well rates based on the BHP limit.
// option 2: stick with the above IPR curve
// we use IPR here
std::vector<double> well_rates_bhp_limit;
computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
const double thp = this->calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, getRefDensity(), deferred_logger);
const double thp_limit = this->getTHPConstraint(summaryState);
if (thp < thp_limit) {
this->operability_status_.obey_thp_limit_under_bhp_limit = false;
}
}
} else {
// defaulted BHP and there is a THP constraint
// default BHP limit is about 1 atm.
// when applied the hydrostatic pressure correction dp,
// most likely we get a negative value (bhp + dp)to search in the VFP table,
// which is not desirable.
// we assume we can operate under defaulted BHP limit and will violate the THP limit
// when operating under defaulted BHP limit.
this->operability_status_.operable_under_only_bhp_limit = true;
this->operability_status_.obey_thp_limit_under_bhp_limit = false;
}
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
updateIPR(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const
{
// TODO: not handling solvent related here for now
// TODO: it only handles the producers for now
// the formular for the injectors are not formulated yet
if (this->isInjector()) {
return;
}
// initialize all the values to be zero to begin with
std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
const int nseg = this->numberOfSegments();
double seg_bhp_press_diff = 0;
double ref_depth = this->ref_depth_;
for (int seg = 0; seg < nseg; ++seg) {
// calculating the perforation rate for each perforation that belongs to this segment
const double segment_depth = this->segmentSet()[seg].depth();
const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth, segment_depth, this->segment_densities_[seg].value(), this->gravity_);
ref_depth = segment_depth;
seg_bhp_press_diff += dp;
for (const int perf : this->segment_perforations_[seg]) {
std::vector<Scalar> mob(this->num_components_, 0.0);
// TODO: mabye we should store the mobility somewhere, so that we only need to calculate it one per iteration
getMobilityScalar(ebos_simulator, perf, mob);
const int cell_idx = this->well_cells_[perf];
const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = int_quantities.fluidState();
// the pressure of the reservoir grid block the well connection is in
// pressure difference between the segment and the perforation
const double perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg].value() * this->perforation_segment_depth_diffs_[perf];
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const double pressure_cell = fs.pressure(FluidSystem::oilPhaseIdx).value();
// calculating the b for the connection
std::vector<double> b_perf(this->num_components_);
for (size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
if (!FluidSystem::phaseIsActive(phase)) {
continue;
}
const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
b_perf[comp_idx] = fs.invB(phase).value();
}
// the pressure difference between the connection and BHP
const double h_perf = cell_perf_press_diff + perf_seg_press_diff + seg_bhp_press_diff;
const double pressure_diff = pressure_cell - h_perf;
// Let us add a check, since the pressure is calculated based on zero value BHP
// it should not be negative anyway. If it is negative, we might need to re-formulate
// to taking into consideration the crossflow here.
if (pressure_diff <= 0.) {
deferred_logger.warning("NON_POSITIVE_DRAWDOWN_IPR",
"non-positive drawdown found when updateIPR for well " + this->name());
}
// the well index associated with the connection
const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
// TODO: there might be some indices related problems here
// phases vs components
// ipr values for the perforation
std::vector<double> ipr_a_perf(this->ipr_a_.size());
std::vector<double> ipr_b_perf(this->ipr_b_.size());
for (int p = 0; p < this->number_of_phases_; ++p) {
const double tw_mob = tw_perf * mob[p] * b_perf[p];
ipr_a_perf[p] += tw_mob * pressure_diff;
ipr_b_perf[p] += tw_mob;
}
// we need to handle the rs and rv when both oil and gas are present
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const double rs = (fs.Rs()).value();
const double rv = (fs.Rv()).value();
const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
ipr_a_perf[gas_comp_idx] += dis_gas_a;
ipr_a_perf[oil_comp_idx] += vap_oil_a;
const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
ipr_b_perf[gas_comp_idx] += dis_gas_b;
ipr_b_perf[oil_comp_idx] += vap_oil_b;
}
for (int p = 0; p < this->number_of_phases_; ++p) {
// TODO: double check the indices here
this->ipr_a_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_a_perf[p];
this->ipr_b_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_b_perf[p];
}
}
}
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& /*well_state*/, DeferredLogger& deferred_logger)
{
const auto& summaryState = ebos_simulator.vanguard().summaryState();
const auto obtain_bhp = computeBhpAtThpLimitProd(ebos_simulator, summaryState, deferred_logger);
if (obtain_bhp) {
this->operability_status_.can_obtain_bhp_with_thp_limit = true;
const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
const double thp_limit = this->getTHPConstraint(summaryState);
if (*obtain_bhp < thp_limit) {
const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
+ " bars is SMALLER than thp limit "
+ std::to_string(unit::convert::to(thp_limit, unit::barsa))
+ " bars as a producer for well " + this->name();
deferred_logger.debug(msg);
}
} else {
// Shutting wells that can not find bhp value from thp
// when under THP control
this->operability_status_.can_obtain_bhp_with_thp_limit = false;
this->operability_status_.obey_bhp_limit_with_thp_limit = false;
if (!this->wellIsStopped()) {
const double thp_limit = this->getTHPConstraint(summaryState);
deferred_logger.debug(" could not find bhp value at thp limit "
+ std::to_string(unit::convert::to(thp_limit, unit::barsa))
+ " bar for well " + this->name() + ", the well might need to be closed ");
}
}
}
template<typename TypeTag>
bool
MultisegmentWell<TypeTag>::
iterateWellEqWithControl(const Simulator& ebosSimulator,
const double dt,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger)
{
if (!this->isOperable() && !this->wellIsStopped()) return true;
const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
const WellState well_state0 = well_state;
const std::vector<Scalar> residuals0 = this->getWellResiduals(Base::B_avg_, deferred_logger);
std::vector<std::vector<Scalar> > residual_history;
std::vector<double> measure_history;
int it = 0;
// relaxation factor
double relaxation_factor = 1.;
const double min_relaxation_factor = 0.6;
bool converged = false;
int stagnate_count = 0;
bool relax_convergence = false;
for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
if (it > this->param_.strict_inner_iter_wells_)
relax_convergence = true;
const auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
if (report.converged()) {
converged = true;
break;
}
residual_history.push_back(this->getWellResiduals(Base::B_avg_, deferred_logger));
measure_history.push_back(this->getResidualMeasureValue(well_state,
residual_history[it],
this->param_.tolerance_wells_,
this->param_.tolerance_pressure_ms_wells_,
deferred_logger) );
bool is_oscillate = false;
bool is_stagnate = false;
this->detectOscillations(measure_history, it, is_oscillate, is_stagnate);
// TODO: maybe we should have more sophiscated strategy to recover the relaxation factor,
// for example, to recover it to be bigger
if (is_oscillate || is_stagnate) {
// HACK!
std::ostringstream sstr;
if (relaxation_factor == min_relaxation_factor) {
// Still stagnating, terminate iterations if 5 iterations pass.
++stagnate_count;
if (stagnate_count == 6) {
sstr << " well " << this->name() << " observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
const auto reportStag = getWellConvergence(well_state, Base::B_avg_, deferred_logger, true);
if (reportStag.converged()) {
converged = true;
sstr << " well " << this->name() << " manages to get converged with relaxed tolerances in " << it << " inner iterations";
deferred_logger.debug(sstr.str());
return converged;
}
}
}
// a factor value to reduce the relaxation_factor
const double reduction_mutliplier = 0.9;
relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
// debug output
if (is_stagnate) {
sstr << " well " << this->name() << " observes stagnation in inner iteration " << it << "\n";
}
if (is_oscillate) {
sstr << " well " << this->name() << " observes oscillation in inner iteration " << it << "\n";
}
sstr << " relaxation_factor is " << relaxation_factor << " now\n";
deferred_logger.debug(sstr.str());
}
updateWellState(dx_well, well_state, deferred_logger, relaxation_factor);
initPrimaryVariablesEvaluation();
}
// TODO: we should decide whether to keep the updated well_state, or recover to use the old well_state
if (converged) {
std::ostringstream sstr;
sstr << " Well " << this->name() << " converged in " << it << " inner iterations.";
if (relax_convergence)
sstr << " (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ << " iterations)";
deferred_logger.debug(sstr.str());
} else {
std::ostringstream sstr;
sstr << " Well " << this->name() << " did not converge in " << it << " inner iterations.";
#define EXTRA_DEBUG_MSW 0
#if EXTRA_DEBUG_MSW
sstr << "***** Outputting the residual history for well " << this->name() << " during inner iterations:";
for (int i = 0; i < it; ++i) {
const auto& residual = residual_history[i];
sstr << " residual at " << i << "th iteration ";
for (const auto& res : residual) {
sstr << " " << res;
}
sstr << " " << measure_history[i] << " \n";
}
#endif
deferred_logger.debug(sstr.str());
}
return converged;
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
const double dt,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger)
{
if (!this->isOperable() && !this->wellIsStopped()) return;
// update the upwinding segments
this->updateUpwindingSegments();
// calculate the fluid properties needed.
computeSegmentFluidProperties(ebosSimulator);
// clear all entries
this->duneB_ = 0.0;
this->duneC_ = 0.0;
this->duneD_ = 0.0;
this->resWell_ = 0.0;
this->duneDSolver_.reset();
auto& ws = well_state.well(this->index_of_well_);
ws.dissolved_gas_rate = 0;
ws.vaporized_oil_rate = 0;
// for the black oil cases, there will be four equations,
// the first three of them are the mass balance equations, the last one is the pressure equations.
//
// but for the top segment, the pressure equation will be the well control equation, and the other three will be the same.
const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
const int nseg = this->numberOfSegments();
for (int seg = 0; seg < nseg; ++seg) {
// calculating the accumulation term
// TODO: without considering the efficiencty factor for now
{
const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);
// Add a regularization_factor to increase the accumulation term
// This will make the system less stiff and help convergence for
// difficult cases
const Scalar regularization_factor = this->param_.regularization_factor_ms_wells_;
// for each component
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx)
- segment_fluid_initial_[seg][comp_idx]) / dt;
this->resWell_[seg][comp_idx] += accumulation_term.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
this->duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq);
}
}
}
// considering the contributions due to flowing out from the segment
{
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) * this->well_efficiency_factor_;
const int seg_upwind = this->upwinding_segments_[seg];
// segment_rate contains the derivatives with respect to GTotal in seg,
// and WFrac and GFrac in seg_upwind
this->resWell_[seg][comp_idx] -= segment_rate.value();
this->duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
this->duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq);
}
// pressure derivative should be zero
}
}
// considering the contributions from the inlet segments
{
for (const int inlet : this->segment_inlets_[seg]) {
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const EvalWell inlet_rate = this->getSegmentRateUpwinding(inlet, comp_idx) * this->well_efficiency_factor_;
const int inlet_upwind = this->upwinding_segments_[inlet];
// inlet_rate contains the derivatives with respect to GTotal in inlet,
// and WFrac and GFrac in inlet_upwind
this->resWell_[seg][comp_idx] += inlet_rate.value();
this->duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
this->duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq);
}
// pressure derivative should be zero
}
}
}
// calculating the perforation rate for each perforation that belongs to this segment
const EvalWell seg_pressure = this->getSegmentPressure(seg);
auto& perf_data = ws.perf_data;
auto& perf_rates = perf_data.phase_rates;
auto& perf_press_state = perf_data.pressure;
for (const int perf : this->segment_perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<EvalWell> mob(this->num_components_, 0.0);
getMobilityEval(ebosSimulator, perf, mob);
const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
const double Tw = this->well_index_[perf] * trans_mult;
std::vector<EvalWell> cq_s(this->num_components_, 0.0);
EvalWell perf_press;
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRateEval(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
// updating the solution gas rate and solution oil rate
if (this->isProducer()) {
ws.dissolved_gas_rate += perf_dis_gas_rate;
ws.vaporized_oil_rate += perf_vap_oil_rate;
}
// store the perf pressure and rates
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
perf_rates[perf*this->number_of_phases_ + this->ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
}
perf_press_state[perf] = perf_press.value();
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
// the cq_s entering mass balance equations need to consider the efficiency factors.
const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
// subtract sum of phase fluxes in the well equations.
this->resWell_[seg][comp_idx] += cq_s_effective.value();
// assemble the jacobians
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
// also need to consider the efficiency factor when manipulating the jacobians.
this->duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq); // intput in transformed matrix
// the index name for the D should be eq_idx / pv_idx
this->duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq);
}
for (int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) {
// also need to consider the efficiency factor when manipulating the jacobians.
this->duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
}
}
}
// the fourth dequation, the pressure drop equation
if (seg == 0) { // top segment, pressure equation is the control equation
const auto& summaryState = ebosSimulator.vanguard().summaryState();
const Schedule& schedule = ebosSimulator.vanguard().schedule();
this->assembleControlEq(well_state,
group_state,
schedule,
summaryState,
inj_controls,
prod_controls,
getRefDensity(),
deferred_logger);
} else {
const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
this->assemblePressureEq(seg, unit_system, well_state, deferred_logger);
}
}
}
template<typename TypeTag>
bool
MultisegmentWell<TypeTag>::
openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const
{
return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
}
template<typename TypeTag>
bool
MultisegmentWell<TypeTag>::
allDrawDownWrongDirection(const Simulator& ebos_simulator) const
{
bool all_drawdown_wrong_direction = true;
const int nseg = this->numberOfSegments();
for (int seg = 0; seg < nseg; ++seg) {
const EvalWell segment_pressure = this->getSegmentPressure(seg);
for (const int perf : this->segment_perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = intQuants.fluidState();
// pressure difference between the segment and the perforation
const EvalWell perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf];
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const double pressure_cell = (fs.pressure(FluidSystem::oilPhaseIdx)).value();
const double perf_press = pressure_cell - cell_perf_press_diff;
// Pressure drawdown (also used to determine direction of flow)
// TODO: not 100% sure about the sign of the seg_perf_press_diff
const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
// for now, if there is one perforation can produce/inject in the correct
// direction, we consider this well can still produce/inject.
// TODO: it can be more complicated than this to cause wrong-signed rates
if ( (drawdown < 0. && this->isInjector()) ||
(drawdown > 0. && this->isProducer()) ) {
all_drawdown_wrong_direction = false;
break;
}
}
}
return all_drawdown_wrong_direction;
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
updateWaterThroughput(const double /*dt*/, WellState& /*well_state*/) const
{
}
template<typename TypeTag>
typename MultisegmentWell<TypeTag>::EvalWell
MultisegmentWell<TypeTag>::
getSegmentSurfaceVolume(const Simulator& ebos_simulator, const int seg_idx) const
{
EvalWell temperature;
EvalWell saltConcentration;
int pvt_region_index;
{
// using the pvt region of first perforated cell
// TODO: it should be a member of the WellInterface, initialized properly
const int cell_idx = this->well_cells_[0];
const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& fs = intQuants.fluidState();
temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
saltConcentration = this->extendEval(fs.saltConcentration());
pvt_region_index = fs.pvtRegionIndex();
}
return this->MSWEval::getSegmentSurfaceVolume(temperature,
saltConcentration,
pvt_region_index,
seg_idx);
}
template<typename TypeTag>
std::optional<double>
MultisegmentWell<TypeTag>::
computeBhpAtThpLimitProd(const Simulator& ebos_simulator,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const
{
// Make the frates() function.
auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
// Not solving the well equations here, which means we are
// calculating at the current Fg/Fw values of the
// well. This does not matter unless the well is
// crossflowing, and then it is likely still a good
// approximation.
std::vector<double> rates(3);
computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
return rates;
};
auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitProd(frates,
summary_state,
maxPerfPress(ebos_simulator),
getRefDensity(),
deferred_logger);
if(bhpAtLimit)
return bhpAtLimit;
auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
// Solver the well iterations to see if we are
// able to get a solution with an update
// solution
std::vector<double> rates(3);
computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
return rates;
};
return this->MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitProd(fratesIter,
summary_state,
maxPerfPress(ebos_simulator),
getRefDensity(),
deferred_logger);
}
template<typename TypeTag>
std::optional<double>
MultisegmentWell<TypeTag>::
computeBhpAtThpLimitInj(const Simulator& ebos_simulator,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const
{
// Make the frates() function.
auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
// Not solving the well equations here, which means we are
// calculating at the current Fg/Fw values of the
// well. This does not matter unless the well is
// crossflowing, and then it is likely still a good
// approximation.
std::vector<double> rates(3);
computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
return rates;
};
auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitInj(frates,
summary_state,
getRefDensity(),
deferred_logger);
if(bhpAtLimit)
return bhpAtLimit;
auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
// Solver the well iterations to see if we are
// able to get a solution with an update
// solution
std::vector<double> rates(3);
computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
return rates;
};
return this->MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitInj(fratesIter, summary_state, getRefDensity(), deferred_logger);
}
template<typename TypeTag>
double
MultisegmentWell<TypeTag>::
maxPerfPress(const Simulator& ebos_simulator) const
{
double max_pressure = 0.0;
const int nseg = this->numberOfSegments();
for (int seg = 0; seg < nseg; ++seg) {
for (const int perf : this->segment_perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = int_quants.fluidState();
double pressure_cell = fs.pressure(FluidSystem::oilPhaseIdx).value();
max_pressure = std::max(max_pressure, pressure_cell);
}
}
return max_pressure;
}
template<typename TypeTag>
std::vector<double>
MultisegmentWell<TypeTag>::
computeCurrentWellRates(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger) const
{
// Calculate the rates that follow from the current primary variables.
std::vector<Scalar> well_q_s(this->num_components_, 0.0);
const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
const int nseg = this->numberOfSegments();
for (int seg = 0; seg < nseg; ++seg) {
// calculating the perforation rate for each perforation that belongs to this segment
const Scalar seg_pressure = getValue(this->getSegmentPressure(seg));
for (const int perf : this->segment_perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<Scalar> mob(this->num_components_, 0.0);
getMobilityScalar(ebosSimulator, perf, mob);
const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
const double Tw = this->well_index_[perf] * trans_mult;
std::vector<Scalar> cq_s(this->num_components_, 0.0);
computePerfRateScalar(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, deferred_logger);
for (int comp = 0; comp < this->num_components_; ++comp) {
well_q_s[comp] += cq_s[comp];
}
}
}
return well_q_s;
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeConnLevelProdInd(const typename MultisegmentWell<TypeTag>::FluidState& fs,
const std::function<double(const double)>& connPICalc,
const std::vector<Scalar>& mobility,
double* connPI) const
{
const auto& pu = this->phaseUsage();
const int np = this->number_of_phases_;
for (int p = 0; p < np; ++p) {
// Note: E100's notion of PI value phase mobility includes
// the reciprocal FVF.
const auto connMob =
mobility[ this->flowPhaseToEbosCompIdx(p) ]
* fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
connPI[p] = connPICalc(connMob);
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
{
const auto io = pu.phase_pos[Oil];
const auto ig = pu.phase_pos[Gas];
const auto vapoil = connPI[ig] * fs.Rv().value();
const auto disgas = connPI[io] * fs.Rs().value();
connPI[io] += vapoil;
connPI[ig] += disgas;
}
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeConnLevelInjInd(const typename MultisegmentWell<TypeTag>::FluidState& fs,
const Phase preferred_phase,
const std::function<double(const double)>& connIICalc,
const std::vector<Scalar>& mobility,
double* connII,
DeferredLogger& deferred_logger) const
{
// Assumes single phase injection
const auto& pu = this->phaseUsage();
auto phase_pos = 0;
if (preferred_phase == Phase::GAS) {
phase_pos = pu.phase_pos[Gas];
}
else if (preferred_phase == Phase::OIL) {
phase_pos = pu.phase_pos[Oil];
}
else if (preferred_phase == Phase::WATER) {
phase_pos = pu.phase_pos[Water];
}
else {
OPM_DEFLOG_THROW(NotImplemented,
"Unsupported Injector Type ("
<< static_cast<int>(preferred_phase)
<< ") for well " << this->name()
<< " during connection I.I. calculation",
deferred_logger);
}
const Scalar mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
}
} // namespace Opm
Reference in New Issue View Git Blame Copy Permalink