removing well_soltutions_ from WellStateFullyImplicitBlackoilDense

adding function setWellSolutions() to StandardWell.

The class WellStateFullyImplicitBlackoilDense is ready to be removed
now, while the only thing can go wrong compred with the original version
is the group control, which is not tested yet.
This commit is contained in:
Kai Bao 2017-08-07 17:49:35 +02:00
parent 8441eb77bd
commit 5af15fa63f
7 changed files with 117 additions and 197 deletions

View File

@ -294,9 +294,6 @@ void wellsToState( const data::Wells& wells,
{
// Set base class variables.
wellsToState(wells, phases, static_cast<WellStateFullyImplicitBlackoil&>(state));
// Set wellSolution() variable.
state.setWellSolutions(phases);
}

View File

@ -104,7 +104,7 @@ namespace Opm
const int num_cells);
virtual void setWellVariables(const WellState& well_state);
virtual void setWellVariables();
// TODO: to check whether all the paramters are required
void computePerfRate(const IntensiveQuantities& intQuants,
@ -158,6 +158,9 @@ namespace Opm
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials) const;
virtual void setWellSolutions(const WellState& well_state) const;
protected:
// protected functions from the Base class

View File

@ -101,9 +101,8 @@ namespace Opm
template<typename TypeTag>
void StandardWell<TypeTag>::
setWellVariables(const WellState& well_state)
setWellVariables()
{
const int nw = well_state.bhp().size();
// TODO: using numComp here is only to make the 2p + dummy phase work
// TODO: in theory, we should use numWellEq here.
// for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) {
@ -1260,89 +1259,7 @@ namespace Opm
break;
} // end of switch
std::vector<double> g = {1.0, 1.0, 0.01};
if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) {
for (int phase = 0; phase < np; ++phase) {
g[phase] = distr[phase];
}
}
// the number of wells
const int nw = xw.bhp().size();
switch (well_controls_iget_type(wc, current)) {
case THP:
case BHP: {
well_solutions_[XvarWell] = 0.0;
if (well_type_ == INJECTOR) {
for (int p = 0; p < np; ++p) {
well_solutions_[XvarWell] += xw.wellRates()[np*well_index + p] * comp_frac_[p];
}
} else {
for (int p = 0; p < np; ++p) {
well_solutions_[XvarWell] += g[p] * xw.wellRates()[np*well_index + p];
}
}
break;
}
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
well_solutions_[XvarWell] = xw.bhp()[well_index];
break;
} // end of switch
double tot_well_rate = 0.0;
for (int p = 0; p < np; ++p) {
tot_well_rate += g[p] * xw.wellRates()[np*well_index + p];
}
if(std::abs(tot_well_rate) > 0) {
if (active()[ Water ]) {
well_solutions_[WFrac] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate;
}
if (active()[ Gas ]) {
well_solutions_[GFrac] = g[Gas] * (xw.wellRates()[np*well_index + Gas] - xw.solventWellRate(well_index)) / tot_well_rate ;
}
if (has_solvent) {
well_solutions_[SFrac] = g[Gas] * xw.solventWellRate(well_index) / tot_well_rate ;
}
} else { // tot_well_rate == 0
if (well_type_ == INJECTOR) {
// only single phase injection handled
if (active()[Water]) {
if (distr[Water] > 0.0) {
well_solutions_[WFrac] = 1.0;
} else {
well_solutions_[WFrac] = 0.0;
}
}
if (active()[Gas]) {
if (distr[Gas] > 0.0) {
well_solutions_[GFrac] = 1.0 - wsolvent();
if (has_solvent) {
well_solutions_[SFrac] = wsolvent();
}
} else {
well_solutions_[GFrac] = 0.0;
}
}
// TODO: it is possible to leave injector as a oil well,
// when F_w and F_g both equals to zero, not sure under what kind of circumstance
// this will happen.
} else if (well_type_ == PRODUCER) { // producers
// TODO: the following are not addressed for the solvent case yet
if (active()[Water]) {
well_solutions_[WFrac] = 1.0 / np;
}
if (active()[Gas]) {
well_solutions_[GFrac] = 1.0 / np;
}
} else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
}
setWellSolutions(xw);
}
@ -2127,4 +2044,102 @@ namespace Opm
}
}
template<typename TypeTag>
void
StandardWell<TypeTag>::
setWellSolutions(const WellState& well_state) const
{
const int np = number_of_phases_;
const int well_index = index_of_well_;
const WellControls* wc = well_controls_;
const double* distr = well_controls_get_current_distr(wc);
std::vector<double> g = {1.0, 1.0, 0.01};
if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
for (int phase = 0; phase < np; ++phase) {
g[phase] = distr[phase];
}
}
// the number of wells
const int nw = well_state.bhp().size();
switch (well_controls_get_current_type(wc)) {
case THP:
case BHP: {
well_solutions_[XvarWell] = 0.0;
if (well_type_ == INJECTOR) {
for (int p = 0; p < np; ++p) {
well_solutions_[XvarWell] += well_state.wellRates()[np*well_index + p] * comp_frac_[p];
}
} else {
for (int p = 0; p < np; ++p) {
well_solutions_[XvarWell] += g[p] * well_state.wellRates()[np*well_index + p];
}
}
break;
}
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
well_solutions_[XvarWell] = well_state.bhp()[well_index];
break;
} // end of switch
double tot_well_rate = 0.0;
for (int p = 0; p < np; ++p) {
tot_well_rate += g[p] * well_state.wellRates()[np*well_index + p];
}
if(std::abs(tot_well_rate) > 0) {
if (active()[ Water ]) {
well_solutions_[WFrac] = g[Water] * well_state.wellRates()[np*well_index + Water] / tot_well_rate;
}
if (active()[ Gas ]) {
well_solutions_[GFrac] = g[Gas] * (well_state.wellRates()[np*well_index + Gas] - well_state.solventWellRate(well_index)) / tot_well_rate ;
}
if (has_solvent) {
well_solutions_[SFrac] = g[Gas] * well_state.solventWellRate(well_index) / tot_well_rate ;
}
} else { // tot_well_rate == 0
if (well_type_ == INJECTOR) {
// only single phase injection handled
if (active()[Water]) {
if (distr[Water] > 0.0) {
well_solutions_[WFrac] = 1.0;
} else {
well_solutions_[WFrac] = 0.0;
}
}
if (active()[Gas]) {
if (distr[Gas] > 0.0) {
well_solutions_[GFrac] = 1.0 - wsolvent();
if (has_solvent) {
well_solutions_[SFrac] = wsolvent();
}
} else {
well_solutions_[GFrac] = 0.0;
}
}
// TODO: it is possible to leave injector as a oil well,
// when F_w and F_g both equals to zero, not sure under what kind of circumstance
// this will happen.
} else if (well_type_ == PRODUCER) { // producers
// TODO: the following are not addressed for the solvent case yet
if (active()[Water]) {
well_solutions_[WFrac] = 1.0 / np;
}
if (active()[Gas]) {
well_solutions_[GFrac] = 1.0 / np;
}
} else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
}
}
}

