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 #4333 from akva2/msw_primaries

Browse Source 
added: MultisegmentWellPrimaryVariables
This commit is contained in:
Bård Skaflestad
2022-12-19 14:01:29 +01:00
committed by GitHub
parent c50cdc2454 bc0610171d
commit fe96a7181d
9 changed files with 892 additions and 725 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@@ -95,6 +95,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/MultisegmentWellEquations.cpp
opm/simulators/wells/MultisegmentWellEval.cpp
opm/simulators/wells/MultisegmentWellGeneric.cpp
opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
opm/simulators/wells/ParallelWellInfo.cpp
opm/simulators/wells/PerfData.cpp
opm/simulators/wells/SegmentState.cpp
@@ -380,6 +381,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/MultisegmentWellEquations.hpp
opm/simulators/wells/MultisegmentWellEval.hpp
opm/simulators/wells/MultisegmentWellGeneric.hpp
opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp
opm/simulators/wells/ParallelWellInfo.hpp
opm/simulators/wells/PerfData.hpp
opm/simulators/wells/PerforationData.hpp

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

@@ -68,11 +68,7 @@ namespace Opm
using typename MSWEval::Equations;
using typename MSWEval::EvalWell;
using typename MSWEval::BVectorWell;
using MSWEval::GFrac;
using MSWEval::WFrac;
using MSWEval::WQTotal;
using MSWEval::SPres;
using MSWEval::numWellEq;
using typename Base::PressureMatrix;
MultisegmentWell(const Well& well,

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

@@ -31,6 +31,7 @@
#include <opm/models/blackoil/blackoiltwophaseindices.hh>
#include <opm/simulators/wells/MultisegmentWellEquations.hpp>
#include <opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp>
#include <opm/simulators/wells/WellAssemble.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
@@ -87,9 +88,7 @@ assembleControlEq(const WellState& well_state,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
const double rho,
const EvalWell& wqTotal,
const EvalWell& bhp,
const std::function<EvalWell(const int)>& getQs,
const PrimaryVariables& primary_variables,
Equations& eqns1,
DeferredLogger& deferred_logger) const
{
@@ -104,19 +103,19 @@ assembleControlEq(const WellState& well_state,
auto getRates = [&]() {
std::vector<EvalWell> rates(3, 0.0);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[Water] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx));
rates[Water] = primary_variables.getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx));
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[Oil] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx));
rates[Oil] = primary_variables.getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx));
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[Gas] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx));
rates[Gas] = primary_variables.getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx));
}
return rates;
};
if (well_.wellIsStopped()) {
control_eq = wqTotal;
control_eq = primary_variables.getWQTotal();
} else if (well_.isInjector() ) {
// Find scaling factor to get injection rate,
const InjectorType injectorType = inj_controls.injector_type;
@@ -141,7 +140,7 @@ assembleControlEq(const WellState& well_state,
default:
throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
}
const EvalWell injection_rate = wqTotal / scaling;
const EvalWell injection_rate = primary_variables.getWQTotal() / scaling;
// Setup function for evaluation of BHP from THP (used only if needed).
std::function<EvalWell()> bhp_from_thp = [&]() {
const auto rates = getRates();
@@ -158,7 +157,7 @@ assembleControlEq(const WellState& well_state,
schedule,
summaryState,
inj_controls,
bhp,
primary_variables.getBhp(),
injection_rate,
bhp_from_thp,
control_eq,
@@ -181,7 +180,7 @@ assembleControlEq(const WellState& well_state,
schedule,
summaryState,
prod_controls,
bhp,
primary_variables.getBhp(),
rates,
bhp_from_thp,
control_eq,

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

@@ -26,14 +26,13 @@
#include <opm/input/eclipse/Schedule/Well/Well.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <functional>
namespace Opm
{
class DeferredLogger;
class GroupState;
template<class Scalar, int numWellEq, int numEq> class MultisegmentWellEquations;
template<class FluidSystem, class Indices, class Scalar> class MultisegmentWellPrimaryVariables;
class Schedule;
class SummaryState;
template<class FluidSystem, class Indices, class Scalar> class WellInterfaceIndices;
@@ -61,7 +60,9 @@ class MultisegmentWellAssemble
public:
static constexpr int numWellEq = Indices::numPhases+1;
using Equations = MultisegmentWellEquations<Scalar,numWellEq,Indices::numEq>;
using PrimaryVariables = MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>;
using EvalWell = DenseAd::Evaluation<Scalar, numWellEq+Indices::numEq>;
//! \brief Constructor initializes reference to well.
MultisegmentWellAssemble(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well)
: well_(well)
@@ -75,9 +76,7 @@ public:
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
const double rho,
const EvalWell& wqTotal,
const EvalWell& bhp,
const std::function<EvalWell(const int)>& getQs,
const PrimaryVariables& primary_variables,
Equations& eqns,
DeferredLogger& deferred_logger) const;

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

@@ -34,9 +34,7 @@
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
#include <opm/simulators/wells/RateConverter.hpp>
#include <opm/simulators/wells/WellAssemble.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellConvergence.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <opm/simulators/wells/WellState.hpp>
@@ -58,6 +56,7 @@ MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif)
: MultisegmentWellGeneric<Scalar>(baseif)
, baseif_(baseif)
, linSys_(*this)
, primary_variables_(baseif)
, upwinding_segments_(this->numberOfSegments(), 0)
, segment_densities_(this->numberOfSegments(), 0.0)
, segment_mass_rates_(this->numberOfSegments(), 0.0)
@@ -77,21 +76,6 @@ initMatrixAndVectors(const int num_cells)
{
linSys_.init(num_cells, baseif_.numPerfs(), baseif_.cells());
primary_variables_.resize(this->numberOfSegments());
primary_variables_evaluation_.resize(this->numberOfSegments());
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
initPrimaryVariablesEvaluation()
{
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
primary_variables_evaluation_[seg][eq_idx] = 0.0;
primary_variables_evaluation_[seg][eq_idx].setValue(primary_variables_[seg][eq_idx]);
primary_variables_evaluation_[seg][eq_idx].setDerivative(eq_idx + Indices::numEq, 1.0);
}
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
@@ -110,7 +94,8 @@ getWellConvergence(const WellState& well_state,
assert(int(B_avg.size()) == baseif_.numComponents());
// checking if any residual is NaN or too large. The two large one is only handled for the well flux
std::vector<std::vector<double>> abs_residual(this->numberOfSegments(), std::vector<double>(numWellEq, 0.0));
std::vector<std::vector<double>> abs_residual(this->numberOfSegments(),
std::vector<double>(numWellEq, 0.0));
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
abs_residual[seg][eq_idx] = std::abs(linSys_.residual()[seg][eq_idx]);
@@ -185,312 +170,6 @@ getWellConvergence(const WellState& well_state,
return report;
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
processFractions(const int seg)
{
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
const PhaseUsage& pu = baseif_.phaseUsage();
std::vector<double> fractions(baseif_.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::waterPhaseIdx) ) {
const int water_pos = pu.phase_pos[Water];
fractions[water_pos] = primary_variables_[seg][WFrac];
fractions[oil_pos] -= fractions[water_pos];
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
const int gas_pos = pu.phase_pos[Gas];
fractions[gas_pos] = primary_variables_[seg][GFrac];
fractions[oil_pos] -= fractions[gas_pos];
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water];
if (fractions[water_pos] < 0.0) {
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]);
}
fractions[oil_pos] /= (1.0 - fractions[water_pos]);
fractions[water_pos] = 0.0;
}
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
if (fractions[gas_pos] < 0.0) {
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
fractions[pu.phase_pos[Water]] /= (1.0 - fractions[gas_pos]);
}
fractions[oil_pos] /= (1.0 - fractions[gas_pos]);
fractions[gas_pos] = 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::gasPhaseIdx) ) {
fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]);
}
fractions[oil_pos] = 0.0;
}
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]];
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]];
}
}
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_;
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_[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_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
}
// handling the overshooting or undershooting of the fractions
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_[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_[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_[seg][WQTotal] = std::max( primary_variables_[seg][WQTotal], 0.0);
} else {
primary_variables_[seg][WQTotal] = std::min( primary_variables_[seg][WQTotal], 0.0);
}
}
}
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
updatePrimaryVariables(const WellState& well_state)
{
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Gas = BlackoilPhases::Vapour;
// TODO: to test using rate conversion coefficients to see if it will be better than
// this default one
if (!baseif_.isOperableAndSolvable() && !baseif_.wellIsStopped()) return;
const Well& well = baseif_.wellEcl();
// the index of the top segment in the WellState
const auto& ws = well_state.well(baseif_.indexOfWell());
const auto& segments = ws.segments;
// maybe a throw for parallel running?
assert(int(segments.size()) == this->numberOfSegments());
const auto& segment_rates = segments.rates;
const auto& segment_pressure = segments.pressure;
const PhaseUsage& pu = baseif_.phaseUsage();
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
// calculate the total rate for each segment
double total_seg_rate = 0.0;
// the segment pressure
primary_variables_[seg][SPres] = segment_pressure[seg];
// TODO: under what kind of circustances, the following will be wrong?
// the definition of g makes the gas phase is always the last phase
for (int p = 0; p < baseif_.numPhases(); p++) {
total_seg_rate += baseif_.scalingFactor(p) * segment_rates[baseif_.numPhases() * seg + p];
}
if (seg == 0) {
if (baseif_.isInjector()) {
total_seg_rate = std::max(total_seg_rate, 0.);
} else {
total_seg_rate = std::min(total_seg_rate, 0.);
}
}
primary_variables_[seg][WQTotal] = total_seg_rate;
if (std::abs(total_seg_rate) > 0.) {
if (has_wfrac_variable) {
const int water_pos = pu.phase_pos[Water];
primary_variables_[seg][WFrac] = baseif_.scalingFactor(water_pos) * segment_rates[baseif_.numPhases() * seg + water_pos] / total_seg_rate;
}
if (has_gfrac_variable) {
const int gas_pos = pu.phase_pos[Gas];
primary_variables_[seg][GFrac] = baseif_.scalingFactor(gas_pos) * segment_rates[baseif_.numPhases() * seg + gas_pos] / total_seg_rate;
}
} else { // total_seg_rate == 0
if (baseif_.isInjector()) {
// only single phase injection handled
auto phase = well.getInjectionProperties().injectorType;
if (has_wfrac_variable) {
if (phase == InjectorType::WATER) {
primary_variables_[seg][WFrac] = 1.0;
} else {
primary_variables_[seg][WFrac] = 0.0;
}
}
if (has_gfrac_variable) {
if (phase == InjectorType::GAS) {
primary_variables_[seg][GFrac] = 1.0;
} else {
primary_variables_[seg][GFrac] = 0.0;
}
}
} else if (baseif_.isProducer()) { // producers
if (has_wfrac_variable) {
primary_variables_[seg][WFrac] = 1.0 / baseif_.numPhases();
}
if (has_gfrac_variable) {
primary_variables_[seg][GFrac] = 1.0 / baseif_.numPhases();
}
}
}
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
volumeFraction(const int seg,
const unsigned compIdx) const
{
if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
return primary_variables_evaluation_[seg][WFrac];
}
if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
return primary_variables_evaluation_[seg][GFrac];
}
// Oil fraction
EvalWell oil_fraction = 1.0;
if (has_wfrac_variable) {
oil_fraction -= primary_variables_evaluation_[seg][WFrac];
}
if (has_gfrac_variable) {
oil_fraction -= primary_variables_evaluation_[seg][GFrac];
}
/* if (has_solvent) {
oil_fraction -= primary_variables_evaluation_[seg][SFrac];
} */
return oil_fraction;
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
volumeFractionScaled(const int seg,
const int comp_idx) const
{
// For reservoir rate control, the distr in well control is used for the
// rate conversion coefficients. For the injection well, only the distr of the injection
// phase is not zero.
const double scale = baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
if (scale > 0.) {
return volumeFraction(seg, comp_idx) / scale;
}
return volumeFraction(seg, comp_idx);
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
surfaceVolumeFraction(const int seg,
const int comp_idx) const
{
EvalWell sum_volume_fraction_scaled = 0.;
for (int idx = 0; idx < baseif_.numComponents(); ++idx) {
sum_volume_fraction_scaled += volumeFractionScaled(seg, idx);
}
assert(sum_volume_fraction_scaled.value() != 0.);
return volumeFractionScaled(seg, comp_idx) / sum_volume_fraction_scaled;
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getSegmentRateUpwinding(const int seg,
const size_t comp_idx) const
{
const int seg_upwind = upwinding_segments_[seg];
// the result will contain the derivative with respect to WQTotal in segment seg,
// and the derivatives with respect to WFrac GFrac in segment seg_upwind.
// the derivative with respect to SPres should be zero.
if (seg == 0 && baseif_.isInjector()) {
const Well& well = baseif_.wellEcl();
auto phase = well.getInjectionProperties().injectorType;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
&& Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == comp_idx
&& phase == InjectorType::WATER)
return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
&& Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == comp_idx
&& phase == InjectorType::OIL)
return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
&& Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == comp_idx
&& phase == InjectorType::GAS)
return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
return 0.0;
}
const EvalWell segment_rate = primary_variables_evaluation_[seg][WQTotal] * volumeFractionScaled(seg_upwind, comp_idx);
assert(segment_rate.derivative(SPres + Indices::numEq) == 0.);
return segment_rate;
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
@@ -527,14 +206,14 @@ computeSegmentFluidProperties(const EvalWell& temperature,
// the compostion of the components inside wellbore under surface condition
std::vector<EvalWell> mix_s(baseif_.numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) {
mix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx);
mix_s[comp_idx] = primary_variables_.surfaceVolumeFraction(seg, comp_idx);
}
std::vector<EvalWell> b(baseif_.numComponents(), 0.0);
std::vector<EvalWell> visc(baseif_.numComponents(), 0.0);
std::vector<EvalWell>& phase_densities = segment_phase_densities_[seg];
const EvalWell seg_pressure = getSegmentPressure(seg);
const EvalWell seg_pressure = primary_variables_.getSegmentPressure(seg);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
b[waterCompIdx] =
@@ -671,61 +350,15 @@ computeSegmentFluidProperties(const EvalWell& temperature,
// calculate the mass rates
segment_mass_rates_[seg] = 0.;
for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) {
const EvalWell rate = getSegmentRateUpwinding(seg, comp_idx);
const int upwind_seg = upwinding_segments_[seg];
const EvalWell rate = primary_variables_.getSegmentRateUpwinding(seg,
upwind_seg,
comp_idx);
this->segment_mass_rates_[seg] += rate * surf_dens[comp_idx];
}
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getSegmentPressure(const int seg) const
{
return primary_variables_evaluation_[seg][SPres];
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getBhp() const
{
return getSegmentPressure(0);
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getSegmentRate(const int seg,
const int comp_idx) const
{
return primary_variables_evaluation_[seg][WQTotal] * volumeFractionScaled(seg, comp_idx);
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getQs(const int comp_idx) const
{
return getSegmentRate(0, comp_idx);
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getSegmentWQTotal(const int seg) const
{
return primary_variables_evaluation_[seg][WQTotal];
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getWQTotal() const
{
return getSegmentWQTotal(0);
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
@@ -972,11 +605,11 @@ getSegmentSurfaceVolume(const EvalWell& temperature,
const int pvt_region_index,
const int seg_idx) const
{
const EvalWell seg_pressure = getSegmentPressure(seg_idx);
const EvalWell seg_pressure = primary_variables_.getSegmentPressure(seg_idx);
std::vector<EvalWell> mix_s(baseif_.numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) {
mix_s[comp_idx] = surfaceVolumeFraction(seg_idx, comp_idx);
mix_s[comp_idx] = primary_variables_.surfaceVolumeFraction(seg_idx, comp_idx);
}
std::vector<EvalWell> b(baseif_.numComponents(), 0.);
@@ -1162,7 +795,7 @@ assembleDefaultPressureEq(const int seg,
assert(seg != 0); // not top segment
// for top segment, the well control equation will be used.
EvalWell pressure_equation = getSegmentPressure(seg);
EvalWell pressure_equation = primary_variables_.getSegmentPressure(seg);
// we need to handle the pressure difference between the two segments
// we only consider the hydrostatic pressure loss first
@@ -1182,7 +815,7 @@ assembleDefaultPressureEq(const int seg,
// contribution from the outlet segment
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);
const EvalWell outlet_pressure = primary_variables_.getSegmentPressure(outlet_segment_index);
const int seg_upwind = upwinding_segments_[seg];
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
@@ -1194,201 +827,6 @@ assembleDefaultPressureEq(const int seg,
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
updateWellStateFromPrimaryVariables(WellState& well_state,
const double rho,
DeferredLogger& deferred_logger) const
{
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;
const auto pvtReg = std::max(this->baseif_.wellEcl().pvt_table_number() - 1, 0);
const PhaseUsage& pu = baseif_.phaseUsage();
assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
const int oil_pos = pu.phase_pos[Oil];
auto& ws = well_state.well(baseif_.indexOfWell());
auto& segments = ws.segments;
auto& segment_rates = segments.rates;
auto& disgas = segments.dissolved_gas_rate;
auto& vapoil = segments.vaporized_oil_rate;
auto& segment_pressure = segments.pressure;
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
std::vector<double> fractions(baseif_.numPhases(), 0.0);
fractions[oil_pos] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water];
fractions[water_pos] = primary_variables_[seg][WFrac];
fractions[oil_pos] -= fractions[water_pos];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
fractions[gas_pos] = primary_variables_[seg][GFrac];
fractions[oil_pos] -= fractions[gas_pos];
}
// convert the fractions to be Q_p / G_total to calculate the phase rates
for (int p = 0; p < baseif_.numPhases(); ++p) {
const double scale = baseif_.scalingFactor(p);
// for injection wells, there should only one non-zero scaling factor
if (scale > 0.) {
fractions[p] /= scale;
} else {
// this should only happens to injection wells
fractions[p] = 0.;
}
}
// calculate the phase rates based on the primary variables
const double g_total = primary_variables_[seg][WQTotal];
for (int p = 0; p < baseif_.numPhases(); ++p) {
const double phase_rate = g_total * fractions[p];
segment_rates[seg*baseif_.numPhases() + p] = phase_rate;
if (seg == 0) { // top segment
ws.surface_rates[p] = phase_rate;
}
}
// update the segment pressure
segment_pressure[seg] = primary_variables_[seg][SPres];
if (seg == 0) { // top segment
ws.bhp = segment_pressure[seg];
}
// Calculate other per-phase dynamic quantities.
const auto temperature = 0.0; // Ignore thermal effects
const auto saltConc = 0.0; // Ignore salt precipitation
const auto Rvw = 0.0; // Ignore vaporised water.
auto rsMax = 0.0;
auto rvMax = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// Both oil and gas active.
rsMax = FluidSystem::oilPvt()
.saturatedGasDissolutionFactor(pvtReg, temperature, segment_pressure[seg]);
rvMax = FluidSystem::gasPvt()
.saturatedOilVaporizationFactor(pvtReg, temperature, segment_pressure[seg]);
}
// 1) Infer phase splitting for oil/gas.
const auto& [Rs, Rv] = this->baseif_.rateConverter().inferDissolvedVaporisedRatio
(rsMax, rvMax, segment_rates.begin() + (seg + 0)*this->baseif_.numPhases());
if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
vapoil[seg] = disgas[seg] = 0.0;
}
else {
const auto* qs = &segment_rates[seg*this->baseif_.numPhases() + 0];
const auto denom = 1.0 - (Rs * Rv);
const auto io = pu.phase_pos[Oil];
const auto ig = pu.phase_pos[Gas];
disgas[seg] = Rs * (qs[io] - Rv*qs[ig]) / denom;
vapoil[seg] = Rv * (qs[ig] - Rs*qs[io]) / denom;
}
// 2) Local condition volume flow rates
{
// Use std::span<> in C++20 and beyond.
const auto rate_start = (seg + 0) * this->baseif_.numPhases();
const auto* surf_rates = segment_rates.data() + rate_start;
auto* resv_rates = segments.phase_resv_rates.data() + rate_start;
this->baseif_.rateConverter().calcReservoirVoidageRates
(pvtReg, segment_pressure[seg],
std::max(0.0, Rs),
std::max(0.0, Rv),
temperature, saltConc, surf_rates, resv_rates);
}
// 3) Local condition holdup fractions.
const auto tot_resv =
std::accumulate(segments.phase_resv_rates.begin() + (seg + 0)*this->baseif_.numPhases(),
segments.phase_resv_rates.begin() + (seg + 1)*this->baseif_.numPhases(),
0.0);
std::transform(segments.phase_resv_rates.begin() + (seg + 0)*this->baseif_.numPhases(),
segments.phase_resv_rates.begin() + (seg + 1)*this->baseif_.numPhases(),
segments.phase_holdup.begin() + (seg + 0)*this->baseif_.numPhases(),
[tot_resv](const auto qr) { return std::clamp(qr / tot_resv, 0.0, 1.0); });
// 4) Local condition flow velocities for segments other than top segment.
if (seg > 0) {
// Possibly poor approximation
// Velocity = Flow rate / cross-sectional area.
// Additionally ignores drift flux.
const auto area = this->baseif_.wellEcl().getSegments()
.getFromSegmentNumber(segments.segment_number()[seg]).crossArea();
const auto velocity = (area > 0.0) ? tot_resv / area : 0.0;
std::transform(segments.phase_holdup.begin() + (seg + 0)*this->baseif_.numPhases(),
segments.phase_holdup.begin() + (seg + 1)*this->baseif_.numPhases(),
segments.phase_velocity.begin() + (seg + 0)*this->baseif_.numPhases(),
[velocity](const auto hf) { return (hf > 0.0) ? velocity : 0.0; });
}
// 5) Local condition phase viscosities.
segments.phase_viscosity[seg*this->baseif_.numPhases() + pu.phase_pos[Oil]] =
FluidSystem::oilPvt().viscosity(pvtReg, temperature, segment_pressure[seg], Rs);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
segments.phase_viscosity[seg*this->baseif_.numPhases() + pu.phase_pos[Water]] =
FluidSystem::waterPvt().viscosity(pvtReg, temperature, segment_pressure[seg], saltConc);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
segments.phase_viscosity[seg*this->baseif_.numPhases() + pu.phase_pos[Gas]] =
FluidSystem::gasPvt().viscosity(pvtReg, temperature, segment_pressure[seg], Rv, Rvw);
}
}
// Segment flow velocity in top segment.
{
const auto np = this->baseif_.numPhases();
auto segVel = [&segments, np](const auto segmentNumber)
{
auto v = 0.0;
const auto* vel = segments.phase_velocity.data() + segmentNumber*np;
for (auto p = 0*np; p < np; ++p) {
if (std::abs(vel[p]) > std::abs(v)) {
v = vel[p];
}
}
return v;
};
const auto seg = 0;
auto maxVel = 0.0;
for (const auto& inlet : this->segmentSet()[seg].inletSegments()) {
const auto v = segVel(this->segmentNumberToIndex(inlet));
if (std::abs(v) > std::abs(maxVel)) {
maxVel = v;
}
}
std::transform(segments.phase_holdup.begin() + (seg + 0)*this->baseif_.numPhases(),
segments.phase_holdup.begin() + (seg + 1)*this->baseif_.numPhases(),
segments.phase_velocity.begin() + (seg + 0)*this->baseif_.numPhases(),
[maxVel](const auto hf) { return (hf > 0.0) ? maxVel : 0.0; });
}
WellBhpThpCalculator(this->baseif_)
.updateThp(rho, [this]() { return this->baseif_.wellEcl().alq_value(); },
{FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)},
well_state, deferred_logger);
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
@@ -1405,7 +843,7 @@ assembleICDPressureEq(const int seg,
(segment.segmentType() == Segment::SegmentType::VALVE) &&
(segment.valve().status() == Opm::ICDStatus::SHUT) ) { // we use a zero rate equation to handle SHUT valve
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assembleTrivialEq(seg, this->primary_variables_evaluation_[seg][WQTotal].value(), linSys_);
assembleTrivialEq(seg, this->primary_variables_.eval(seg)[WQTotal].value(), linSys_);
auto& ws = well_state.well(baseif_.indexOfWell());
ws.segments.pressure_drop_friction[seg] = 0.;
@@ -1416,7 +854,7 @@ assembleICDPressureEq(const int seg,
// p_seg - deltaP - p_outlet = 0.
// the major part is how to calculate the deltaP
EvalWell pressure_equation = getSegmentPressure(seg);
EvalWell pressure_equation = primary_variables_.getSegmentPressure(seg);
EvalWell icd_pressure_drop;
switch(this->segmentSet()[seg].segmentType()) {
@@ -1440,7 +878,7 @@ assembleICDPressureEq(const int seg,
// contribution from the outlet segment
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);
const EvalWell outlet_pressure = primary_variables_.getSegmentPressure(outlet_segment_index);
const int seg_upwind = upwinding_segments_[seg];
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
@@ -1620,31 +1058,6 @@ getResidualMeasureValue(const WellState& well_state,
return sum;
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
updateUpwindingSegments()
{
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
// special treatment is needed for segment 0
if (seg == 0) {
// we are not supposed to have injecting producers and producing injectors
assert( ! (baseif_.isProducer() && primary_variables_evaluation_[seg][WQTotal] > 0.) );
assert( ! (baseif_.isInjector() && primary_variables_evaluation_[seg][WQTotal] < 0.) );
upwinding_segments_[seg] = seg;
continue;
}
// for other normal segments
if (primary_variables_evaluation_[seg][WQTotal] <= 0.) {
upwinding_segments_[seg] = seg;
} else {
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
upwinding_segments_[seg] = outlet_segment_index;
}
}
}
#define INSTANCE(...) \
template class MultisegmentWellEval<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;

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

@@ -24,6 +24,7 @@
#include <opm/simulators/wells/MultisegmentWellEquations.hpp>
#include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
#include <opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp>
#include <opm/material/densead/Evaluation.hpp>
@@ -52,36 +53,10 @@ template<typename FluidSystem, typename Indices, typename Scalar>
class MultisegmentWellEval : public MultisegmentWellGeneric<Scalar>
{
protected:
// TODO: for now, not considering the polymer, solvent and so on to simplify the development process.
// TODO: we need to have order for the primary variables and also the order for the well equations.
// sometimes, they are similar, while sometimes, they can have very different forms.
// Table showing the primary variable indices, depending on what phases are present:
//
// 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
// 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 int WQTotal = 0;
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;
// the number of well equations TODO: it should have a more general strategy for it
static constexpr int numWellEq = Indices::numPhases + 1;
using PrimaryVariables = MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>;
static constexpr int numWellEq = PrimaryVariables::numWellEq;
static constexpr int SPres = PrimaryVariables::SPres;
static constexpr int WQTotal = PrimaryVariables::WQTotal;
using Equations = MultisegmentWellEquations<Scalar,numWellEq,Indices::numEq>;
@@ -91,7 +66,7 @@ protected:
// TODO: for more efficient implementation, we should have EvalReservoir, EvalWell, and EvalRerservoirAndWell
// EvalR (Eval), EvalW, EvalRW
// TODO: for now, we only use one type to save some implementation efforts, while improve later.
using EvalWell = DenseAd::Evaluation<double, /*size=*/Indices::numEq + numWellEq>;
using EvalWell = typename PrimaryVariables::EvalWell;
using Eval = DenseAd::Evaluation<Scalar, /*size=*/Indices::numEq>;
public:
@@ -103,7 +78,6 @@ protected:
MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif);
void initMatrixAndVectors(const int num_cells);
void initPrimaryVariablesEvaluation();
void assembleDefaultPressureEq(const int seg,
WellState& well_state);
@@ -131,39 +105,17 @@ protected:
const double relaxed_inner_tolerance_pressure_ms_well,
const bool relax_tolerance) const;
// handling the overshooting and undershooting of the fractions
void processFractions(const int seg);
void updatePrimaryVariables(const WellState& well_state);
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,
DeferredLogger& deferred_logger);
EvalWell getBhp() const;
EvalWell getFrictionPressureLoss(const int seg) const;
EvalWell getHydroPressureLoss(const int seg) const;
EvalWell getQs(const int comp_idx) const;
EvalWell getSegmentWQTotal(const int seg) const;
EvalWell getSegmentPressure(const int seg) const;
EvalWell getSegmentRate(const int seg, const int comp_idx) const;
EvalWell getSegmentRateUpwinding(const int seg,
const size_t comp_idx) const;
EvalWell getSegmentSurfaceVolume(const EvalWell& temperature,
const EvalWell& saltConcentration,
const int pvt_region_index,
const int seg_idx) const;
EvalWell getWQTotal() const;
std::pair<bool, std::vector<Scalar> >
getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
@@ -193,23 +145,6 @@ protected:
// pressure drop for sub-critical valve (WSEGVALV)
EvalWell pressureDropValve(const int seg) const;
void updateWellStateFromPrimaryVariables(WellState& well_state,
const double rho,
DeferredLogger& deferred_logger) const;
// fraction value of the primary variables
// should we just use member variables to store them instead of calculating them again and again
EvalWell volumeFraction(const int seg,
const unsigned compIdx) const;
// F_p / g_p, the basic usage of this value is because Q_p = G_t * F_p / G_p
EvalWell volumeFractionScaled(const int seg,
const int comp_idx) const;
// basically Q_p / \sigma_p Q_p
EvalWell surfaceVolumeFraction(const int seg,
const int comp_idx) const;
// convert a Eval from reservoir to contain the derivative related to wells
EvalWell extendEval(const Eval& in) const;
@@ -217,12 +152,7 @@ protected:
Equations linSys_; //!< The equation system
// the values for the primary varibles
// based on different solutioin strategies, the wells can have different primary variables
std::vector<std::array<double, numWellEq> > primary_variables_;
// the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
std::vector<std::array<EvalWell, numWellEq> > primary_variables_evaluation_;
PrimaryVariables primary_variables_; //!< The primary variables
// the upwinding segment for each segment based on the flow direction
std::vector<int> upwinding_segments_;

658
opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp Normal file
View File

@@ -0,0 +1,658 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2016 - 2017 IRIS AS.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp>
#include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/models/blackoil/blackoilindices.hh>
#include <opm/models/blackoil/blackoilonephaseindices.hh>
#include <opm/models/blackoil/blackoiltwophaseindices.hh>
#include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
#include <opm/simulators/wells/RateConverter.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <algorithm>
namespace Opm {
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
resize(const int numSegments)
{
value_.resize(numSegments);
evaluation_.resize(numSegments);
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
init()
{
for (size_t seg = 0; seg < value_.size(); ++seg) {
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
evaluation_[seg][eq_idx] = 0.0;
evaluation_[seg][eq_idx].setValue(value_[seg][eq_idx]);
evaluation_[seg][eq_idx].setDerivative(eq_idx + Indices::numEq, 1.0);
}
}
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
update(const WellState& well_state)
{
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Gas = BlackoilPhases::Vapour;
// TODO: to test using rate conversion coefficients to see if it will be better than
// this default one
if (!well_.isOperableAndSolvable() && !well_.wellIsStopped())
return;
const Well& well = well_.wellEcl();
// the index of the top segment in the WellState
const auto& ws = well_state.well(well_.indexOfWell());
const auto& segments = ws.segments;
// maybe a throw for parallel running?
assert(segments.size() == value_.size());
const auto& segment_rates = segments.rates;
const auto& segment_pressure = segments.pressure;
const PhaseUsage& pu = well_.phaseUsage();
for (size_t seg = 0; seg < value_.size(); ++seg) {
// calculate the total rate for each segment
double total_seg_rate = 0.0;
// the segment pressure
value_[seg][SPres] = segment_pressure[seg];
// TODO: under what kind of circustances, the following will be wrong?
// the definition of g makes the gas phase is always the last phase
for (int p = 0; p < well_.numPhases(); p++) {
total_seg_rate += well_.scalingFactor(p) * segment_rates[well_.numPhases() * seg + p];
}
if (seg == 0) {
if (well_.isInjector()) {
total_seg_rate = std::max(total_seg_rate, 0.);
} else {
total_seg_rate = std::min(total_seg_rate, 0.);
}
}
value_[seg][WQTotal] = total_seg_rate;
if (std::abs(total_seg_rate) > 0.) {
if (has_wfrac_variable) {
const int water_pos = pu.phase_pos[Water];
value_[seg][WFrac] = well_.scalingFactor(water_pos) * segment_rates[well_.numPhases() * seg + water_pos] / total_seg_rate;
}
if (has_gfrac_variable) {
const int gas_pos = pu.phase_pos[Gas];
value_[seg][GFrac] = well_.scalingFactor(gas_pos) * segment_rates[well_.numPhases() * seg + gas_pos] / total_seg_rate;
}
} else { // total_seg_rate == 0
if (well_.isInjector()) {
// only single phase injection handled
auto phase = well.getInjectionProperties().injectorType;
if (has_wfrac_variable) {
if (phase == InjectorType::WATER) {
value_[seg][WFrac] = 1.0;
} else {
value_[seg][WFrac] = 0.0;
}
}
if (has_gfrac_variable) {
if (phase == InjectorType::GAS) {
value_[seg][GFrac] = 1.0;
} else {
value_[seg][GFrac] = 0.0;
}
}
} else if (well_.isProducer()) { // producers
if (has_wfrac_variable) {
value_[seg][WFrac] = 1.0 / well_.numPhases();
}
if (has_gfrac_variable) {
value_[seg][GFrac] = 1.0 / well_.numPhases();
}
}
}
}
}
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>::
copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
const double rho,
WellState& well_state,
DeferredLogger& deferred_logger) const
{
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;
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());
auto& segments = ws.segments;
auto& segment_rates = segments.rates;
auto& disgas = segments.dissolved_gas_rate;
auto& vapoil = segments.vaporized_oil_rate;
auto& segment_pressure = segments.pressure;
for (size_t seg = 0; seg < value_.size(); ++seg) {
std::vector<double> 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::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
fractions[gas_pos] = value_[seg][GFrac];
fractions[oil_pos] -= fractions[gas_pos];
}
// convert the fractions to be Q_p / G_total to calculate the phase rates
for (int p = 0; p < well_.numPhases(); ++p) {
const double scale = well_.scalingFactor(p);
// for injection wells, there should only one non-zero scaling factor
if (scale > 0.) {
fractions[p] /= scale;
} else {
// this should only happens to injection wells
fractions[p] = 0.;
}
}
// calculate the phase rates based on the primary variables
const double g_total = value_[seg][WQTotal];
for (int p = 0; p < well_.numPhases(); ++p) {
const double phase_rate = g_total * fractions[p];
segment_rates[seg * well_.numPhases() + p] = phase_rate;
if (seg == 0) { // top segment
ws.surface_rates[p] = phase_rate;
}
}
// update the segment pressure
segment_pressure[seg] = value_[seg][SPres];
if (seg == 0) { // top segment
ws.bhp = segment_pressure[seg];
}
// Calculate other per-phase dynamic quantities.
const auto temperature = 0.0; // Ignore thermal effects
const auto saltConc = 0.0; // Ignore salt precipitation
const auto Rvw = 0.0; // Ignore vaporised water.
auto rsMax = 0.0;
auto rvMax = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// Both oil and gas active.
rsMax = FluidSystem::oilPvt()
.saturatedGasDissolutionFactor(pvtReg, temperature, segment_pressure[seg]);
rvMax = FluidSystem::gasPvt()
.saturatedOilVaporizationFactor(pvtReg, temperature, segment_pressure[seg]);
}
// 1) Infer phase splitting for oil/gas.
const auto& [Rs, Rv] = well_.rateConverter().inferDissolvedVaporisedRatio
(rsMax, rvMax, segment_rates.begin() + seg * well_.numPhases());
if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
vapoil[seg] = disgas[seg] = 0.0;
}
else {
const auto* qs = &segment_rates[seg * well_.numPhases()];
const auto denom = 1.0 - (Rs * Rv);
const auto io = pu.phase_pos[Oil];
const auto ig = pu.phase_pos[Gas];
disgas[seg] = Rs * (qs[io] - Rv*qs[ig]) / denom;
vapoil[seg] = Rv * (qs[ig] - Rs*qs[io]) / denom;
}
// 2) Local condition volume flow rates
{
// Use std::span<> in C++20 and beyond.
const auto rate_start = seg * well_.numPhases();
const auto* surf_rates = segment_rates.data() + rate_start;
auto* resv_rates = segments.phase_resv_rates.data() + rate_start;
well_.rateConverter().calcReservoirVoidageRates
(pvtReg, segment_pressure[seg],
std::max(0.0, Rs),
std::max(0.0, Rv),
temperature, saltConc, surf_rates, resv_rates);
}
// 3) Local condition holdup fractions.
const auto tot_resv =
std::accumulate(segments.phase_resv_rates.begin() + (seg + 0) * well_.numPhases(),
segments.phase_resv_rates.begin() + (seg + 1) * well_.numPhases(),
0.0);
std::transform(segments.phase_resv_rates.begin() + (seg + 0) * well_.numPhases(),
segments.phase_resv_rates.begin() + (seg + 1) * well_.numPhases(),
segments.phase_holdup.begin() + (seg + 0) * well_.numPhases(),
[tot_resv](const auto qr) { return std::clamp(qr / tot_resv, 0.0, 1.0); });
// 4) Local condition flow velocities for segments other than top segment.
if (seg > 0) {
// Possibly poor approximation
// Velocity = Flow rate / cross-sectional area.
// Additionally ignores drift flux.
const auto area = well_.wellEcl().getSegments()
.getFromSegmentNumber(segments.segment_number()[seg]).crossArea();
const auto velocity = (area > 0.0) ? tot_resv / area : 0.0;
std::transform(segments.phase_holdup.begin() + (seg + 0) * well_.numPhases(),
segments.phase_holdup.begin() + (seg + 1) * well_.numPhases(),
segments.phase_velocity.begin() + (seg + 0) * well_.numPhases(),
[velocity](const auto hf) { return (hf > 0.0) ? velocity : 0.0; });
}
// 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::waterPhaseIdx)) {
segments.phase_viscosity[seg * well_.numPhases() + pu.phase_pos[Water]] =
FluidSystem::waterPvt().viscosity(pvtReg, temperature, segment_pressure[seg], saltConc);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
segments.phase_viscosity[seg * well_.numPhases() + pu.phase_pos[Gas]] =
FluidSystem::gasPvt().viscosity(pvtReg, temperature, segment_pressure[seg], Rv, Rvw);
}
}
// Segment flow velocity in top segment.
{
const auto np = well_.numPhases();
auto segVel = [&segments, np](const auto segmentNumber)
{
auto v = 0.0;
const auto* vel = segments.phase_velocity.data() + segmentNumber*np;
for (auto p = 0*np; p < np; ++p) {
if (std::abs(vel[p]) > std::abs(v)) {
v = vel[p];
}
}
return v;
};
const auto seg = 0;
auto maxVel = 0.0;
for (const auto& inlet : mswell.segmentSet()[seg].inletSegments()) {
const auto v = segVel(mswell.segmentNumberToIndex(inlet));
if (std::abs(v) > std::abs(maxVel)) {
maxVel = v;
}
}
std::transform(segments.phase_holdup.begin() + (seg + 0) * well_.numPhases(),
segments.phase_holdup.begin() + (seg + 1) * well_.numPhases(),
segments.phase_velocity.begin() + (seg + 0) * well_.numPhases(),
[maxVel](const auto hf) { return (hf > 0.0) ? maxVel : 0.0; });
}
WellBhpThpCalculator(well_)
.updateThp(rho, [this]() { return well_.wellEcl().alq_value(); },
{FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)},
well_state, deferred_logger);
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
processFractions(const int seg)
{
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
const PhaseUsage& pu = well_.phaseUsage();
std::vector<double> 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::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];
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water];
if (fractions[water_pos] < 0.0) {
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]);
}
fractions[oil_pos] /= (1.0 - fractions[water_pos]);
fractions[water_pos] = 0.0;
}
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
if (fractions[gas_pos] < 0.0) {
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
fractions[pu.phase_pos[Water]] /= (1.0 - fractions[gas_pos]);
}
fractions[oil_pos] /= (1.0 - fractions[gas_pos]);
fractions[gas_pos] = 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::gasPhaseIdx) ) {
fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]);
}
fractions[oil_pos] = 0.0;
}
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
value_[seg][WFrac] = fractions[pu.phase_pos[Water]];
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
value_[seg][GFrac] = fractions[pu.phase_pos[Gas]];
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
volumeFraction(const int seg,
const unsigned compIdx) const
{
if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
return evaluation_[seg][WFrac];
}
if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
return evaluation_[seg][GFrac];
}
// Oil fraction
EvalWell oil_fraction = 1.0;
if (has_wfrac_variable) {
oil_fraction -= evaluation_[seg][WFrac];
}
if (has_gfrac_variable) {
oil_fraction -= evaluation_[seg][GFrac];
}
/* if (has_solvent) {
oil_fraction -= evaluation_[seg][SFrac];
} */
return oil_fraction;
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
volumeFractionScaled(const int seg,
const int comp_idx) const
{
// For reservoir rate control, the distr in well control is used for the
// rate conversion coefficients. For the injection well, only the distr of the injection
// phase is not zero.
const double scale = well_.scalingFactor(well_.ebosCompIdxToFlowCompIdx(comp_idx));
if (scale > 0.) {
return this->volumeFraction(seg, comp_idx) / scale;
}
return this->volumeFraction(seg, comp_idx);
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
surfaceVolumeFraction(const int seg,
const int comp_idx) const
{
EvalWell sum_volume_fraction_scaled = 0.;
for (int idx = 0; idx < well_.numComponents(); ++idx) {
sum_volume_fraction_scaled += this->volumeFractionScaled(seg, idx);
}
assert(sum_volume_fraction_scaled.value() != 0.);
return this->volumeFractionScaled(seg, comp_idx) / sum_volume_fraction_scaled;
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
getSegmentRateUpwinding(const int seg,
const int seg_upwind,
const size_t comp_idx) const
{
// the result will contain the derivative with respect to WQTotal in segment seg,
// and the derivatives with respect to WFrac GFrac in segment seg_upwind.
// the derivative with respect to SPres should be zero.
if (seg == 0 && well_.isInjector()) {
const Well& well = well_.wellEcl();
auto phase = well.getInjectionProperties().injectorType;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
&& Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == comp_idx
&& phase == InjectorType::WATER)
return evaluation_[seg][WQTotal] / well_.scalingFactor(well_.ebosCompIdxToFlowCompIdx(comp_idx));
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
&& Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == comp_idx
&& phase == InjectorType::OIL)
return evaluation_[seg][WQTotal] / well_.scalingFactor(well_.ebosCompIdxToFlowCompIdx(comp_idx));
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
&& Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == comp_idx
&& phase == InjectorType::GAS)
return evaluation_[seg][WQTotal] / well_.scalingFactor(well_.ebosCompIdxToFlowCompIdx(comp_idx));
return 0.0;
}
const EvalWell segment_rate = evaluation_[seg][WQTotal] *
this->volumeFractionScaled(seg_upwind, comp_idx);
assert(segment_rate.derivative(SPres + Indices::numEq) == 0.);
return segment_rate;
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
getSegmentPressure(const int seg) const
{
return evaluation_[seg][SPres];
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
getBhp() const
{
return this->getSegmentPressure(0);
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
getSegmentRate(const int seg,
const int comp_idx) const
{
return evaluation_[seg][WQTotal] * this->volumeFractionScaled(seg, comp_idx);
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
getQs(const int comp_idx) const
{
return this->getSegmentRate(0, comp_idx);
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
getWQTotal() const
{
return evaluation_[0][WQTotal];
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
updateUpwindingSegments(const MultisegmentWellGeneric<Scalar>& mswell,
std::vector<int>& upwinding_segments) const
{
for (size_t seg = 0; seg < value_.size(); ++seg) {
// special treatment is needed for segment 0
if (seg == 0) {
// we are not supposed to have injecting producers and producing injectors
assert(!(well_.isProducer() && evaluation_[seg][WQTotal] > 0.));
assert(!(well_.isInjector() && evaluation_[seg][WQTotal] < 0.));
upwinding_segments[seg] = seg;
continue;
}
// for other normal segments
if (evaluation_[seg][WQTotal] <= 0.) {
upwinding_segments[seg] = seg;
} else {
const int outlet_segment_index = mswell.segmentNumberToIndex(mswell.segmentSet()[seg].outletSegment());
upwinding_segments[seg] = outlet_segment_index;
}
}
}
#define INSTANCE(...) \
template class MultisegmentWellPrimaryVariables<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;
// One phase
INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>)
INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>)
INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>)
// Two phase
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,0u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>)
// Blackoil
INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,false,0u,0u>)
INSTANCE(BlackOilIndices<1u,0u,0u,0u,false,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,1u,0u,0u,false,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,1u,0u,false,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,0u,true,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,true,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,false,1u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,2u,0u>)
}

