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

suppport for gas water multisegment well

Browse Source
This commit is contained in:
Kai Bao 2024-12-19 15:38:18 +01:00
parent cbc12a915b
commit cf1b52e52f
8 changed files with 105 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
opm/simulators/wells/MultisegmentWellAssemble.cpp
View File

@ -279,9 +279,7 @@ assemblePressureEq(const int seg,
const int outlet_segment_index,
const EvalWell& pressure_equation,
const EvalWell& outlet_pressure,
Equations& eqns1,
bool wfrac,
bool gfrac) const
Equations& eqns1) const
{
/*
This method does *not* need communication.
@ -290,10 +288,10 @@ assemblePressureEq(const int seg,
eqns.residual()[seg][SPres] += pressure_equation.value();
eqns.D()[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq);
eqns.D()[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq);
if (wfrac) {
if (has_wfrac_variable) {
eqns.D()[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq);
}
if (gfrac) {
if (has_gfrac_variable) {
eqns.D()[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq);
}
@ -355,10 +353,10 @@ assembleOutflowTerm(const int seg,
MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
eqns.residual()[seg][comp_idx] -= segment_rate.value();
eqns.D()[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if constexpr (has_wfrac_variable) {
eqns.D()[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if constexpr (has_gfrac_variable) {
eqns.D()[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq);
}
// pressure derivative should be zero
@ -380,10 +378,10 @@ assembleInflowTerm(const int seg,
MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
eqns.residual()[seg][comp_idx] += inlet_rate.value();
eqns.D()[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if constexpr (has_wfrac_variable) {
eqns.D()[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if constexpr (has_gfrac_variable) {
eqns.D()[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq);
}
// pressure derivative should be zero

16
opm/simulators/wells/MultisegmentWellAssemble.hpp
View File

@ -42,17 +42,9 @@ template<class Scalar> class WellState;
template<class FluidSystem, class Indices>
class MultisegmentWellAssemble
{
static constexpr bool has_water = (Indices::waterSwitchIdx >= 0);
static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0);
static constexpr bool has_oil = (Indices::numPhases - has_gas - has_water) > 0;
// In the implementation, one should use has_wfrac_variable
// rather than has_water to check if you should do something
// with the variable at the WFrac location, similar for GFrac.
static constexpr bool has_wfrac_variable = has_water && Indices::numPhases > 1;
static constexpr bool has_gfrac_variable = has_gas && has_oil;
static constexpr int WQTotal = 0;
static constexpr bool has_wfrac_variable = Indices::waterEnabled && Indices::oilEnabled;
static constexpr bool has_gfrac_variable = Indices::gasEnabled && Indices::numPhases > 1;
static constexpr int WFrac = 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;
@ -107,9 +99,7 @@ public:
const int outlet_segment_index,
const EvalWell& pressure_equation,
const EvalWell& outlet_pressure,
Equations& eqns,
bool wfrac = has_wfrac_variable,
bool gfrac = has_gfrac_variable) const;
Equations& eqns) const;
//! \brief Assembles a trivial equation.
void assembleTrivialEq(const int seg,

4
opm/simulators/wells/MultisegmentWellEval.cpp
View File

@ -368,9 +368,7 @@ assembleICDPressureEq(const int seg,
MultisegmentWellAssemble(baseif_).
assemblePressureEq(seg, seg_upwind, outlet_segment_index,
pressure_equation, outlet_pressure,
linSys_,
FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
linSys_);
assembleAccelerationAndHydroPressureLosses(seg, well_state, use_average_density);
}

115
opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
View File

@ -225,7 +225,6 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
const auto pvtReg = std::max(well_.wellEcl().pvt_table_number() - 1, 0);
const PhaseUsage& pu = well_.phaseUsage();
assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
const int oil_pos = pu.phase_pos[Oil];
auto& ws = well_state.well(well_.indexOfWell());
@ -236,18 +235,34 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
auto& segment_pressure = segments.pressure;
for (std::size_t seg = 0; seg < value_.size(); ++seg) {
std::vector<Scalar> fractions(well_.numPhases(), 0.0);
fractions[oil_pos] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water];
fractions[water_pos] = value_[seg][WFrac];
fractions[oil_pos] -= fractions[water_pos];
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
fractions[oil_pos] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water];
fractions[water_pos] = value_[seg][WFrac];
fractions[oil_pos] -= fractions[water_pos];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
fractions[gas_pos] = value_[seg][GFrac];
fractions[oil_pos] -= fractions[gas_pos];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water];
fractions[water_pos] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
fractions[gas_pos] = value_[seg][GFrac];
fractions[water_pos] -= fractions[gas_pos];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
fractions[gas_pos] = value_[seg][GFrac];
fractions[oil_pos] -= fractions[gas_pos];
fractions[gas_pos] = 1.0;
}
// convert the fractions to be Q_p / G_total to calculate the phase rates
@ -287,7 +302,7 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
Scalar rsMax = 0.0;
Scalar rvMax = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
// Both oil and gas active.
rsMax = FluidSystem::oilPvt()
.saturatedGasDissolutionFactor(pvtReg, temperature, segment_pressure[seg]);
@ -300,7 +315,7 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
const auto& [Rs, Rv] = well_.rateConverter().inferDissolvedVaporisedRatio
(rsMax, rvMax, segment_rates.begin() + seg * well_.numPhases());
if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (! (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))) ) {
vapoil[seg] = disgas[seg] = 0.0;
}
else {
@ -355,8 +370,10 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
}
// 5) Local condition phase viscosities.
segments.phase_viscosity[seg * well_.numPhases() + pu.phase_pos[Oil]] =
FluidSystem::oilPvt().viscosity(pvtReg, temperature, segment_pressure[seg], Rs);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
segments.phase_viscosity[seg * well_.numPhases() + pu.phase_pos[Oil]] =
FluidSystem::oilPvt().viscosity(pvtReg, temperature, segment_pressure[seg], Rs);
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
segments.phase_viscosity[seg * well_.numPhases() + pu.phase_pos[Water]] =
@ -424,59 +441,71 @@ processFractions(const int seg)
std::vector<Scalar> fractions(well_.numPhases(), 0.0);
assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
const int oil_pos = pu.phase_pos[Oil];
fractions[oil_pos] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
fractions[pu.phase_pos[Oil]] = 1.0;
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
const int water_pos = pu.phase_pos[Water];
fractions[water_pos] = value_[seg][WFrac];
fractions[oil_pos] -= fractions[water_pos];
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
fractions[pu.phase_pos[Water]] = value_[seg][WFrac];
fractions[pu.phase_pos[Oil]] -= fractions[pu.phase_pos[Water]];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
fractions[pu.phase_pos[Gas]] = value_[seg][GFrac];
fractions[pu.phase_pos[Oil]] -= fractions[pu.phase_pos[Gas]];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
fractions[pu.phase_pos[Water]] = 1.0;
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
const int gas_pos = pu.phase_pos[Gas];
fractions[gas_pos] = value_[seg][GFrac];
fractions[oil_pos] -= fractions[gas_pos];
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
fractions[pu.phase_pos[Gas]] = value_[seg][GFrac];
fractions[pu.phase_pos[Water]] -= fractions[pu.phase_pos[Gas]];
}
} else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
fractions[pu.phase_pos[Gas]] = 1.0;
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water];
if (fractions[water_pos] < 0.0) {
if (fractions[pu.phase_pos[Water]] < 0.0) {
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]);
fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[pu.phase_pos[Water]]);
}
fractions[oil_pos] /= (1.0 - fractions[water_pos]);
fractions[water_pos] = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
fractions[pu.phase_pos[Oil]] /= (1.0 - fractions[pu.phase_pos[Water]]);
}
fractions[pu.phase_pos[Water]] = 0.0;
}
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
if (fractions[gas_pos] < 0.0) {
if (fractions[pu.phase_pos[Gas]] < 0.0) {
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
fractions[pu.phase_pos[Water]] /= (1.0 - fractions[gas_pos]);
fractions[pu.phase_pos[Water]] /= (1.0 - fractions[pu.phase_pos[Gas]]);
}
fractions[oil_pos] /= (1.0 - fractions[gas_pos]);
fractions[gas_pos] = 0.0;
if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ) {
fractions[pu.phase_pos[Oil]] /= (1.0 - fractions[pu.phase_pos[Gas]]);
}
fractions[pu.phase_pos[Gas]] = 0.0;
}
}
if (fractions[oil_pos] < 0.0) {
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
fractions[pu.phase_pos[Water]] /= (1.0 - fractions[oil_pos]);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
if (fractions[pu.phase_pos[Oil]] < 0.0) {
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
fractions[pu.phase_pos[Water]] /= (1.0 - fractions[pu.phase_pos[Oil]]);
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[pu.phase_pos[Oil]]);
}
fractions[pu.phase_pos[Oil]] = 0.0;
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]);
}
fractions[oil_pos] = 0.0;
}
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
if constexpr (has_wfrac_variable) {
value_[seg][WFrac] = fractions[pu.phase_pos[Water]];
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
if constexpr (has_gfrac_variable) {
value_[seg][GFrac] = fractions[pu.phase_pos[Gas]];
}
}

15
opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp
View File

@ -53,19 +53,12 @@ public:
//
// WOG OG WG WO W/O/G (single phase)
// WQTotal 0 0 0 0 0
// WFrac 1 -1000 1 1 -1000
// GFrac 2 1 -1000 -1000 -1000
// WFrac 1 -1000 -1000 1 -1000
// GFrac 2 1 1 -1000 -1000
// Spres 3 2 2 2 1
static constexpr bool has_water = (Indices::waterSwitchIdx >= 0);
static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0);
static constexpr bool has_oil = (Indices::numPhases - has_gas - has_water) > 0;
// In the implementation, one should use has_wfrac_variable
// rather than has_water to check if you should do something
// with the variable at the WFrac location, similar for GFrac.
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_wfrac_variable = Indices::waterEnabled && Indices::oilEnabled;
static constexpr bool has_gfrac_variable = Indices::gasEnabled && Indices::numPhases > 1;
static constexpr int WQTotal = 0;
static constexpr int WFrac = has_wfrac_variable ? 1 : -1000;

34
opm/simulators/wells/MultisegmentWellSegments.cpp
View File

@ -514,7 +514,7 @@ getFrictionPressureLoss(const int seg,
// at segment node while fraction derivatives are given at upwind node.
if (seg != seg_upwind) {
if (!extra_reverse_flow_derivatives){
if (!extra_reverse_flow_derivatives) {
constexpr int WQTotal = Indices::numEq + PrimaryVariables::WQTotal;
constexpr int SPres = Indices::numEq + PrimaryVariables::SPres;
density.setDerivative(WQTotal, 0.0);
@ -522,12 +522,12 @@ getFrictionPressureLoss(const int seg,
visc.setDerivative(WQTotal, 0.0);
visc.setDerivative(SPres, 0.0);
} else {
if (PrimaryVariables::has_water){
if (PrimaryVariables::has_wfrac_variable) {
constexpr int WFrac = Indices::numEq + PrimaryVariables::WFrac;
density.setDerivative(WFrac, 0.0);
visc.setDerivative(WFrac, 0.0);
}
if (PrimaryVariables::has_gas){
if (PrimaryVariables::has_gfrac_variable) {
constexpr int GFrac = Indices::numEq + PrimaryVariables::GFrac;
density.setDerivative(GFrac, 0.0);
visc.setDerivative(GFrac, 0.0);
@ -600,10 +600,10 @@ pressureDropSpiralICD(const int seg,
zero_mask[PrimaryVariables::WQTotal] = true;
zero_mask[PrimaryVariables::SPres] = true;
} else {
if constexpr (PrimaryVariables::has_water){
if constexpr (PrimaryVariables::has_wfrac_variable) {
zero_mask[PrimaryVariables::WFrac] = true;
}
if constexpr (PrimaryVariables::has_gas){
if constexpr (PrimaryVariables::has_gfrac_variable) {
zero_mask[PrimaryVariables::GFrac] = true;
}
// mass_rate has no extra derivatives (they are organized as in equations)
@ -712,10 +712,10 @@ pressureDropAutoICD(const int seg,
zero_mask[PrimaryVariables::WQTotal] = true;
zero_mask[PrimaryVariables::SPres] = true;
} else {
if (PrimaryVariables::has_water){
if (PrimaryVariables::has_wfrac_variable) {
zero_mask[PrimaryVariables::WFrac] = true;
}
if (PrimaryVariables::has_gas){
if (PrimaryVariables::has_gfrac_variable) {
zero_mask[PrimaryVariables::GFrac] = true;
}
// mass_rate has no extra derivatives (they are organized as in equations)
@ -787,7 +787,7 @@ pressureDropValve(const int seg,
// For reference: the pressure equation assumes pressure/flow derivatives are given
// at segment node while fraction derivatives are given at upwind node.
if (seg != seg_upwind) {
if (!extra_reverse_flow_derivatives){
if (!extra_reverse_flow_derivatives) {
constexpr int WQTotal = Indices::numEq + PrimaryVariables::WQTotal;
constexpr int SPres = Indices::numEq + PrimaryVariables::SPres;
density.setDerivative(WQTotal, 0.0);
@ -795,12 +795,12 @@ pressureDropValve(const int seg,
visc.setDerivative(WQTotal, 0.0);
visc.setDerivative(SPres, 0.0);
} else {
if (PrimaryVariables::has_water){
if (PrimaryVariables::has_wfrac_variable) {
constexpr int WFrac = Indices::numEq + PrimaryVariables::WFrac;
density.setDerivative(WFrac, 0.0);
visc.setDerivative(WFrac, 0.0);
}
if (PrimaryVariables::has_gas){
if (PrimaryVariables::has_gfrac_variable) {
constexpr int GFrac = Indices::numEq + PrimaryVariables::GFrac;
density.setDerivative(GFrac, 0.0);
visc.setDerivative(GFrac, 0.0);
@ -841,17 +841,17 @@ accelerationPressureLossContribution(const int seg,
const int seg_upwind = upwinding_segments_[seg];
EvalWell density = densities_[seg_upwind];
if (seg != seg_upwind) {
if (!extra_reverse_flow_derivatives){
if (!extra_reverse_flow_derivatives) {
constexpr int WQTotal = Indices::numEq + PrimaryVariables::WQTotal;
constexpr int SPres = Indices::numEq + PrimaryVariables::SPres;
density.setDerivative(WQTotal, 0.0);
density.setDerivative(SPres, 0.0);
} else {
if (PrimaryVariables::has_water){
if (PrimaryVariables::has_wfrac_variable) {
constexpr int WFrac = Indices::numEq + PrimaryVariables::WFrac;
density.setDerivative(WFrac, 0.0);
}
if (PrimaryVariables::has_gas){
if (PrimaryVariables::has_gfrac_variable) {
constexpr int GFrac = Indices::numEq + PrimaryVariables::GFrac;
density.setDerivative(GFrac, 0.0);
}
@ -933,10 +933,10 @@ mixtureDensity(const int seg) const
continue;
}
const auto compIdx = Indices::
canonicalToActiveComponentIndex(phIdx);
const auto active_comp_index = Indices::
canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phIdx));
mixDens += q[compIdx].value() * rho[compIdx].value();
mixDens += q[active_comp_index].value() * rho[active_comp_index].value();
}
return mixDens;
@ -977,7 +977,7 @@ mixtureDensityWithExponents(const AutoICD& aicd, const int seg) const
for (const auto& [fsPhaseIdx, densityExponent] : densityExponents) {
if (FluidSystem::phaseIsActive(fsPhaseIdx)) {
const auto compIdx = Indices::
canonicalToActiveComponentIndex(fsPhaseIdx);
canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(fsPhaseIdx));
// exp = (aicd.*densityExponent)() in native syntax.
const auto exp = std::invoke(densityExponent, aicd);

3
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@ -103,9 +103,6 @@ namespace Opm
"dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet."
" \n See (WCONINJE item 10 / WCONHIST item 8)");
}
if constexpr (!Indices::oilEnabled && Indices::numPhases > 1) {
OPM_THROW(std::runtime_error, "water + gas case not supported by multisegment well yet");
}
this->thp_update_iterations = true;
}

2
opm/simulators/wells/StandardWellPrimaryVariables.cpp
View File

@ -382,7 +382,7 @@ copyToWellState(WellState<Scalar>& well_state,
if (scal > 0) {
F[p] /= scal ;
} else {
// this should only happens to injection wells
// this should only happen to injection wells
F[p] = 0.;
}
}