Merge pull request #4166 from GitPaean/removing_GTotal

removing the usage of GTotal in the well code
This commit is contained in:
Bård Skaflestad 2022-10-13 23:01:43 +02:00 committed by GitHub
commit dd61e79e86
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 31 additions and 31 deletions

View File

@ -71,7 +71,7 @@ namespace Opm
using typename MSWEval::OffDiagMatrixBlockWellType; using typename MSWEval::OffDiagMatrixBlockWellType;
using MSWEval::GFrac; using MSWEval::GFrac;
using MSWEval::WFrac; using MSWEval::WFrac;
using MSWEval::GTotal; using MSWEval::WQTotal;
using MSWEval::SPres; using MSWEval::SPres;
using MSWEval::numWellEq; using MSWEval::numWellEq;
using typename Base::PressureMatrix; using typename Base::PressureMatrix;

View File

@ -417,14 +417,14 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells,
// update the total rate // TODO: should we have a limitation of the total rate change? // update the total rate // TODO: should we have a limitation of the total rate change?
{ {
primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - relaxation_factor * dwells[seg][GTotal]; primary_variables_[seg][WQTotal] = old_primary_variables[seg][WQTotal] - relaxation_factor * dwells[seg][WQTotal];
// make sure that no injector produce and no producer inject // make sure that no injector produce and no producer inject
if (seg == 0) { if (seg == 0) {
if (baseif_.isInjector()) { if (baseif_.isInjector()) {
primary_variables_[seg][GTotal] = std::max( primary_variables_[seg][GTotal], 0.0); primary_variables_[seg][WQTotal] = std::max( primary_variables_[seg][WQTotal], 0.0);
} else { } else {
primary_variables_[seg][GTotal] = std::min( primary_variables_[seg][GTotal], 0.0); primary_variables_[seg][WQTotal] = std::min( primary_variables_[seg][WQTotal], 0.0);
} }
} }
} }
@ -470,7 +470,7 @@ updatePrimaryVariables(const WellState& well_state) const
total_seg_rate = std::min(total_seg_rate, 0.); total_seg_rate = std::min(total_seg_rate, 0.);
} }
} }
primary_variables_[seg][GTotal] = total_seg_rate; primary_variables_[seg][WQTotal] = total_seg_rate;
if (std::abs(total_seg_rate) > 0.) { if (std::abs(total_seg_rate) > 0.) {
if (has_wfrac_variable) { if (has_wfrac_variable) {
const int water_pos = pu.phase_pos[Water]; const int water_pos = pu.phase_pos[Water];
@ -597,7 +597,7 @@ getSegmentRateUpwinding(const int seg,
const size_t comp_idx) const const size_t comp_idx) const
{ {
const int seg_upwind = upwinding_segments_[seg]; const int seg_upwind = upwinding_segments_[seg];
// the result will contain the derivative with resepct to GTotal in segment seg, // the result will contain the derivative with respect to WQTotal in segment seg,
// and the derivatives with respect to WFrac GFrac in segment seg_upwind. // and the derivatives with respect to WFrac GFrac in segment seg_upwind.
// the derivative with respect to SPres should be zero. // the derivative with respect to SPres should be zero.
if (seg == 0 && baseif_.isInjector()) { if (seg == 0 && baseif_.isInjector()) {
@ -607,23 +607,23 @@ getSegmentRateUpwinding(const int seg,
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
&& Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == comp_idx && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == comp_idx
&& phase == InjectorType::WATER) && phase == InjectorType::WATER)
return primary_variables_evaluation_[seg][GTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
&& Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == comp_idx && Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == comp_idx
&& phase == InjectorType::OIL) && phase == InjectorType::OIL)
return primary_variables_evaluation_[seg][GTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
&& Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == comp_idx && Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == comp_idx
&& phase == InjectorType::GAS) && phase == InjectorType::GAS)
return primary_variables_evaluation_[seg][GTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
return 0.0; return 0.0;
} }
const EvalWell segment_rate = primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg_upwind, comp_idx); const EvalWell segment_rate = primary_variables_evaluation_[seg][WQTotal] * volumeFractionScaled(seg_upwind, comp_idx);
assert(segment_rate.derivative(SPres + Indices::numEq) == 0.); assert(segment_rate.derivative(SPres + Indices::numEq) == 0.);
@ -838,7 +838,7 @@ MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getSegmentRate(const int seg, getSegmentRate(const int seg,
const int comp_idx) const const int comp_idx) const
{ {
return primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg, comp_idx); return primary_variables_evaluation_[seg][WQTotal] * volumeFractionScaled(seg, comp_idx);
} }
template<typename FluidSystem, typename Indices, typename Scalar> template<typename FluidSystem, typename Indices, typename Scalar>
@ -852,9 +852,9 @@ getQs(const int comp_idx) const
template<typename FluidSystem, typename Indices, typename Scalar> template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>:: MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getSegmentGTotal(const int seg) const getSegmentWQTotal(const int seg) const
{ {
return primary_variables_evaluation_[seg][GTotal]; return primary_variables_evaluation_[seg][WQTotal];
} }
template<typename FluidSystem, typename Indices, typename Scalar> template<typename FluidSystem, typename Indices, typename Scalar>
@ -862,7 +862,7 @@ typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>:: MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getWQTotal() const getWQTotal() const
{ {
return getSegmentGTotal(0); return getSegmentWQTotal(0);
} }
template<typename FluidSystem, typename Indices, typename Scalar> template<typename FluidSystem, typename Indices, typename Scalar>
@ -1424,7 +1424,7 @@ handleAccelerationPressureLoss(const int seg,
resWell_[seg][SPres] -= accelerationPressureLoss.value(); resWell_[seg][SPres] -= accelerationPressureLoss.value();
duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq); duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq);
duneD_[seg][seg][SPres][GTotal] -= accelerationPressureLoss.derivative(GTotal + Indices::numEq); duneD_[seg][seg][SPres][WQTotal] -= accelerationPressureLoss.derivative(WQTotal + Indices::numEq);
if (has_wfrac_variable) { if (has_wfrac_variable) {
duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq); duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq);
} }
@ -1463,7 +1463,7 @@ assembleDefaultPressureEq(const int seg,
resWell_[seg][SPres] = pressure_equation.value(); resWell_[seg][SPres] = pressure_equation.value();
const int seg_upwind = upwinding_segments_[seg]; const int seg_upwind = upwinding_segments_[seg];
duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq);
duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + Indices::numEq); duneD_[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq);
if (has_wfrac_variable) { if (has_wfrac_variable) {
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq); duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq);
} }
@ -1533,7 +1533,7 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
} }
// calculate the phase rates based on the primary variables // calculate the phase rates based on the primary variables
const double g_total = primary_variables_[seg][GTotal]; const double g_total = primary_variables_[seg][WQTotal];
for (int p = 0; p < baseif_.numPhases(); ++p) { for (int p = 0; p < baseif_.numPhases(); ++p) {
const double phase_rate = g_total * fractions[p]; const double phase_rate = g_total * fractions[p];
segment_rates[seg*baseif_.numPhases() + p] = phase_rate; segment_rates[seg*baseif_.numPhases() + p] = phase_rate;
@ -1566,8 +1566,8 @@ assembleICDPressureEq(const int seg,
if (const auto& segment = this->segmentSet()[seg]; if (const auto& segment = this->segmentSet()[seg];
(segment.segmentType() == Segment::SegmentType::VALVE) && (segment.segmentType() == Segment::SegmentType::VALVE) &&
(segment.valve().status() == Opm::ICDStatus::SHUT) ) { // we use a zero rate equation to handle SHUT valve (segment.valve().status() == Opm::ICDStatus::SHUT) ) { // we use a zero rate equation to handle SHUT valve
resWell_[seg][SPres] = this->primary_variables_evaluation_[seg][GTotal].value(); resWell_[seg][SPres] = this->primary_variables_evaluation_[seg][WQTotal].value();
duneD_[seg][seg][SPres][GTotal] = 1.; duneD_[seg][seg][SPres][WQTotal] = 1.;
auto& ws = well_state.well(baseif_.indexOfWell()); auto& ws = well_state.well(baseif_.indexOfWell());
ws.segments.pressure_drop_friction[seg] = 0.; ws.segments.pressure_drop_friction[seg] = 0.;
@ -1603,7 +1603,7 @@ assembleICDPressureEq(const int seg,
const int seg_upwind = upwinding_segments_[seg]; const int seg_upwind = upwinding_segments_[seg];
resWell_[seg][SPres] = pressure_equation.value(); resWell_[seg][SPres] = pressure_equation.value();
duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq);
duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + Indices::numEq); duneD_[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq); duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq);
} }
@ -1799,14 +1799,14 @@ updateUpwindingSegments()
// special treatment is needed for segment 0 // special treatment is needed for segment 0
if (seg == 0) { if (seg == 0) {
// we are not supposed to have injecting producers and producing injectors // we are not supposed to have injecting producers and producing injectors
assert( ! (baseif_.isProducer() && primary_variables_evaluation_[seg][GTotal] > 0.) ); assert( ! (baseif_.isProducer() && primary_variables_evaluation_[seg][WQTotal] > 0.) );
assert( ! (baseif_.isInjector() && primary_variables_evaluation_[seg][GTotal] < 0.) ); assert( ! (baseif_.isInjector() && primary_variables_evaluation_[seg][WQTotal] < 0.) );
upwinding_segments_[seg] = seg; upwinding_segments_[seg] = seg;
continue; continue;
} }
// for other normal segments // for other normal segments
if (primary_variables_evaluation_[seg][GTotal] <= 0.) { if (primary_variables_evaluation_[seg][WQTotal] <= 0.) {
upwinding_segments_[seg] = seg; upwinding_segments_[seg] = seg;
} else { } else {
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());

View File

@ -68,7 +68,7 @@ protected:
// Table showing the primary variable indices, depending on what phases are present: // Table showing the primary variable indices, depending on what phases are present:
// //
// WOG OG WG WO W/O/G (single phase) // WOG OG WG WO W/O/G (single phase)
// GTotal 0 0 0 0 0 // WQTotal 0 0 0 0 0
// WFrac 1 -1000 1 1 -1000 // WFrac 1 -1000 1 1 -1000
// GFrac 2 1 -1000 -1000 -1000 // GFrac 2 1 -1000 -1000 -1000
// Spres 3 2 2 2 1 // Spres 3 2 2 2 1
@ -83,7 +83,7 @@ protected:
static constexpr bool has_wfrac_variable = has_water && Indices::numPhases > 1; static constexpr bool has_wfrac_variable = has_water && Indices::numPhases > 1;
static constexpr bool has_gfrac_variable = has_gas && has_oil; static constexpr bool has_gfrac_variable = has_gas && has_oil;
static constexpr int GTotal = 0; static constexpr int WQTotal = 0;
static constexpr int WFrac = has_wfrac_variable ? 1 : -1000; static constexpr int WFrac = has_wfrac_variable ? 1 : -1000;
static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000; static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000;
static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1; static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1;
@ -190,7 +190,7 @@ protected:
EvalWell getFrictionPressureLoss(const int seg) const; EvalWell getFrictionPressureLoss(const int seg) const;
EvalWell getHydroPressureLoss(const int seg) const; EvalWell getHydroPressureLoss(const int seg) const;
EvalWell getQs(const int comp_idx) const; EvalWell getQs(const int comp_idx) const;
EvalWell getSegmentGTotal(const int seg) const; EvalWell getSegmentWQTotal(const int seg) const;
EvalWell getSegmentPressure(const int seg) const; EvalWell getSegmentPressure(const int seg) const;
EvalWell getSegmentRate(const int seg, const int comp_idx) const; EvalWell getSegmentRate(const int seg, const int comp_idx) const;
EvalWell getSegmentRateUpwinding(const int seg, EvalWell getSegmentRateUpwinding(const int seg,

View File

@ -1655,10 +1655,10 @@ namespace Opm
const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) * this->well_efficiency_factor_; const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) * this->well_efficiency_factor_;
const int seg_upwind = this->upwinding_segments_[seg]; const int seg_upwind = this->upwinding_segments_[seg];
// segment_rate contains the derivatives with respect to GTotal in seg, // segment_rate contains the derivatives with respect to WQTotal in seg,
// and WFrac and GFrac in seg_upwind // and WFrac and GFrac in seg_upwind
this->resWell_[seg][comp_idx] -= segment_rate.value(); this->resWell_[seg][comp_idx] -= segment_rate.value();
this->duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + Indices::numEq); this->duneD_[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq); this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq);
} }
@ -1676,10 +1676,10 @@ namespace Opm
const EvalWell inlet_rate = this->getSegmentRateUpwinding(inlet, comp_idx) * this->well_efficiency_factor_; const EvalWell inlet_rate = this->getSegmentRateUpwinding(inlet, comp_idx) * this->well_efficiency_factor_;
const int inlet_upwind = this->upwinding_segments_[inlet]; const int inlet_upwind = this->upwinding_segments_[inlet];
// inlet_rate contains the derivatives with respect to GTotal in inlet, // inlet_rate contains the derivatives with respect to WQTotal in inlet,
// and WFrac and GFrac in inlet_upwind // and WFrac and GFrac in inlet_upwind
this->resWell_[seg][comp_idx] += inlet_rate.value(); this->resWell_[seg][comp_idx] += inlet_rate.value();
this->duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + Indices::numEq); this->duneD_[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq); this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq);
} }

View File

@ -71,7 +71,7 @@ protected:
// Table showing the primary variable indices, depending on what phases are present: // Table showing the primary variable indices, depending on what phases are present:
// //
// WOG OG WG WO W/O/G (single phase) // WOG OG WG WO W/O/G (single phase)
// GTotal 0 0 0 0 0 // WQTotal 0 0 0 0 0
// WFrac 1 -1000 -1000 1 -1000 // WFrac 1 -1000 -1000 1 -1000
// GFrac 2 1 1 -1000 -1000 // GFrac 2 1 1 -1000 -1000
// Spres 3 2 2 2 1 // Spres 3 2 2 2 1