165
opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp Normal file
View File

@@ -0,0 +1,165 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_MULTISEGMENTWELL_PRIMARY_VARIABLES_HEADER_INCLUDED
#define OPM_MULTISEGMENTWELL_PRIMARY_VARIABLES_HEADER_INCLUDED
#include <opm/material/densead/Evaluation.hpp>
#include <opm/simulators/wells/MultisegmentWellEquations.hpp>
#include <array>
#include <vector>
namespace Opm
{
class DeferredLogger;
template<class Scalar> class MultisegmentWellGeneric;
template<class FluidSystem, class Indices, class Scalar> class WellInterfaceIndices;
class WellState;
template<class FluidSystem, class Indices, class Scalar>
class MultisegmentWellPrimaryVariables
{
public:
// TODO: for now, not considering the polymer, solvent and so on to simplify the development process.
// TODO: we need to have order for the primary variables and also the order for the well equations.
// sometimes, they are similar, while sometimes, they can have very different forms.
// Table showing the primary variable indices, depending on what phases are present:
//
// 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
// 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 int WQTotal = 0;
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;
// the number of well equations TODO: it should have a more general strategy for it
static constexpr int numWellEq = Indices::numPhases + 1;
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)
{}
//! \brief Resize values and evaluations.
void resize(const int numSegments);
//! \brief Initialize evaluations from values.
void init();
//! \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);
//! \brief Copy values to well state.
void copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
const double rho,
WellState& well_state,
DeferredLogger& deferred_logger) const;
//! \brief Returns scaled volume fraction for a component in a segment.
//! \details F_p / g_p, the basic usage of this value is because Q_p = G_t * F_p / G_p
EvalWell volumeFractionScaled(const int seg,
const int compIdx) const;
//! \brief Returns surface volume fraction for a component in a segment.
//! \details basically Q_p / \sigma_p Q_p
EvalWell surfaceVolumeFraction(const int seg,
const int compIdx) const;
//! \brief Returns upwinding rate for a component in a segment.
EvalWell getSegmentRateUpwinding(const int seg,
const int seg_upwind,
const size_t comp_idx) const;
//! \brief Get bottomhole pressure.
EvalWell getBhp() const;
//! \brief Get pressure for a segment.
EvalWell getSegmentPressure(const int seg) const;
//! \brief Get rate for a component in a segment.
EvalWell getSegmentRate(const int seg,
const int comp_idx) const;
//! \brief Returns scaled rate for a component.
EvalWell getQs(const int comp_idx) const;
//! \brief Get WQTotal.
EvalWell getWQTotal() const;
//! \brief Update upwinding segments.
void updateUpwindingSegments(const MultisegmentWellGeneric<Scalar>& mswell,
std::vector<int>& upwinding_segments) const;
//! \brief Returns a const ref to an evaluation.
const std::array<EvalWell,numWellEq>& eval(const int idx) const
{ return evaluation_[idx]; }
private:
//! \brief Handle non-reasonable fractions due to numerical overshoot.
void processFractions(const int seg);
//! \brief Returns volume fraction for component in a segment.
EvalWell volumeFraction(const int seg,
const unsigned compIdx) const;
//! \brief The values for the primary variables
//! \details Based on different solution strategies, the wells can have different primary variables
std::vector<std::array<double, numWellEq>> value_;
//! \brief The Evaluation for the well primary variables.
//! \details Contains derivatives and are used in AD calculation
std::vector<std::array<EvalWell, numWellEq>> evaluation_;
const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well_; //!< Reference to well interface
};
}
#endif // OPM_MULTISEGMENTWELL_PRIMARY_VARIABLES_HEADER_INCLUDED

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