View File

@ -172,7 +172,7 @@ namespace Opm {
/// return true if wells are available on this process
bool localWellsActive() const;
void setWellVariables(const WellState& xw);
void setWellVariables();
void computeAccumWells();

View File

@ -174,7 +174,7 @@ namespace Opm {
updateWellControls(well_state);
updateGroupControls(well_state);
// Set the primary variables for the wells
setWellVariables(well_state);
setWellVariables();
if (iterationIdx == 0) {
computeWellConnectionPressures(ebosSimulator, well_state);
@ -430,10 +430,10 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
setWellVariables(const WellState& xw)
setWellVariables()
{
for (auto& well : well_container_) {
well->setWellVariables(xw);
well->setWellVariables();
}
}
@ -499,7 +499,7 @@ namespace Opm {
{
updateWellControls(well_state);
updateGroupControls(well_state);
setWellVariables(well_state);
setWellVariables();
}
} while (it < 15);
@ -758,7 +758,11 @@ namespace Opm {
if (well_collection_->requireWellPotentials()) {
// calculate the well potentials
setWellVariables(well_state);
// TODO: for the purpose of group control, not tested yet
for (const auto& well : well_container_) {
well->setWellSolutions(well_state);
}
setWellVariables();
computeWellConnectionPressures(ebos_simulator, well_state);
// To store well potentials for each well

View File

@ -104,7 +104,7 @@ namespace Opm
const double gravity_arg,
const int num_cells);
virtual void setWellVariables(const WellState& well_state) = 0;
virtual void setWellVariables() = 0;
virtual bool getWellConvergence(Simulator& ebosSimulator,
const std::vector<double>& B_avg,
@ -153,6 +153,8 @@ namespace Opm
virtual void updateWellControl(WellState& xw) const = 0;
virtual void setWellSolutions(const WellState& well_state) const = 0;
protected:
const std::vector<bool>& active() const;

View File

@ -108,101 +108,6 @@ namespace Opm
}
}
}
// TODO: the reason to keep this is to avoid getting defaulted value BHP
// some facilities needed from opm-parser or opm-core
// It is a little tricky, since sometimes before applying group control, the only
// available constraints in the well_controls is the defaulted BHP value, and it
// is really not desirable to use this value to enter the Newton iterations.
setWellSolutions(pu);
}
/// Set wellSolutions() based on the base class members.
void setWellSolutions(const PhaseUsage& pu)
{
// Set nw and np, or return if no wells.
if (wells_.get() == nullptr) {
return;
}
const int nw = wells_->number_of_wells;
if (nw == 0) {
return;
}
const int np = wells_->number_of_phases;
const int numComp = pu.has_solvent? np+1:np;
well_solutions_.clear();
well_solutions_.resize(nw * numComp, 0.0);
std::vector<double> g = {1.0,1.0,0.01};
for (int w = 0; w < nw; ++w) {
WellControls* wc = wells_->ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
const int current = currentControls()[w];
well_controls_set_current( wc, current);
const WellType& well_type = wells_->type[w];
switch (well_controls_iget_type(wc, current)) {
case THP: // Intentional fall-through
case BHP:
if (well_type == INJECTOR) {
for (int p = 0; p < np; ++p) {
well_solutions_[w] += wellRates()[np*w + p] * wells_->comp_frac[np*w + p];
}
} else {
for (int p = 0; p < np; ++p) {
well_solutions_[w] += g[p] * wellRates()[np*w + p];
}
}
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
wellSolutions()[w] = bhp()[w];
break;
}
double total_rates = 0.0;
for (int p = 0; p < np; ++p) {
total_rates += g[p] * wellRates()[np*w + p];
}
const int waterpos = pu.phase_pos[Water];
const int gaspos = pu.phase_pos[Gas];
assert(np > 2 || (np == 2 && !pu.phase_used[Gas]));
// assumes the gas fractions are stored after water fractions
if(std::abs(total_rates) > 0) {
if( pu.phase_used[Water] ) {
wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + waterpos] / total_rates;
}
if( pu.phase_used[Gas] ) {
wellSolutions()[2*nw + w] = g[Gas] * (wellRates()[np*w + gaspos] - solventWellRate(w)) / total_rates ;
}
if( pu.has_solvent) {
wellSolutions()[3*nw + w] = g[Gas] * solventWellRate(w) / total_rates;
}
} else {
if( pu.phase_used[Water] ) {
wellSolutions()[nw + w] = wells_->comp_frac[np*w + waterpos];
}
if( pu.phase_used[Gas] ) {
wellSolutions()[2*nw + w] = wells_->comp_frac[np*w + gaspos];
}
if (pu.has_solvent) {
wellSolutions()[3*nw + w] = 0;
}
}
}
}
@ -212,11 +117,6 @@ namespace Opm
init(wells, state.pressure(), dummy_state, pu) ;
}
/// One rate per phase and well connection.
std::vector<double>& wellSolutions() { return well_solutions_; }
const std::vector<double>& wellSolutions() const { return well_solutions_; }
/// One rate pr well connection.
std::vector<double>& perfRateSolvent() { return perfRateSolvent_; }
const std::vector<double>& perfRateSolvent() const { return perfRateSolvent_; }
@ -252,7 +152,6 @@ namespace Opm
private:
std::vector<double> well_solutions_;
std::vector<double> perfRateSolvent_;
};