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

Merge pull request #4640 from steink/ms_pressure_eq_derivatives

Browse Source 
Update of multisegment well pressure equations - testing
This commit is contained in:
Kai Bao
2023-05-26 22:57:23 +02:00
committed by GitHub
parent 84cb77a870 3beef758e4
commit 23c0126095
9 changed files with 136 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
opm/simulators/flow/BlackoilModelParametersEbos.hpp
View File

@@ -164,6 +164,10 @@ template<class TypeTag, class MyTypeTag>
struct MaximumNumberOfWellSwitches {
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct UseAverageDensityMsWells {
using type = UndefinedProperty;
};
template<class TypeTag>
struct DbhpMaxRel<TypeTag, TTag::FlowModelParameters> {
@@ -308,6 +312,14 @@ template<class TypeTag>
struct MaximumNumberOfWellSwitches<TypeTag, TTag::FlowModelParameters> {
static constexpr int value = 3;
};
template<class TypeTag>
struct UseAverageDensityMsWells<TypeTag, TTag::FlowModelParameters> {
static constexpr bool value = false;
};
// if openMP is available, determine the number threads per process automatically.
#if _OPENMP
@@ -421,6 +433,9 @@ namespace Opm
/// Maximum number of times a well can switch to the same controt
int max_number_of_well_switches_;
/// Whether to approximate segment densities by averaging over segment and its outlet
bool use_average_density_ms_wells_;
/// Construct from user parameters or defaults.
@@ -457,6 +472,7 @@ namespace Opm
check_well_operability_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheck);
check_well_operability_iter_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter);
max_number_of_well_switches_ = EWOMS_GET_PARAM(TypeTag, int, MaximumNumberOfWellSwitches);
use_average_density_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseAverageDensityMsWells);
deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
}
@@ -495,6 +511,7 @@ namespace Opm
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheck, "Enable the well operability checking");
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter, "Enable the well operability checking during iterations");
EWOMS_REGISTER_PARAM(TypeTag, int, MaximumNumberOfWellSwitches, "Maximum number of times a well can switch to the same control");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseAverageDensityMsWells, "Approximate segment densitities by averaging over segment and its outlet");
}
};
} // namespace Opm

3
opm/simulators/flow/FlowMainEbos.hpp
View File

