MultisegmentWell: move updatePrimaryVariablesNewton to MultisegmentWellPrimaryVariables

This commit is contained in:
Arne Morten Kvarving 2022-12-19 09:52:48 +01:00
parent 6404d69201
commit b112a793c5
5 changed files with 64 additions and 60 deletions

View File

@ -171,55 +171,6 @@ getWellConvergence(const WellState& well_state,
return report;
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
updatePrimaryVariablesNewton(const BVectorWell& dwells,
const double relaxation_factor,
const double dFLimit,
const double max_pressure_change)
{
const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_.value_;
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
if (has_wfrac_variable) {
const int sign = dwells[seg][WFrac] > 0. ? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]) * relaxation_factor, dFLimit);
primary_variables_.value_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
}
if (has_gfrac_variable) {
const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]) * relaxation_factor, dFLimit);
primary_variables_.value_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
}
// handling the overshooting or undershooting of the fractions
primary_variables_.processFractions(seg);
// update the segment pressure
{
const int sign = dwells[seg][SPres] > 0.? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change);
primary_variables_.value_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5);
}
// update the total rate // TODO: should we have a limitation of the total rate change?
{
primary_variables_.value_[seg][WQTotal] = old_primary_variables[seg][WQTotal] - relaxation_factor * dwells[seg][WQTotal];
// make sure that no injector produce and no producer inject
if (seg == 0) {
if (baseif_.isInjector()) {
primary_variables_.value_[seg][WQTotal] = std::max( primary_variables_.value_[seg][WQTotal], 0.0);
} else {
primary_variables_.value_[seg][WQTotal] = std::min( primary_variables_.value_[seg][WQTotal], 0.0);
}
}
}
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::

View File

@ -133,12 +133,6 @@ protected:
void updateUpwindingSegments();
// updating the well_state based on well solution dwells
void updatePrimaryVariablesNewton(const BVectorWell& dwells,
const double relaxation_factor,
const double DFLimit,
const double max_pressure_change);
void computeSegmentFluidProperties(const EvalWell& temperature,
const EvalWell& saltConcentration,
int pvt_region_index,

View File

@ -142,6 +142,54 @@ update(const WellState& well_state)
}
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
updateNewton(const BVectorWell& dwells,
const double relaxation_factor,
const double dFLimit,
const double max_pressure_change)
{
const std::vector<std::array<double, numWellEq>> old_primary_variables = value_;
for (size_t seg = 0; seg < value_.size(); ++seg) {
if (has_wfrac_variable) {
const int sign = dwells[seg][WFrac] > 0. ? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]) * relaxation_factor, dFLimit);
value_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
}
if (has_gfrac_variable) {
const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]) * relaxation_factor, dFLimit);
value_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
}
// handling the overshooting or undershooting of the fractions
this->processFractions(seg);
// update the segment pressure
{
const int sign = dwells[seg][SPres] > 0.? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change);
value_[seg][SPres] = std::max(old_primary_variables[seg][SPres] - dx_limited, 1e5);
}
// update the total rate // TODO: should we have a limitation of the total rate change?
{
value_[seg][WQTotal] = old_primary_variables[seg][WQTotal] - relaxation_factor * dwells[seg][WQTotal];
// make sure that no injector produce and no producer inject
if (seg == 0) {
if (well_.isInjector()) {
value_[seg][WQTotal] = std::max(value_[seg][WQTotal], 0.0);
} else {
value_[seg][WQTotal] = std::min(value_[seg][WQTotal], 0.0);
}
}
}
}
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
processFractions(const int seg)

View File

@ -24,6 +24,8 @@
#include <opm/material/densead/Evaluation.hpp>
#include <opm/simulators/wells/MultisegmentWellEquations.hpp>
#include <array>
#include <vector>
@ -70,6 +72,9 @@ public:
using EvalWell = DenseAd::Evaluation<double, /*size=*/Indices::numEq + numWellEq>;
using Equations = MultisegmentWellEquations<Scalar,numWellEq,Indices::numEq>;
using BVectorWell = typename Equations::BVectorWell;
MultisegmentWellPrimaryVariables(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well)
: well_(well)
{}
@ -83,6 +88,12 @@ public:
//! \brief Copy values from well state.
void update(const WellState& well_state);
//! \brief Update values from newton update vector.
void updateNewton(const BVectorWell& dwells,
const double relaxation_factor,
const double DFLimit,
const double max_pressure_change);
// the values for the primary varibles
// based on different solutioin strategies, the wells can have different primary variables
std::vector<std::array<double, numWellEq> > value_;
@ -90,10 +101,10 @@ public:
// the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
std::vector<std::array<EvalWell, numWellEq> > evaluation_;
private:
//! \brief Handle non-reasonable fractions due to numerical overshoot.
void processFractions(const int seg);
private:
const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well_; //!< Reference to well interface
};

View File

@ -616,10 +616,10 @@ namespace Opm
const double dFLimit = this->param_.dwell_fraction_max_;
const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
this->MSWEval::updatePrimaryVariablesNewton(dwells,
relaxation_factor,
dFLimit,
max_pressure_change);
this->primary_variables_.updateNewton(dwells,
relaxation_factor,
dFLimit,
max_pressure_change);
this->updateWellStateFromPrimaryVariables(well_state, getRefDensity(), deferred_logger);
Base::calculateReservoirRates(well_state.well(this->index_of_well_));