diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index db4064e9e..1fbe7a60e 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -283,19 +283,31 @@ assembleICDPressureEq(const int seg, // the pressure equation is something like // p_seg - deltaP - p_outlet = 0. // the major part is how to calculate the deltaP + const int seg_upwind = segments_.upwinding_segment(seg); + const bool reverseFlow = seg != seg_upwind; // special treatment for reverse flow EvalWell pressure_equation = primary_variables_.getSegmentPressure(seg); EvalWell icd_pressure_drop; + EvalWell extra_derivatives; switch(this->segmentSet()[seg].segmentType()) { case Segment::SegmentType::SICD : - icd_pressure_drop = segments_.pressureDropSpiralICD(seg); + icd_pressure_drop = segments_.pressureDropSpiralICD(seg, /*extra derivatives*/false); + if (reverseFlow){ + extra_derivatives = -segments_.pressureDropSpiralICD(seg, /*extra derivatives*/true); + } break; case Segment::SegmentType::AICD : - icd_pressure_drop = segments_.pressureDropAutoICD(seg, unit_system); + icd_pressure_drop = segments_.pressureDropAutoICD(seg, unit_system, /*extra derivatives*/false); + if (reverseFlow){ + extra_derivatives = -segments_.pressureDropAutoICD(seg, unit_system, /*extra derivatives*/true); + } break; case Segment::SegmentType::VALVE : - icd_pressure_drop = segments_.pressureDropValve(seg); + icd_pressure_drop = segments_.pressureDropValve(seg, /*extra derivatives*/false); + if (reverseFlow){ + extra_derivatives = -segments_.pressureDropValve(seg, /*extra derivatives*/true); + } break; default: { OPM_DEFLOG_THROW(std::runtime_error, @@ -305,6 +317,12 @@ assembleICDPressureEq(const int seg, deferred_logger); } } + if (reverseFlow){ + MultisegmentWellAssemble(baseif_). + assemblePressureEqExtraDerivatives(seg, seg_upwind, extra_derivatives, linSys_); + + } + pressure_equation = pressure_equation - icd_pressure_drop; auto& ws = well_state.well(baseif_.indexOfWell()); ws.segments.pressure_drop_friction[seg] = icd_pressure_drop.value(); @@ -313,7 +331,6 @@ assembleICDPressureEq(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(baseif_). assemblePressureEq(seg, seg_upwind, outlet_segment_index, pressure_equation, outlet_pressure, diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index a307bb6a9..720e05e47 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -549,7 +549,8 @@ getFrictionPressureLoss(const int seg, const bool return_extra_derivatives) cons template typename MultisegmentWellSegments::EvalWell MultisegmentWellSegments:: -pressureDropSpiralICD(const int seg) const +pressureDropSpiralICD(const int seg, + const bool return_extra_derivatives) const { const auto& segment_set = well_.wellEcl().getSegments(); const SICD& sicd = segment_set[seg].spiralICD(); @@ -583,17 +584,39 @@ pressureDropSpiralICD(const int seg) const } EvalWell density = densities_[seg_upwind]; - // WARNING - // We disregard the derivatives from the upwind density to make sure derivatives - // wrt. to different segments dont get mixed. + EvalWell mass_rate = mass_rates_[seg]; + // 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 option, 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) { - water_fraction.clearDerivatives(); - water_viscosity.clearDerivatives(); - oil_fraction.clearDerivatives(); - oil_viscosity.clearDerivatives(); - gas_fraction.clearDerivatives(); - gas_viscosity.clearDerivatives(); - density.clearDerivatives(); + constexpr int nvar = FluidSystem::numPhases + 1; + std::vector zero_mask(nvar, false); + if (!return_extra_derivatives){ + zero_mask[PrimaryVariables::WQTotal] = true; + zero_mask[PrimaryVariables::SPres] = true; + } else { + if (PrimaryVariables::has_water){ + zero_mask[PrimaryVariables::WFrac] = true; + } + if (PrimaryVariables::has_gas){ + zero_mask[PrimaryVariables::GFrac] = true; + } + // mass_rate has no extra derivatives (they are organized as in equations) + mass_rate.clearDerivatives(); + } + for (int ii = 0; ii < nvar; ++ii) { + if (zero_mask[ii]) { + water_fraction.setDerivative(Indices::numEq + ii, 0.0); + water_viscosity.setDerivative(Indices::numEq + ii, 0.0); + oil_fraction.setDerivative(Indices::numEq + ii, 0.0); + oil_viscosity.setDerivative(Indices::numEq + ii, 0.0); + gas_fraction.setDerivative(Indices::numEq + ii, 0.0); + gas_viscosity.setDerivative(Indices::numEq + ii, 0.0); + density.setDerivative(Indices::numEq + ii, 0.0); + } + } } const EvalWell liquid_fraction = water_fraction + oil_fraction; @@ -604,7 +627,7 @@ pressureDropSpiralICD(const int seg) const const EvalWell mixture_viscosity = liquid_viscosity_fraction + gas_fraction * gas_viscosity; - const EvalWell reservoir_rate = mass_rates_[seg] / density; + const EvalWell reservoir_rate = mass_rate / density; const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor(); @@ -630,7 +653,8 @@ template typename MultisegmentWellSegments::EvalWell MultisegmentWellSegments:: pressureDropAutoICD(const int seg, - const UnitSystem& unit_system) const + const UnitSystem& unit_system, + const bool return_extra_derivatives) const { const auto& segment_set = well_.wellEcl().getSegments(); const AutoICD& aicd = segment_set[seg].autoICD(); @@ -671,20 +695,42 @@ pressureDropAutoICD(const int seg, } EvalWell density = densities_[seg_upwind]; - // WARNING - // We disregard the derivatives from the upwind density to make sure derivatives - // wrt. to different segments dont get mixed. + EvalWell mass_rate = mass_rates_[seg]; + // 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 option, 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) { - water_fraction.clearDerivatives(); - water_viscosity.clearDerivatives(); - water_density.clearDerivatives(); - oil_fraction.clearDerivatives(); - oil_viscosity.clearDerivatives(); - oil_density.clearDerivatives(); - gas_fraction.clearDerivatives(); - gas_viscosity.clearDerivatives(); - gas_density.clearDerivatives(); - density.clearDerivatives(); + constexpr int nvar = FluidSystem::numPhases + 1; + std::vector zero_mask(nvar, false); + if (!return_extra_derivatives){ + zero_mask[PrimaryVariables::WQTotal] = true; + zero_mask[PrimaryVariables::SPres] = true; + } else { + if (PrimaryVariables::has_water){ + zero_mask[PrimaryVariables::WFrac] = true; + } + if (PrimaryVariables::has_gas){ + zero_mask[PrimaryVariables::GFrac] = true; + } + // mass_rate has no extra derivatives (they are organized as in equations) + mass_rate.clearDerivatives(); + } + for (int ii = 0; ii < nvar; ++ii) { + if (zero_mask[ii]) { + water_fraction.setDerivative(Indices::numEq + ii, 0.0); + water_viscosity.setDerivative(Indices::numEq + ii, 0.0); + water_density.setDerivative(Indices::numEq + ii, 0.0); + oil_fraction.setDerivative(Indices::numEq + ii, 0.0); + oil_viscosity.setDerivative(Indices::numEq + ii, 0.0); + oil_density.setDerivative(Indices::numEq + ii, 0.0); + gas_fraction.setDerivative(Indices::numEq + ii, 0.0); + gas_viscosity.setDerivative(Indices::numEq + ii, 0.0); + gas_density.setDerivative(Indices::numEq + ii, 0.0); + density.setDerivative(Indices::numEq + ii, 0.0); + } + } } using MathTool = MathToolbox; @@ -698,7 +744,7 @@ pressureDropAutoICD(const int seg, const double rho_reference = aicd.densityCalibration(); const double visc_reference = aicd.viscosityCalibration(); - const auto volume_rate_icd = mass_rates_[seg] * aicd.scalingFactor() / mixture_density; + const auto volume_rate_icd = mass_rate * aicd.scalingFactor() / mixture_density; const double sign = volume_rate_icd <= 0. ? 1.0 : -1.0; // convert 1 unit volume rate using M = UnitSystem::measure; @@ -715,21 +761,42 @@ pressureDropAutoICD(const int seg, template typename MultisegmentWellSegments::EvalWell MultisegmentWellSegments:: -pressureDropValve(const int seg) const +pressureDropValve(const int seg, + const bool return_extra_derivatives) const { const auto& segment_set = well_.wellEcl().getSegments(); const Valve& valve = segment_set[seg].valve(); - const EvalWell& mass_rate = mass_rates_[seg]; + EvalWell mass_rate = mass_rates_[seg]; const int seg_upwind = upwinding_segments_[seg]; EvalWell visc = viscosities_[seg_upwind]; EvalWell density = densities_[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) { - visc.clearDerivatives(); - density.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 double additional_length = valve.pipeAdditionalLength(); diff --git a/opm/simulators/wells/MultisegmentWellSegments.hpp b/opm/simulators/wells/MultisegmentWellSegments.hpp index 282fddb9b..731835d7e 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.hpp +++ b/opm/simulators/wells/MultisegmentWellSegments.hpp @@ -74,15 +74,18 @@ public: EvalWell getFrictionPressureLoss(const int seg, const bool return_upwind_derivatives) const; // pressure drop for Spiral ICD segment (WSEGSICD) - EvalWell pressureDropSpiralICD(const int seg) const; + EvalWell pressureDropSpiralICD(const int seg, + const bool return_extra_derivatives) const; // pressure drop for Autonomous ICD segment (WSEGAICD) EvalWell pressureDropAutoICD(const int seg, - const UnitSystem& unit_system) const; + const UnitSystem& unit_system, + const bool return_extra_derivatives) const; // pressure drop for sub-critical valve (WSEGVALV) - EvalWell pressureDropValve(const int seg) const; - + EvalWell pressureDropValve(const int seg, + const bool return_extra_derivatives) const; + // pressure loss due to acceleration EvalWell accelerationPressureLoss(const int seg) const;