@@ -130,7 +130,7 @@ namespace Opm
MultisegmentWell<TypeTag>::
initPrimaryVariablesEvaluation()
{
this->MSWEval::initPrimaryVariablesEvaluation();
this->primary_variables_.init();
}
@@ -142,7 +142,7 @@ namespace Opm
MultisegmentWell<TypeTag>::
updatePrimaryVariables(const WellState& well_state, DeferredLogger& /* deferred_logger */)
{
this->MSWEval::updatePrimaryVariables(well_state);
this->primary_variables_.update(well_state);
}
@@ -453,7 +453,7 @@ namespace Opm
well_flux.clear();
well_flux.resize(np, 0.0);
for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
const EvalWell rate = well_copy.getQs(compIdx);
const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
well_flux[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
}
debug_cost_counter_ += well_copy.debug_cost_counter_;
@@ -595,7 +595,7 @@ namespace Opm
// TODO: trying to reduce the times for the surfaceVolumeFraction calculation
const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value();
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
segment_fluid_initial_[seg][comp_idx] = surface_volume * this->surfaceVolumeFraction(seg, comp_idx).value();
segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value();
}
}
}
@@ -616,12 +616,13 @@ 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);
this->primary_variables_.copyToWellState(*this, getRefDensity(),
well_state, deferred_logger);
Base::calculateReservoirRates(well_state.well(this->index_of_well_));
}
@@ -745,7 +746,7 @@ namespace Opm
this->SPres,
well_state);
}
template<typename TypeTag>
template<class Value>
@@ -919,7 +920,7 @@ namespace Opm
std::vector<EvalWell> cmix_s(this->numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
cmix_s[comp_idx] = this->surfaceVolumeFraction(seg, comp_idx);
cmix_s[comp_idx] = this->primary_variables_.surfaceVolumeFraction(seg, comp_idx);
}
this->computePerfRate(pressure_cell,
@@ -976,7 +977,7 @@ namespace Opm
std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
cmix_s[comp_idx] = getValue(this->surfaceVolumeFraction(seg, comp_idx));
cmix_s[comp_idx] = getValue(this->primary_variables_.surfaceVolumeFraction(seg, comp_idx));
}
Scalar perf_dis_gas_rate = 0.0;
@@ -1525,11 +1526,11 @@ namespace Opm
const GroupState& group_state,
DeferredLogger& deferred_logger)
{
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
// update the upwinding segments
this->updateUpwindingSegments();
this->primary_variables_.updateUpwindingSegments(*this,
this->upwinding_segments_);
// calculate the fluid properties needed.
computeSegmentFluidProperties(ebosSimulator, deferred_logger);
@@ -1563,7 +1564,7 @@ namespace Opm
const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
// for each component
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx)
const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx)
- segment_fluid_initial_[seg][comp_idx]) / dt;
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
assembleAccumulationTerm(seg, comp_idx, accumulation_term, this->linSys_);
@@ -1573,8 +1574,11 @@ namespace Opm
{
const int seg_upwind = this->upwinding_segments_[seg];
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) *
this->well_efficiency_factor_;
const EvalWell segment_rate =
this->primary_variables_.getSegmentRateUpwinding(seg,
seg_upwind,
comp_idx) *
this->well_efficiency_factor_;
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_);
}
@@ -1583,9 +1587,13 @@ namespace Opm
// considering the contributions from the inlet segments
{
for (const int inlet : this->segment_inlets_[seg]) {
const int inlet_upwind = this->upwinding_segments_[inlet];
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const EvalWell inlet_rate = this->getSegmentRateUpwinding(inlet, comp_idx) * this->well_efficiency_factor_;
const int inlet_upwind = this->upwinding_segments_[inlet];
const EvalWell inlet_rate =
this->primary_variables_.getSegmentRateUpwinding(inlet,
inlet_upwind,
comp_idx) *
this->well_efficiency_factor_;
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_);
}
@@ -1593,7 +1601,7 @@ namespace Opm
}
// calculating the perforation rate for each perforation that belongs to this segment
const EvalWell seg_pressure = this->getSegmentPressure(seg);
const EvalWell seg_pressure = this->primary_variables_.getSegmentPressure(seg);
auto& perf_data = ws.perf_data;
auto& perf_rates = perf_data.phase_rates;
auto& perf_press_state = perf_data.pressure;
@@ -1637,7 +1645,6 @@ namespace Opm
if (seg == 0) { // top segment, pressure equation is the control equation
const auto& summaryState = ebosSimulator.vanguard().summaryState();
const Schedule& schedule = ebosSimulator.vanguard().schedule();
std::function<EvalWell(const int)> gQ = [this](int a) { return this->getQs(a); };
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
assembleControlEq(well_state,
group_state,
@@ -1646,9 +1653,7 @@ namespace Opm
inj_controls,
prod_controls,
getRefDensity(),
this->getWQTotal(),
this->getBhp(),
gQ,
this->primary_variables_,
this->linSys_,
deferred_logger);
} else {
@@ -1681,7 +1686,7 @@ namespace Opm
const int nseg = this->numberOfSegments();
for (int seg = 0; seg < nseg; ++seg) {
const EvalWell segment_pressure = this->getSegmentPressure(seg);
const EvalWell segment_pressure = this->primary_variables_.getSegmentPressure(seg);
for (const int perf : this->segment_perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
@@ -1910,7 +1915,7 @@ namespace Opm
const int nseg = this->numberOfSegments();
for (int seg = 0; seg < nseg; ++seg) {
// calculating the perforation rate for each perforation that belongs to this segment
const Scalar seg_pressure = getValue(this->getSegmentPressure(seg));
const Scalar seg_pressure = getValue(this->primary_variables_.getSegmentPressure(seg));
for (const int perf : this->segment_perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));