@@ -220,6 +220,9 @@ void handleExtraConvergenceOutput(SimulatorReport& report,
EWOMS_HIDE_PARAM(TypeTag, VtkWriteTortuosities);
EWOMS_HIDE_PARAM(TypeTag, VtkWriteDiffusionCoefficients);
EWOMS_HIDE_PARAM(TypeTag, VtkWriteEffectiveDiffusionCoefficients);
// hide average density option
EWOMS_HIDE_PARAM(TypeTag, UseAverageDensityMsWells);
EWOMS_END_PARAM_REGISTRATION(TypeTag);

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

@@ -216,6 +216,36 @@ assemblePressureLoss(const int seg,
}
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
assembleHydroPressureLoss(const int seg,
const int seg_density,
const EvalWell& hydro_pressure_drop_seg,
Equations& eqns1) const
{
MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
eqns.residual()[seg][SPres] -= hydro_pressure_drop_seg.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
eqns.D()[seg][seg_density][SPres][pv_idx] -= hydro_pressure_drop_seg.derivative(pv_idx + Indices::numEq);
}
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
assemblePressureEqExtraDerivatives(const int seg,
const int seg_upwind,
const EvalWell& extra_derivatives,
Equations& eqns1) const
{
MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
// diregard residual
// Frac - derivatives are zero (they belong to upwind^2)
eqns.D()[seg][seg_upwind][SPres][SPres] += extra_derivatives.derivative(SPres + Indices::numEq);
eqns.D()[seg][seg_upwind][SPres][WQTotal] += extra_derivatives.derivative(WQTotal + Indices::numEq);
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
assemblePressureEq(const int seg,
@@ -228,7 +258,7 @@ assemblePressureEq(const int seg,
bool gfrac) const
{
MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
eqns.residual()[seg][SPres] = pressure_equation.value();
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) {

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

@@ -87,6 +87,17 @@ public:
const EvalWell& accelerationPressureLoss,
Equations& eqns) const;
void assembleHydroPressureLoss(const int seg,
const int seg_density,
const EvalWell& hydro_pressure_drop_seg,
Equations& eqns1) const;
void assemblePressureEqExtraDerivatives(const int seg,
const int seg_upwind,
const EvalWell& extra_derivatives,
Equations& eqns1) const;
//! \brief Assemble pressure terms.
void assemblePressureEq(const int seg,
const int seg_upwind,

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

@@ -214,25 +214,33 @@ template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
assembleDefaultPressureEq(const int seg,
WellState& well_state)
WellState& well_state,
const bool use_average_density)
{
assert(seg != 0); // not top segment
const int seg_upwind = segments_.upwinding_segment(seg);
const bool reverseFlow = seg != seg_upwind; // special treatment for reverse flow
// for top segment, the well control equation will be used.
EvalWell pressure_equation = primary_variables_.getSegmentPressure(seg);
EvalWell extra_derivatives;
// we need to handle the pressure difference between the two segments
// we only consider the hydrostatic pressure loss first
// hydrostatic pressure loss is assembled seperately at the end
// TODO: we might be able to add member variables to store these values, then we update well state
// after converged
const auto hydro_pressure_drop = segments_.getHydroPressureLoss(seg);
auto& ws = well_state.well(baseif_.indexOfWell());
auto& segments = ws.segments;
segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop.value();
pressure_equation -= hydro_pressure_drop;
if (this->frictionalPressureLossConsidered()) {
const auto friction_pressure_drop = segments_.getFrictionPressureLoss(seg);
const auto friction_pressure_drop = segments_.getFrictionPressureLoss(seg, false);
if (reverseFlow){
// call function once again to obtain/assemble remaining derivatives
extra_derivatives = -segments_.getFrictionPressureLoss(seg, true);
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assemblePressureEqExtraDerivatives(seg, seg_upwind, extra_derivatives, linSys_);
}
pressure_equation -= friction_pressure_drop;
segments.pressure_drop_friction[seg] = friction_pressure_drop.value();
}
@@ -241,7 +249,6 @@ assembleDefaultPressureEq(const int seg,
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const EvalWell outlet_pressure = primary_variables_.getSegmentPressure(outlet_segment_index);
const int seg_upwind = segments_.upwinding_segment(seg);
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assemblePressureEq(seg, seg_upwind, outlet_segment_index,
pressure_equation, outlet_pressure, linSys_);
@@ -249,6 +256,23 @@ assembleDefaultPressureEq(const int seg,
if (this->accelerationalPressureLossConsidered()) {
handleAccelerationPressureLoss(seg, well_state);
}
// Since density derivatives are organized differently than what is required for assemblePressureEq,
// this part needs to be assembled separately. Optionally use average density variant.
const auto hydro_pressure_drop_seg = segments_.getHydroPressureLoss(seg, seg);
if (!use_average_density){
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assembleHydroPressureLoss(seg, seg, hydro_pressure_drop_seg, linSys_);
segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop_seg.value();
} else {
const int seg_outlet = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const auto hydro_pressure_drop_outlet = segments_.getHydroPressureLoss(seg, seg_outlet);
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assembleHydroPressureLoss(seg, seg, 0.5*hydro_pressure_drop_seg, linSys_);
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assembleHydroPressureLoss(seg, seg_outlet, 0.5*hydro_pressure_drop_outlet, linSys_);
segments.pressure_drop_hydrostatic[seg] = 0.5*hydro_pressure_drop_seg.value() + 0.5*hydro_pressure_drop_outlet.value();
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
@@ -322,6 +346,7 @@ MultisegmentWellEval<FluidSystem,Indices,Scalar>::
assemblePressureEq(const int seg,
const UnitSystem& unit_system,
WellState& well_state,
const bool use_average_density,
DeferredLogger& deferred_logger)
{
switch(this->segmentSet()[seg].segmentType()) {
@@ -332,7 +357,7 @@ assemblePressureEq(const int seg,
break;
}
default :
assembleDefaultPressureEq(seg, well_state);
assembleDefaultPressureEq(seg, well_state, use_average_density);
}
}

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

@@ -76,7 +76,8 @@ protected:
void initMatrixAndVectors(const int num_cells);
void assembleDefaultPressureEq(const int seg,
WellState& well_state);
WellState& well_state,
const bool use_average_density);
// assemble pressure equation for ICD segments
void assembleICDPressureEq(const int seg,
@@ -88,6 +89,7 @@ protected:
void assemblePressureEq(const int seg,
const UnitSystem& unit_system,
WellState& well_state,
const bool use_average_density,
DeferredLogger& deferred_logger);
/// check whether the well equations get converged for this well

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

@@ -321,11 +321,13 @@ updateUpwindingSegments(const PrimaryVariables& primary_variables)
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
getHydroPressureLoss(const int seg) const
{
return densities_[seg] * well_.gravity() * depth_diffs_[seg];
getHydroPressureLoss(const int seg,
const int seg_density) const
{
return densities_[seg_density] * well_.gravity() * depth_diffs_[seg];
}
template<class FluidSystem, class Indices, class Scalar>
Scalar MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
getPressureDiffSegPerf(const int seg,
@@ -484,19 +486,41 @@ getSurfaceVolume(const EvalWell& temperature,
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
getFrictionPressureLoss(const int seg) const
getFrictionPressureLoss(const int seg, const bool return_extra_derivatives) const
{
const EvalWell mass_rate = mass_rates_[seg];
EvalWell mass_rate = mass_rates_[seg];
const int seg_upwind = upwinding_segments_[seg];
EvalWell density = densities_[seg_upwind];
EvalWell visc = viscosities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
// In the reverse flow case, we don't have enough slots for all derivatives, e.g.,
// upwind pressure and flow. We amend this by a second function call optioin, where
// only these remaining derivatives are considered.
// 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) {
density.clearDerivatives();
visc.clearDerivatives();
if (!return_extra_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);
visc.setDerivative(WQTotal, 0.0);
visc.setDerivative(SPres, 0.0);
} else {
if (PrimaryVariables::has_water){
constexpr int WFrac = Indices::numEq + PrimaryVariables::WFrac;
density.setDerivative(WFrac, 0.0);
visc.setDerivative(WFrac, 0.0);
}
if (PrimaryVariables::has_gas){
constexpr int GFrac = Indices::numEq + PrimaryVariables::GFrac;
density.setDerivative(GFrac, 0.0);
visc.setDerivative(GFrac, 0.0);
}
mass_rate.clearDerivatives();
}
}
const auto& segment_set = well_.wellEcl().getSegments();
const int outlet_segment_index = segment_set.segmentNumberToIndex(segment_set[seg].outletSegment());
const double length = segment_set[seg].totalLength() - segment_set[outlet_segment_index].totalLength();

5
opm/simulators/wells/MultisegmentWellSegments.hpp
View File

@@ -51,7 +51,8 @@ public:
//! \brief Update upwinding segments.
void updateUpwindingSegments(const PrimaryVariables& primary_variables);
EvalWell getHydroPressureLoss(const int seg) const;
EvalWell getHydroPressureLoss(const int seg,
const int seg_side) const;
//! Pressure difference between segment and perforation.
Scalar getPressureDiffSegPerf(const int seg,
@@ -63,7 +64,7 @@ public:
const int pvt_region_index,
const int seg_idx) const;
EvalWell getFrictionPressureLoss(const int seg) const;
EvalWell getFrictionPressureLoss(const int seg, const bool return_upwind_derivatives) const;
// pressure drop for Spiral ICD segment (WSEGSICD)
EvalWell pressureDropSpiralICD(const int seg) const;

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

@@ -1589,7 +1589,7 @@ namespace Opm
deferred_logger);
} else {
const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
this->assemblePressureEq(seg, unit_system, well_state, deferred_logger);
this->assemblePressureEq(seg, unit_system, well_state, this->param_.use_average_density_ms_wells_, deferred_logger);
}
}