OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-26 01:01:00 -06:00
Code Issues Packages Projects Releases Wiki Activity

adding assembleControlEq to StandardWell

Browse Source 
solvent model needs to be handled for the injector control.
This commit is contained in:
Kai Bao 2018-05-09 12:16:33 +02:00
parent 62947c85b9
commit 3b70f5004f
2 changed files with 132 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
opm/autodiff/StandardWell.hpp
View File

@ -76,7 +76,6 @@ namespace Opm
static const int Bhp = gasoil? 2 : 3;
static const int SFrac = !has_solvent ? -1000 : Bhp + 1;
using typename Base::Scalar;
using typename Base::ConvergenceReport;
@ -245,6 +244,8 @@ namespace Opm
// TODO: it is also possible to be moved to the base class.
EvalWell getQs(const int comp_idx) const;
const EvalWell& getGTotal() const;
EvalWell wellVolumeFractionScaled(const int phase) const;
EvalWell wellVolumeFraction(const unsigned compIdx) const;
@ -328,6 +329,8 @@ namespace Opm
const WellState& well_state) const;
void updateWellStateFromPrimaryVariables(WellState& well_state) const;
void assembleControlEq();
};
}

144
opm/autodiff/StandardWell_impl.hpp
View File

@ -35,7 +35,7 @@ namespace Opm
, perf_pressure_diffs_(number_of_perforations_)
, primary_variables_(numWellEq, 0.0)
, primary_variables_evaluation_(numWellEq) // the number of the primary variables
, F0_(numWellEq)
, F0_(num_components_)
{
duneB_.setBuildMode( OffDiagMatWell::row_wise );
duneC_.setBuildMode( OffDiagMatWell::row_wise );
@ -133,6 +133,18 @@ namespace Opm
template<typename TypeTag>
const typename StandardWell<TypeTag>::EvalWell&
StandardWell<TypeTag>::
getGTotal() const
{
return primary_variables_evaluation_[GTotal];
}
template<typename TypeTag>
typename StandardWell<TypeTag>::EvalWell
StandardWell<TypeTag>::
@ -146,23 +158,24 @@ namespace Opm
assert(comp_idx < num_components_);
if (well_type_ == INJECTOR) { // only single phase injection
// TODO: using comp_frac here is dangerous, it should be changed later
// Most likely, it should be changed to use distr, or at least, we need to update comp_frac_ based on distr
// while solvent might complicate the situation
EvalWell qs = 0.0;
const auto pu = phaseUsage();
const int legacyCompIdx = ebosCompIdxToFlowCompIdx(comp_idx);
if (has_solvent) {
// TODO: investigate whether the use of the comp_frac is justified.
// The usage of the comp_frac is not correct, which should be changed later.
double comp_frac = 0.0;
if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent
comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent();
} else if (legacyCompIdx == pu.phase_pos[ Gas ]) {
comp_frac = comp_frac_[legacyCompIdx] * (1.0 - wsolvent());
} else {
comp_frac = comp_frac_[legacyCompIdx];
double comp_frac = 0.0;
if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent
comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent();
} else if (legacyCompIdx == pu.phase_pos[ Gas ]) {
comp_frac = comp_frac_[legacyCompIdx];
if (has_solvent) {
comp_frac *= (1.0 - wsolvent());
}
return comp_frac * primary_variables_evaluation_[GTotal];
} else {
comp_frac = comp_frac_[legacyCompIdx];
}
return primary_variables_evaluation_[GTotal];
return comp_frac * primary_variables_evaluation_[GTotal];
} else {
return primary_variables_evaluation_[GTotal] * wellVolumeFractionScaled(comp_idx);
}
@ -605,6 +618,8 @@ namespace Opm
resWell_[0][componentIdx] += resWell_loc.value();
}
assembleControlEq();
// do the local inversion of D.
// we do this manually with invertMatrix to always get our
// specializations in for 3x3 and 4x4 matrices.
@ -618,6 +633,101 @@ namespace Opm
template <typename TypeTag>
void
StandardWell<TypeTag>::
assembleControlEq()
{
EvalWell control_eq(0.0);
switch (well_controls_get_current_type(well_controls_)) {
case THP:
{
OPM_THROW(std::runtime_error, "Not handling THP control for Multisegment wells for now");
// TODO: it will be a function based on BHP <-> THP relation
break;
}
case BHP:
{
const double target_bhp = well_controls_get_current_target(well_controls_);
control_eq = getBhp() - target_bhp;
break;
}
case SURFACE_RATE:
{
int number_phases_under_control = 0;
const double* distr = well_controls_get_current_distr(well_controls_);
for (int phase = 0; phase < number_of_phases_; ++phase) {
if (distr[phase] > 0.0) {
++number_phases_under_control;
}
}
assert(number_phases_under_control > 0);
const double target_rate = well_controls_get_current_target(well_controls_);
if (well_type_ == INJECTOR) {
assert(number_phases_under_control == 1); // only handles single phase injection now
// TODO: considering the solvent part here
control_eq = getGTotal() - target_rate;
} else if (well_type_ == PRODUCER) {
EvalWell rate_for_control(0.);
const EvalWell& g_total = getGTotal();
for (int phase = 0; phase < number_of_phases_; ++phase) {
if (distr[phase] > 0.) {
rate_for_control += g_total * wellVolumeFractionScaled(flowPhaseToEbosCompIdx(phase));
}
}
control_eq = rate_for_control - target_rate;
}
break;
}
case RESERVOIR_RATE:
{
// TODO: repeated code here
int number_phases_under_control = 0;
const double* distr = well_controls_get_current_distr(well_controls_);
for (int phase = 0; phase < number_of_phases_; ++phase) {
if (distr[phase] > 0.0) {
++number_phases_under_control;
}
}
assert(number_phases_under_control > 0);
const double target_rate = well_controls_get_current_target(well_controls_); // reservoir rate target
if (well_type_ == INJECTOR) {
assert(number_phases_under_control == 1);
for (int phase = 0; phase < number_of_phases_; ++phase) {
if (distr[phase] > 0.0) {
control_eq = getGTotal() * scalingFactor(phase) - target_rate;
break;
}
}
} else {
const EvalWell& g_total = getGTotal();
EvalWell rate_for_control(0.0); // reservoir rate
for (int phase = 0; phase < number_of_phases_; ++phase) {
rate_for_control += g_total * wellVolumeFraction( flowPhaseToEbosCompIdx(phase) );
}
control_eq = rate_for_control - target_rate;
}
break;
}
default:
OPM_THROW(std::runtime_error, "Unknown well control control types for well " << name());
}
// using control_eq to update the matrix and residuals
// TODO: we should use a different index system for the well equations
resWell_[0][Bhp] = control_eq.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
invDuneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + numEq);
}
}
template<typename TypeTag>
bool
StandardWell<TypeTag>::
@ -992,9 +1102,13 @@ namespace Opm
well_state.bhp()[well_index] = target;
// TODO: similar to the way below to handle THP
// we should not something related to thp here when there is thp constraint
// or when can calculate the THP (table avaiable or requested for output?)
break;
case THP: {
// TODO: this will be the big task here.
// p_bhp = BHP(THP, rates(p_bhp))
// more sophiscated techniques is required to obtain the bhp and rates here
well_state.thp()[well_index] = target;
const Opm::PhaseUsage& pu = phaseUsage();
@ -1499,9 +1613,7 @@ namespace Opm
StandardWell<TypeTag>::
computeAccumWell()
{
// TODO: it should be num_comp, while it also bring problem for
// the polymer case.
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
for (int eq_idx = 0; eq_idx < num_components_; ++eq_idx) {
F0_[eq_idx] = wellSurfaceVolumeFraction(eq_idx).value();
}
}