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

Use robustSolveBhpAtThpLimitProd() for potentials and initialization.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2019-10-03 16:42:56 +02:00
committed by Tor Harald Sandve
parent 7581275b03
commit 2b9f5fb4fb
2 changed files with 37 additions and 182 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
opm/simulators/wells/StandardWell.hpp
View File

@@ -361,10 +361,7 @@ namespace Opm
Opm::DeferredLogger& deferred_logger);
std::vector<double> computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential,
Opm::DeferredLogger& deferred_logger);
Opm::DeferredLogger& deferred_logger) const;
template <class ValueType>
ValueType calculateBhpFromThp(const std::vector<ValueType>& rates, const Well2& well, const SummaryState& summaryState, Opm::DeferredLogger& deferred_logger) const;

214
opm/simulators/wells/StandardWell_impl.hpp
View File

@@ -239,10 +239,6 @@ namespace Opm
StandardWell<TypeTag>::
wellVolumeFraction(const unsigned compIdx) const
{
if (FluidSystem::numActivePhases() == 1) {
return EvalWell(numWellEq_ + numEq, 1.0);
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
return primary_variables_evaluation_[WFrac];
}
@@ -738,7 +734,6 @@ namespace Opm
template <typename TypeTag>
void
StandardWell<TypeTag>::
@@ -1328,6 +1323,7 @@ namespace Opm
processFractions();
}
// updating the total rates Q_t
const double relaxation_factor_rate = relaxationFactorRate(old_primary_variables, dwells);
primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
@@ -1501,7 +1497,6 @@ namespace Opm
F[p] = 0.;
}
}
// F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent.
// More testing is needed to make sure this is correct for well groups and THP.
if (has_solvent){
@@ -1777,24 +1772,10 @@ namespace Opm
}
case Well2::ProducerCMode::THP:
{
std::vector<double> rates(3, 0.0);
auto eq = [this, &ebos_simulator, &rates, &deferred_logger, &well, &summaryState](const double bhp) {
computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
const double bhp_calculated = calculateBhpFromThp(rates, well, summaryState, deferred_logger);
return bhp_calculated - bhp;
};
try {
int iteration = 0;
const double solved_bhp = RegulaFalsi<>::solve(eq,
controls.bhp_limit,
400e5,
100,
1000.0,
iteration);
well_state.bhp()[well_index] = solved_bhp;
}
catch (...) {
auto bhp = robustSolveBhpAtThpLimitProd(ebos_simulator, summaryState, deferred_logger);
if (bhp) {
well_state.bhp()[well_index] = *bhp;
} else {
deferred_logger.warning("FAILURE_GETTING_CONVERGED_BHP",
"Failed to find BHP when switching to THP control for well " + name());
}
@@ -2087,9 +2068,8 @@ namespace Opm
const int cell_idx = well_cells_[perf];
const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& fs = intQuants.fluidState();
Eval perf_pressure = getPerfCellPressure(fs);
const double pressure = perf_pressure.value();
const double pressure = (fs.pressure(FluidSystem::oilPhaseIdx)).value();
const double bhp = getBhp().value();
// Pressure drawdown (also used to determine direction of flow)
@@ -2808,130 +2788,22 @@ namespace Opm
template<typename TypeTag>
std::vector<double>
StandardWell<TypeTag>::
computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential,
Opm::DeferredLogger& deferred_logger)
computeWellPotentialWithTHP(const Simulator& ebos_simulator,
Opm::DeferredLogger& deferred_logger) const
{
// TODO: pay attention to the situation that finally the potential is calculated based on the bhp control
// TODO: should we consider the bhp constraints during the iterative process?
const int np = number_of_phases_;
std::vector<double> potentials(number_of_phases_, 0.0);
const auto& summary_state = ebos_simulator.vanguard().summaryState();
assert( np == int(initial_potential.size()) );
std::vector<double> potentials = initial_potential;
std::vector<double> old_potentials = potentials; // keeping track of the old potentials
double bhp = initial_bhp;
double old_bhp = bhp;
bool converged = false;
const int max_iteration = 1000;
const double bhp_tolerance = 1000.; // 1000 pascal
int iteration = 0;
const auto& well = Base::wellEcl();
const auto& summaryState = ebosSimulator.vanguard().summaryState();
std::vector<double> rates(3, 0.0);
auto eq = [this, &ebosSimulator, &rates, &deferred_logger, &well, &summaryState](const double bhp) {
computeWellRatesWithBhp(ebosSimulator, bhp, rates, deferred_logger);
const double bhp_calculated = calculateBhpFromThp(rates, well, summaryState, deferred_logger);
return bhp_calculated - bhp;
};
try {
const double solved_bhp = RegulaFalsi<>::solve(eq,
initial_bhp,
400e5,
max_iteration,
bhp_tolerance,
iteration);
computeWellRatesWithBhp(ebosSimulator, solved_bhp, potentials, deferred_logger);
}
catch (...) {
auto bhp_at_thp_limit = robustSolveBhpAtThpLimitProd(ebos_simulator, summary_state, deferred_logger);
if (bhp_at_thp_limit) {
const auto& controls = well_ecl_.productionControls(summary_state);
const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
} else {
deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
"Failed in getting converged for the potential calculation for well " + name());
"Failed in getting converged for the potential calculation for well " + name());
}
/*
while ( !converged && iteration < max_iteration ) {
// for each iteration, we calculate the bhp based on the rates/potentials with thp constraints
// with considering the bhp value from the bhp limits. At the beginning of each iteration,
// we initialize the bhp to be the bhp value from the bhp limits. Then based on the bhp values calculated
// from the thp constraints, we decide the effective bhp value for well potential calculation.
bhp = initial_bhp;
// The number of the well controls/constraints
const int nwc = well_controls_get_num(well_controls_);
const Opm::PhaseUsage& pu = phaseUsage();
std::vector<double> rates(3, 0.0);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = potentials[pu.phase_pos[ Water ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = potentials[pu.phase_pos[ Oil ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = potentials[pu.phase_pos[ Gas ] ];
}
const double bhp_calculated = calculateBhpFromThp(rates, well, summaryState, deferred_logger);
if (well_type_ == INJECTOR && bhp_calculated < bhp ) {
bhp = bhp_calculated;
}
if (well_type_ == PRODUCER && bhp_calculated > bhp) {
bhp = bhp_calculated;
}
// there should be always some available bhp/thp constraints there
if (std::isinf(bhp) || std::isnan(bhp)) {
OPM_DEFLOG_THROW(std::runtime_error, "Unvalid bhp value obtained during the potential calculation for well " << name(), deferred_logger);
}
converged = std::abs(old_bhp - bhp) < bhp_tolerance;
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, potentials, deferred_logger);
// checking whether the potentials have valid values
for (const double value : potentials) {
if (std::isinf(value) || std::isnan(value)) {
OPM_DEFLOG_THROW(std::runtime_error, "Unvalid potential value obtained during the potential calculation for well " << name(), deferred_logger);
}
}
if (!converged) {
old_bhp = bhp;
for (int p = 0; p < np; ++p) {
// TODO: improve the interpolation, will it always be valid with the way below?
// TODO: finding better paramters, better iteration strategy for better convergence rate.
const double potential_update_damping_factor = 0.9;
potentials[p] = potential_update_damping_factor * potentials[p] + (1.0 - potential_update_damping_factor) * old_potentials[p];
old_potentials[p] = potentials[p];
}
}
++iteration;
}
if (!converged) {
// TODO: temporary approach to recover the running of prediction case
// more sophisicated treatment should be developed to handle the situation in a general way
// OPM_DEFLOG_THROW(std::runtime_error, "Failed in getting converged for the potential calculation for well " << name(), deferred_logger);
deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
"Failed in getting converged for the potential calculation for well " + name());
}
*/
return potentials;
}
@@ -2963,35 +2835,17 @@ namespace Opm
const int np = number_of_phases_;
well_potentials.resize(np, 0.0);
// get the bhp value based on the bhp constraints
// does the well have a THP related constraint?
const auto& summaryState = ebosSimulator.vanguard().summaryState();
const double bhp = well.mostStrictBhpFromBhpLimits(summaryState);
// does the well have a THP related constraint?
if ( !well.Base::wellHasTHPConstraints(summaryState) ) {
// get the bhp value based on the bhp constraints
const double bhp = well.mostStrictBhpFromBhpLimits(summaryState);
assert(std::abs(bhp) != std::numeric_limits<double>::max());
well.computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
} else {
// the well has a THP related constraint
// checking whether a well is newly added, it only happens at the beginning of the report step
if ( !well_state.effectiveEventsOccurred(index_of_well_) ) {
for (int p = 0; p < np; ++p) {
// This is dangerous for new added well
// since we are not handling the initialization correctly for now
well_potentials[p] = well_state.wellRates()[index_of_well_ * np + p];
}
} else {
// We need to generate a reasonable rates to start the iteration process
well.computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
for (double& value : well_potentials) {
// make the value a little safer in case the BHP limits are default ones
// TODO: a better way should be a better rescaling based on the investigation of the VFP table.
const double rate_safety_scaling_factor = 0.00001;
value *= rate_safety_scaling_factor;
}
}
well_potentials = well.computeWellPotentialWithTHP(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
well_potentials = well.computeWellPotentialWithTHP(ebosSimulator, deferred_logger);
}
}
@@ -3099,8 +2953,6 @@ namespace Opm
template<typename TypeTag>
template<class ValueType>
ValueType
@@ -3122,10 +2974,6 @@ namespace Opm
const ValueType liquid = rates[Oil];
const ValueType vapour = rates[Gas];
const int vfp = well_controls_iget_vfp(well_controls_, control_index);
const double& thp = well_controls_iget_target(well_controls_, control_index);
const double& alq = well_controls_iget_alq(well_controls_, control_index);
// pick the density in the top layer
// TODO: it is possible it should be a Evaluation
const double rho = perf_densities_[0];
@@ -3146,6 +2994,9 @@ namespace Opm
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger);
}
}
@@ -3296,7 +3147,7 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian) const
StandardWell<TypeTag>::addWellContributions(Mat& mat) const
{
// We need to change matrx A as follows
// A -= C^T D^-1 B
@@ -3304,17 +3155,24 @@ namespace Opm
// B and C have 1 row, nc colums and nonzero
// at (0,j) only if this well has a perforation at cell j.
typename SparseMatrixAdapter::MatrixBlock tmpMat;
Dune::DynamicMatrix<Scalar> tmp;
for ( auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC )
{
const auto row_index = colC.index();
auto& row = mat[row_index];
auto col = row.begin();
for ( auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB )
{
const auto col_index = colB.index();
// Move col to index col_index
while ( col != row.end() && col.index() < col_index ) ++col;
assert(col != row.end() && col.index() == col_index);
Dune::DynamicMatrix<Scalar> tmp;
Detail::multMatrix(invDuneD_[0][0], (*colB), tmp);
Detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat);
jacobian.addToBlock( row_index, colB.index(), tmpMat );
typename Mat::block_type tmp1;
Detail::multMatrixTransposed((*colC), tmp, tmp1);
(*col) -= tmp1;
}
}
}