OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-02 05:49:09 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5040 from steink/fix_acceleration_pressure_loss

Browse Source 
Include all derivatives for acceleration term in MS wells
This commit is contained in:
Kai Bao 2024-01-09 16:35:03 +01:00 committed by GitHub
commit e405f49def
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 103 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -199,22 +199,28 @@ assembleControlEq(const WellState& well_state,
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
assemblePressureLoss(const int seg,
const int seg_upwind,
const EvalWell& accelerationPressureLoss,
Equations& eqns1) const
{
assembleAccelerationTerm(const int seg_target,
const int seg,
const int seg_upwind,
const EvalWell& accelerationTerm,
Equations& eqns1) const
{ // seg_target: segment for which we are assembling the acc term
// seg: segment for wich we have computed the term
// seg_upwind: upwind segment to seg
// acceleration term shold be
// * velocity head for seg_target if seg = seg_target
// * negative velocity head for seg if seg != seg_target
MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
eqns.residual()[seg][SPres] -= accelerationPressureLoss.value();
eqns.D()[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq);
eqns.D()[seg][seg][SPres][WQTotal] -= accelerationPressureLoss.derivative(WQTotal + Indices::numEq);
eqns.residual()[seg_target][SPres] -= accelerationTerm.value();
eqns.D()[seg_target][seg][SPres][SPres] -= accelerationTerm.derivative(SPres + Indices::numEq);
eqns.D()[seg_target][seg][SPres][WQTotal] -= accelerationTerm.derivative(WQTotal + Indices::numEq);
if constexpr (has_wfrac_variable) {
eqns.D()[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq);
eqns.D()[seg_target][seg_upwind][SPres][WFrac] -= accelerationTerm.derivative(WFrac + Indices::numEq);
}
if constexpr (has_gfrac_variable) {
eqns.D()[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + Indices::numEq);
eqns.D()[seg_target][seg_upwind][SPres][GFrac] -= accelerationTerm.derivative(GFrac + Indices::numEq);
}
}
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
@ -239,7 +245,7 @@ assemblePressureEqExtraDerivatives(const int seg,
Equations& eqns1) const
{
MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
// diregard residual
// disregard 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);

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

@ -80,19 +80,20 @@ public:
Equations& eqns,
DeferredLogger& deferred_logger) const;
//! \brief Assemble piece of the acceleration term
void assembleAccelerationTerm(const int seg_target,
const int seg,
const int seg_upwing,
const EvalWell& accelerationTerm,
Equations& eqns1) const;
//! \brief Assemble pressure loss term.
void assemblePressureLoss(const int seg,
const int seg_upwind,
const EvalWell& accelerationPressureLoss,
Equations& eqns) const;
//! \brief Assemble hydraulic pressure term
void assembleHydroPressureLoss(const int seg,
const int seg_density,
const EvalWell& hydro_pressure_drop_seg,
Equations& eqns1) const;
//! \brief Assemble additional derivatives due to reverse flow
void assemblePressureEqExtraDerivatives(const int seg,
const int seg_upwind,
const EvalWell& extra_derivatives,

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

@ -198,20 +198,51 @@ extendEval(const Eval& in) const
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
handleAccelerationPressureLoss(const int seg,
WellState& well_state)
assembleAccelerationPressureLoss(const int seg,
WellState& well_state)
{
const EvalWell accelerationPressureLoss = segments_.accelerationPressureLoss(seg);
// Computes and assembles p-drop due to acceleration
assert(seg != 0); // top segment can not enter here
auto& segments = well_state.well(baseif_.indexOfWell()).segments;
segments.pressure_drop_accel[seg] = accelerationPressureLoss.value();
const auto& segment_set = this->segmentSet();
// Add segment head
const double seg_area = segment_set[seg].crossArea();
const EvalWell signed_velocity_head = segments_.accelerationPressureLossContribution(seg, seg_area);
segments.pressure_drop_accel[seg] = signed_velocity_head.value();
const int seg_upwind = segments_.upwinding_segment(seg);
// acceleration term is *subtracted* from pressure equation
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assemblePressureLoss(seg,
segments_.upwinding_segment(seg),
accelerationPressureLoss, linSys_);
assembleAccelerationTerm(seg, seg, seg_upwind, signed_velocity_head, linSys_);
if (seg != seg_upwind) {// special treatment for reverse flow
// extra derivatives are *added* to Jacobian (hence minus)
const EvalWell extra_derivatives = -segments_.accelerationPressureLossContribution(seg, seg_area, /*extra_derivatives*/ true);
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assemblePressureEqExtraDerivatives(seg, seg_upwind, extra_derivatives, linSys_);
}
// Subtract inlet head(s), opposite signs from above
for (const int inlet : segments_.inlets(seg)) {
// area used in formula is max of areas
const double inlet_area = std::max(seg_area, segment_set[inlet].crossArea());
const EvalWell signed_velocity_head_inlet = segments_.accelerationPressureLossContribution(inlet, inlet_area);
segments.pressure_drop_accel[seg] -= signed_velocity_head_inlet.value();
const int inlet_upwind = segments_.upwinding_segment(inlet);
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assembleAccelerationTerm(seg, inlet, inlet_upwind, -signed_velocity_head_inlet, linSys_);
if (inlet != inlet_upwind) {// special treatment for reverse flow
// extra derivatives are *added* to Jacobian (hence minus minus)
const EvalWell extra_derivatives_inlet = segments_.accelerationPressureLossContribution(inlet, inlet_area, /*extra_derivatives*/ true);
// in this case inlet_upwind = seg
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assemblePressureEqExtraDerivatives(seg, inlet_upwind, extra_derivatives_inlet, linSys_);
}
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
@ -352,7 +383,7 @@ assembleAccelerationAndHydroPressureLosses(const int seg,
const bool use_average_density)
{
if (this->accelerationalPressureLossConsidered()) {
handleAccelerationPressureLoss(seg, well_state);
assembleAccelerationPressureLoss(seg, well_state);
}
// Since density derivatives are organized differently than what is required for assemblePressureEq,

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

@ -128,8 +128,8 @@ protected:
const double tolerance_pressure_ms_wells,
DeferredLogger& deferred_logger) const;
void handleAccelerationPressureLoss(const int seg,
WellState& well_state);
void assembleAccelerationPressureLoss(const int seg,
WellState& well_state);
EvalWell pressureDropAutoICD(const int seg,
const UnitSystem& unit_system) const;

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

@ -829,43 +829,36 @@ pressureDropValve(const int seg,
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
accelerationPressureLoss(const int seg) const
accelerationPressureLossContribution(const int seg,
const double area,
const bool extra_reverse_flow_derivatives /*false*/) const
{
const auto& segment_set = well_.wellEcl().getSegments();
const double area = segment_set[seg].crossArea();
const EvalWell mass_rate = mass_rates_[seg];
// Compute the *signed* velocity head for given segment (sign is positive for flow towards surface, i.e., negative rate)
// Optionally return derivatives for reversed flow case
EvalWell mass_rate = mass_rates_[seg];
const int seg_upwind = upwinding_segments_[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 density = densities_[seg_upwind];
if (seg != seg_upwind) {
density.clearDerivatives();
}
EvalWell accelerationPressureLoss = mswellhelpers::velocityHead(area, mass_rate, density);
// handling the velocity head of intlet segments
for (const int inlet : inlets_[seg]) {
const int seg_upwind_inlet = upwinding_segments_[inlet];
const double inlet_area = segment_set[inlet].crossArea();
EvalWell inlet_density = densities_[seg_upwind_inlet];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (inlet != seg_upwind_inlet) {
inlet_density.clearDerivatives();
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){
constexpr int WFrac = Indices::numEq + PrimaryVariables::WFrac;
density.setDerivative(WFrac, 0.0);
}
if (PrimaryVariables::has_gas){
constexpr int GFrac = Indices::numEq + PrimaryVariables::GFrac;
density.setDerivative(GFrac, 0.0);
}
mass_rate.clearDerivatives();
}
const EvalWell inlet_mass_rate = mass_rates_[inlet];
accelerationPressureLoss -= mswellhelpers::velocityHead(std::max(inlet_area, area), inlet_mass_rate, inlet_density);
}
// We change the sign of the accelerationPressureLoss for injectors.
// Is this correct? Testing indicates that this is what the reference simulator does
const double sign = mass_rate < 0. ? 1.0 : - 1.0;
accelerationPressureLoss *= sign;
return accelerationPressureLoss;
}
const double sign = mass_rate > 0 ? -1.0 : 1.0;
return sign*mswellhelpers::velocityHead(area, mass_rate, density);
}
template <class FluidSystem, class Indices, class Scalar>
void

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

@ -88,14 +88,22 @@ public:
EvalWell pressureDropValve(const int seg,
const SummaryState& st,
const bool extra_reverse_flow_derivatives = false) const;
// pressure loss due to acceleration
EvalWell accelerationPressureLoss(const int seg) const;
// pressure loss contribution due to acceleration
EvalWell accelerationPressureLossContribution(const int seg,
const double area,
const bool extra_reverse_flow_derivatives = false) const;
const std::vector<std::vector<int>>& inlets() const
{
return inlets_;
}
const std::vector<int>& inlets(const int seg) const
{
return inlets_[seg];
}
const std::vector<std::vector<int>>& perforations() const
{
return perforations_;