Refactor to put updateWellStateRates() in base class.

This commit is contained in:
Atgeirr Flø Rasmussen 2020-10-15 14:15:05 +02:00
parent 52c695937b
commit 7e87ea3200
8 changed files with 124 additions and 75 deletions

View File

@ -1290,16 +1290,6 @@ namespace Opm {
}
}
}
// Output debug info.
if (terminal_output_) {
std::ostringstream oss;
oss << "Node pressures in network nodes, in bar:\n";
for (const auto& [name, pressure] : node_pressures_) {
oss << name << " " << unit::convert::to(pressure, unit::barsa) << '\n';
}
OpmLog::debug(oss.str());
}
}

View File

@ -176,6 +176,9 @@ namespace Opm
int numberOfPerforations() const;
virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger) const override;
protected:
int number_segments_;

View File

@ -3819,4 +3819,40 @@ namespace Opm
const double sign = mass_rate <= 0. ? 1.0 : -1.0;
return sign * (friction_pressure_loss + constriction_pressure_loss);
}
template<typename TypeTag>
std::vector<double>
MultisegmentWell<TypeTag>::
computeCurrentWellRates(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger) const
{
#if 0
// Calculate the rates that follow from the current primary variables.
std::vector<EvalWell> well_q_s(num_components_, 0.0);
const EvalWell& bhp = getBhp();
const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const int cell_idx = well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<EvalWell> mob(num_components_, 0.0);
getMobility(ebosSimulator, perf, mob);
std::vector<EvalWell> cq_s(num_components_, 0.0);
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
const double Tw = well_index_[perf] * trans_mult;
// computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
// cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
for (int comp = 0; comp < num_components_; ++comp) {
well_q_s[comp] += cq_s[comp];
}
}
#endif
}
} // namespace Opm

View File

@ -296,9 +296,8 @@ namespace Opm
using Base::phaseUsage;
using Base::vfp_properties_;
virtual void updateWellStateRates(const Simulator& ebosSimulator,
WellState& well_state,
DeferredLogger& deferred_logger) const override;
virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger) const override;
protected:

View File

@ -4031,30 +4031,11 @@ namespace Opm
template<typename TypeTag>
void
std::vector<double>
StandardWell<TypeTag>::
updateWellStateRates(const Simulator& ebosSimulator,
WellState& well_state,
DeferredLogger& deferred_logger) const
computeCurrentWellRates(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger) const
{
// Check if the rates of this well only are single-phase, do nothing
// if more than one nonzero rate.
int nonzero_rate_index = -1;
for (int p = 0; p < number_of_phases_; ++p) {
if (well_state.wellRates()[index_of_well_ * number_of_phases_ + p] != 0.0) {
if (nonzero_rate_index == -1) {
nonzero_rate_index = p;
} else {
// More than one nonzero rate.
return;
}
}
}
if (nonzero_rate_index == -1) {
// No nonzero rates.
return;
}
// Calculate the rates that follow from the current primary variables.
std::vector<EvalWell> well_q_s(num_components_, {numWellEq_ + numEq, 0.});
const EvalWell& bhp = getBhp();
@ -4075,39 +4056,11 @@ namespace Opm
well_q_s[comp] += cq_s[comp];
}
}
// We must keep in mind that component and phase indices are different here.
// Therefore we set up a mapping to make the last code block simpler.
const auto& pu = phaseUsage();
const int wpi = Water;
const int opi = Oil;
const int gpi = Gas;
const unsigned int wci = FluidSystem::waterCompIdx;
const unsigned int oci = FluidSystem::oilCompIdx;
const unsigned int gci = FluidSystem::gasCompIdx;
std::pair<int, unsigned int> phase_comp[3] = { {wpi, wci},
{opi, oci},
{gpi, gci} };
std::vector<int> phase_to_comp(number_of_phases_, -1);
for (const auto& pc : phase_comp) {
if (pu.phase_used[pc.first]) {
const int phase_idx = pu.phase_pos[pc.first];
const int comp_idx = Indices::canonicalToActiveComponentIndex(pc.second);
phase_to_comp[phase_idx] = comp_idx;
}
}
// Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
const double initial_nonzero_rate = well_state.wellRates()[index_of_well_ * number_of_phases_ + nonzero_rate_index];
const int comp_idx_nz = phase_to_comp[nonzero_rate_index];
for (int p = 0; p < number_of_phases_; ++p) {
if (p != nonzero_rate_index) {
const int comp_idx = phase_to_comp[p];
double& rate = well_state.wellRates()[index_of_well_ * number_of_phases_ + p];
rate = (initial_nonzero_rate/well_q_s[comp_idx_nz].value()) * (well_q_s[comp_idx].value());
}
std::vector<double> well_q_s_noderiv(well_q_s.size());
for (int comp = 0; comp < num_components_; ++comp) {
well_q_s_noderiv[comp] = well_q_s[comp].value();
}
return well_q_s_noderiv;
}

View File

@ -679,7 +679,7 @@ namespace WellGroupHelpers
rates[BlackoilPhases::Vapour],
up_press,
alq);
#define EXTRA_DEBUG_NETWORK 1
#define EXTRA_DEBUG_NETWORK 0
#if EXTRA_DEBUG_NETWORK
std::ostringstream oss;
oss << "parent: " << (*upbranch).uptree_node() << " child: " << node

View File

@ -275,13 +275,17 @@ namespace Opm
// update perforation water throughput based on solved water rate
virtual void updateWaterThroughput(const double dt, WellState& well_state) const = 0;
virtual void updateWellStateRates(const Simulator& ebosSimulator,
WellState& well_state,
DeferredLogger& deferred_logger) const
{
// Default: do nothing.
// TODO: make this a pure virtual.
}
/// Compute well rates based on current reservoir conditions and well variables.
/// Used in updateWellStateRates().
virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger) const = 0;
/// Modify the well_state's rates if there is only one nonzero rate.
/// If so, that rate is kept as is, but the others are set proportionally
/// to the rates returned by computeCurrentWellRates().
void updateWellStateRates(const Simulator& ebosSimulator,
WellState& well_state,
DeferredLogger& deferred_logger) const;
void stopWell() {
wellIsStopped_ = true;

View File

@ -2287,4 +2287,68 @@ namespace Opm
}
template <typename TypeTag>
void
WellInterface<TypeTag>::
updateWellStateRates(const Simulator& ebosSimulator,
WellState& well_state,
DeferredLogger& deferred_logger) const
{
// Check if the rates of this well only are single-phase, do nothing
// if more than one nonzero rate.
int nonzero_rate_index = -1;
for (int p = 0; p < number_of_phases_; ++p) {
if (well_state.wellRates()[index_of_well_ * number_of_phases_ + p] != 0.0) {
if (nonzero_rate_index == -1) {
nonzero_rate_index = p;
} else {
// More than one nonzero rate.
return;
}
}
}
if (nonzero_rate_index == -1) {
// No nonzero rates.
return;
}
// Calculate the rates that follow from the current primary variables.
std::vector<double> well_q_s = computeCurrentWellRates(ebosSimulator, deferred_logger);
// We must keep in mind that component and phase indices are different here.
// Therefore we set up a mapping to make the last code block simpler.
const auto& pu = phaseUsage();
const int wpi = Water;
const int opi = Oil;
const int gpi = Gas;
const unsigned int wci = FluidSystem::waterCompIdx;
const unsigned int oci = FluidSystem::oilCompIdx;
const unsigned int gci = FluidSystem::gasCompIdx;
std::pair<int, unsigned int> phase_comp[3] = { {wpi, wci},
{opi, oci},
{gpi, gci} };
std::vector<int> phase_to_comp(number_of_phases_, -1);
for (const auto& pc : phase_comp) {
if (pu.phase_used[pc.first]) {
const int phase_idx = pu.phase_pos[pc.first];
const int comp_idx = Indices::canonicalToActiveComponentIndex(pc.second);
phase_to_comp[phase_idx] = comp_idx;
}
}
// Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
const double initial_nonzero_rate = well_state.wellRates()[index_of_well_ * number_of_phases_ + nonzero_rate_index];
const int comp_idx_nz = phase_to_comp[nonzero_rate_index];
for (int p = 0; p < number_of_phases_; ++p) {
if (p != nonzero_rate_index) {
const int comp_idx = phase_to_comp[p];
double& rate = well_state.wellRates()[index_of_well_ * number_of_phases_ + p];
rate = (initial_nonzero_rate/well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
}
}
}
} // namespace Opm