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 #4279 from akva2/stdwell_primaries

Browse Source 
added: StandardWellPrimaryVariables
This commit is contained in:
Bård Skaflestad 2022-11-22 13:54:00 +01:00 committed by GitHub
commit ab3da9815d
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 992 additions and 842 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -100,6 +100,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/StandardWellEquations.cpp
opm/simulators/wells/StandardWellEval.cpp
opm/simulators/wells/StandardWellGeneric.cpp
opm/simulators/wells/StandardWellPrimaryVariables.cpp
opm/simulators/wells/TargetCalculator.cpp
opm/simulators/wells/VFPHelpers.cpp
opm/simulators/wells/VFPProdProperties.cpp
@ -386,6 +387,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/StandardWellEquations.hpp
opm/simulators/wells/StandardWellEval.hpp
opm/simulators/wells/StandardWellGeneric.hpp
opm/simulators/wells/StandardWellPrimaryVariables.hpp
opm/simulators/wells/TargetCalculator.hpp
opm/simulators/wells/VFPHelpers.hpp
opm/simulators/wells/VFPInjProperties.hpp

7
opm/simulators/wells/StandardWell.hpp
View File

@ -109,6 +109,8 @@ namespace Opm
// TODO: we should have indices for the well equations and well primary variables separately
static constexpr int Bhp = numStaticWellEq - numWellControlEq;
using StdWellEval::WQTotal;
using typename Base::Scalar;
@ -122,7 +124,6 @@ namespace Opm
using Eval = typename StdWellEval::Eval;
using EvalWell = typename StdWellEval::EvalWell;
using BVectorWell = typename StdWellEval::BVectorWell;
using DiagMatrixBlockWellType = typename StdWellEval::DiagMatrixBlockWellType;
StandardWell(const Well& well,
const ParallelWellInfo& pw_info,
@ -356,10 +357,6 @@ namespace Opm
const WellState& well_state,
DeferredLogger& deferred_logger) const;
// update extra primary vriables if there are any
void updateExtraPrimaryVariables(const BVectorWell& dwells) const;
void updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger) const;
virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,

595
opm/simulators/wells/StandardWellEval.cpp
View File

@ -51,531 +51,31 @@ StandardWellEval<FluidSystem,Indices,Scalar>::
StandardWellEval(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif)
: StandardWellGeneric<Scalar>(baseif)
, baseif_(baseif)
, primary_variables_(baseif_)
, F0_(numWellConservationEq)
, linSys_(baseif_.parallelWellInfo())
{
}
template<class FluidSystem, class Indices, class Scalar>
void StandardWellEval<FluidSystem,Indices,Scalar>::
initPrimaryVariablesEvaluation() const
{
for (int eqIdx = 0; eqIdx < numWellEq_; ++eqIdx) {
primary_variables_evaluation_[eqIdx] =
EvalWell::createVariable(numWellEq_ + Indices::numEq, primary_variables_[eqIdx], Indices::numEq + eqIdx);
}
}
template<class FluidSystem, class Indices, class Scalar>
typename StandardWellEval<FluidSystem,Indices,Scalar>::EvalWell
StandardWellEval<FluidSystem,Indices,Scalar>::
extendEval(const Eval& in) const
{
EvalWell out(numWellEq_ + Indices::numEq, in.value());
EvalWell out(primary_variables_.numWellEq() + Indices::numEq, in.value());
for(int eqIdx = 0; eqIdx < Indices::numEq;++eqIdx) {
out.setDerivative(eqIdx, in.derivative(eqIdx));
}
return out;
}
template<class FluidSystem, class Indices, class Scalar>
double
StandardWellEval<FluidSystem,Indices,Scalar>::
relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
const BVectorWell& dwells)
{
// TODO: not considering solvent yet
// 0.95 is a experimental value, which remains to be optimized
double relaxation_factor = 1.0;
if (FluidSystem::numActivePhases() > 1) {
if constexpr (has_wfrac_variable) {
const double relaxation_factor_w = StandardWellGeneric<Scalar>::
relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]);
relaxation_factor = std::min(relaxation_factor, relaxation_factor_w);
}
if constexpr (has_gfrac_variable) {
const double relaxation_factor_g = StandardWellGeneric<Scalar>::
relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]);
relaxation_factor = std::min(relaxation_factor, relaxation_factor_g);
}
if constexpr (has_wfrac_variable && has_gfrac_variable) {
// We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will
// not be negative oil fraction later
const double original_sum = primary_variables[WFrac] + primary_variables[GFrac];
const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor;
const double possible_updated_sum = original_sum - relaxed_update;
// We only relax if fraction is above 1.
// The newton solver should handle the rest
const double epsilon = 0.001;
if (possible_updated_sum > 1.0 + epsilon) {
// since the orignal sum <= 1.0 the epsilon asserts that
// the relaxed_update is non trivial.
assert(relaxed_update != 0.);
const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95;
relaxation_factor *= further_relaxation_factor;
}
}
assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0);
}
return relaxation_factor;
}
template<class FluidSystem, class Indices, class Scalar>
typename StandardWellEval<FluidSystem,Indices,Scalar>::EvalWell
StandardWellEval<FluidSystem,Indices,Scalar>::
wellVolumeFraction(const unsigned compIdx) const
{
if (FluidSystem::numActivePhases() == 1) {
return EvalWell(numWellEq_ + Indices::numEq, 1.0);
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
return primary_variables_evaluation_[WFrac];
}
if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
return primary_variables_evaluation_[GFrac];
}
if (Indices::enableSolvent && compIdx == (unsigned)Indices::contiSolventEqIdx) {
return primary_variables_evaluation_[SFrac];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
return primary_variables_evaluation_[GFrac];
}
}
// Oil or WATER fraction
EvalWell well_fraction(numWellEq_ + Indices::numEq, 1.0);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
well_fraction -= primary_variables_evaluation_[WFrac];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
well_fraction -= primary_variables_evaluation_[GFrac];
}
if (Indices::enableSolvent) {
well_fraction -= primary_variables_evaluation_[SFrac];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))) {
well_fraction -= primary_variables_evaluation_[GFrac];
}
return well_fraction;
}
template<class FluidSystem, class Indices, class Scalar>
typename StandardWellEval<FluidSystem,Indices,Scalar>::EvalWell
StandardWellEval<FluidSystem,Indices,Scalar>::
getQs(const int comp_idx) const
{
// Note: currently, the WQTotal definition is still depends on Injector/Producer.
assert(comp_idx < baseif_.numComponents());
if (baseif_.isInjector()) { // only single phase injection
double inj_frac = 0.0;
switch (baseif_.wellEcl().injectorType()) {
case InjectorType::WATER:
if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx))) {
inj_frac = 1.0;
}
break;
case InjectorType::GAS:
if (Indices::enableSolvent && comp_idx == Indices::contiSolventEqIdx) { // solvent
inj_frac = baseif_.wsolvent();
} else if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx))) {
inj_frac = 1.0 - baseif_.rsRvInj();
if (Indices::enableSolvent) {
inj_frac -= baseif_.wsolvent();
}
} else if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx))) {
inj_frac = baseif_.rsRvInj();
}
break;
case InjectorType::OIL:
if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx))) {
inj_frac = 1.0 - baseif_.rsRvInj();
} else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx))) {
inj_frac = baseif_.rsRvInj();
}
break;
case InjectorType::MULTI:
// Not supported.
// deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED",
// "Multi phase injectors are not supported, requested for well " + name());
break;
}
return inj_frac * primary_variables_evaluation_[WQTotal];
} else { // producers
return primary_variables_evaluation_[WQTotal] * wellVolumeFractionScaled(comp_idx);
}
}
template<class FluidSystem, class Indices, class Scalar>
typename StandardWellEval<FluidSystem,Indices,Scalar>::EvalWell
StandardWellEval<FluidSystem,Indices,Scalar>::
wellVolumeFractionScaled(const int compIdx) const
{
const int legacyCompIdx = baseif_.ebosCompIdxToFlowCompIdx(compIdx);
const double scal = baseif_.scalingFactor(legacyCompIdx);
if (scal > 0)
return wellVolumeFraction(compIdx) / scal;
// the scaling factor may be zero for RESV controlled wells.
return wellVolumeFraction(compIdx);
}
template<class FluidSystem, class Indices, class Scalar>
typename StandardWellEval<FluidSystem,Indices,Scalar>::EvalWell
StandardWellEval<FluidSystem,Indices,Scalar>::
wellSurfaceVolumeFraction(const int compIdx) const
{
EvalWell sum_volume_fraction_scaled(numWellEq_ + Indices::numEq, 0.);
for (int idx = 0; idx < baseif_.numComponents(); ++idx) {
sum_volume_fraction_scaled += wellVolumeFractionScaled(idx);
}
assert(sum_volume_fraction_scaled.value() != 0.);
return wellVolumeFractionScaled(compIdx) / sum_volume_fraction_scaled;
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_logger) const
{
static constexpr int Gas = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Gas;
static constexpr int Oil = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Oil;
static constexpr int Water = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Water;
if (!baseif_.isOperableAndSolvable() && !baseif_.wellIsStopped()) return;
const int well_index = baseif_.indexOfWell();
const int np = baseif_.numPhases();
const auto& pu = baseif_.phaseUsage();
const auto& ws = well_state.well(well_index);
// the weighted total well rate
double total_well_rate = 0.0;
for (int p = 0; p < np; ++p) {
total_well_rate += baseif_.scalingFactor(p) * ws.surface_rates[p];
}
// Not: for the moment, the first primary variable for the injectors is not G_total. The injection rate
// under surface condition is used here
if (baseif_.isInjector()) {
switch (baseif_.wellEcl().injectorType()) {
case InjectorType::WATER:
primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Water]];
break;
case InjectorType::GAS:
primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Gas]];
break;
case InjectorType::OIL:
primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Oil]];
break;
case InjectorType::MULTI:
// Not supported.
deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED",
"Multi phase injectors are not supported, requested for well " + baseif_.name());
break;
}
} else {
primary_variables_[WQTotal] = total_well_rate;
}
if (std::abs(total_well_rate) > 0.) {
if constexpr (has_wfrac_variable) {
primary_variables_[WFrac] = baseif_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate;
}
if constexpr (has_gfrac_variable) {
primary_variables_[GFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * (ws.surface_rates[pu.phase_pos[Gas]]
- (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0) ) / total_well_rate ;
}
if constexpr (Indices::enableSolvent) {
primary_variables_[SFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * ws.sum_solvent_rates() / total_well_rate ;
}
} else { // total_well_rate == 0
if (baseif_.isInjector()) {
// only single phase injection handled
if constexpr (has_wfrac_variable) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
auto phase = baseif_.wellEcl().getInjectionProperties().injectorType;
if (phase == InjectorType::WATER) {
primary_variables_[WFrac] = 1.0;
} else {
primary_variables_[WFrac] = 0.0;
}
}
}
if constexpr (has_gfrac_variable) {
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
auto phase = baseif_.wellEcl().getInjectionProperties().injectorType;
if (phase == InjectorType::GAS) {
primary_variables_[GFrac] = (1.0 - baseif_.rsRvInj());
if constexpr (Indices::enableSolvent) {
primary_variables_[GFrac] = 1.0 - baseif_.rsRvInj() - baseif_.wsolvent();
primary_variables_[SFrac] = baseif_.wsolvent();
}
} else {
primary_variables_[GFrac] = 0.0;
}
}
}
// TODO: it is possible to leave injector as a oil well,
// when F_w and F_g both equals to zero, not sure under what kind of circumstance
// this will happen.
} else if (baseif_.isProducer()) { // producers
// TODO: the following are not addressed for the solvent case yet
if constexpr (has_wfrac_variable) {
primary_variables_[WFrac] = 1.0 / np;
}
if constexpr (has_gfrac_variable) {
primary_variables_[GFrac] = 1.0 / np;
}
} else {
OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
}
}
// BHP
primary_variables_[Bhp] = ws.bhp;
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
updatePrimaryVariablesPolyMW(const BVectorWell& dwells) const
{
if (baseif_.isInjector()) {
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
const int wat_vel_index = Bhp + 1 + perf;
const int pskin_index = Bhp + 1 + baseif_.numPerfs() + perf;
const double relaxation_factor = 0.9;
const double dx_wat_vel = dwells[0][wat_vel_index];
primary_variables_[wat_vel_index] -= relaxation_factor * dx_wat_vel;
const double dx_pskin = dwells[0][pskin_index];
primary_variables_[pskin_index] -= relaxation_factor * dx_pskin;
}
}
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
processFractions() const
{
static constexpr int Gas = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Gas;
static constexpr int Oil = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Oil;
static constexpr int Water = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Water;
const auto pu = baseif_.phaseUsage();
std::vector<double> F(baseif_.numPhases(), 0.0);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
F[pu.phase_pos[Oil]] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] = primary_variables_[WFrac];
F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] = primary_variables_[GFrac];
F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] = primary_variables_[GFrac];
F[pu.phase_pos[Water]] -= F[pu.phase_pos[Gas]];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] = 1.0;
}
[[maybe_unused]] double F_solvent;
if constexpr (Indices::enableSolvent) {
F_solvent = primary_variables_[SFrac];
F[pu.phase_pos[Oil]] -= F_solvent;
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if (F[Water] < 0.0) {
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]);
}
if constexpr (Indices::enableSolvent) {
F_solvent /= (1.0 - F[pu.phase_pos[Water]]);
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]);
}
F[pu.phase_pos[Water]] = 0.0;
}
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (F[pu.phase_pos[Gas]] < 0.0) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]);
}
if constexpr (Indices::enableSolvent) {
F_solvent /= (1.0 - F[pu.phase_pos[Gas]]);
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]);
}
F[pu.phase_pos[Gas]] = 0.0;
}
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
if (F[pu.phase_pos[Oil]] < 0.0) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]);
}
if constexpr (Indices::enableSolvent) {
F_solvent /= (1.0 - F[pu.phase_pos[Oil]]);
}
F[pu.phase_pos[Oil]] = 0.0;
}
}
if constexpr (has_wfrac_variable) {
primary_variables_[WFrac] = F[pu.phase_pos[Water]];
}
if constexpr (has_gfrac_variable) {
primary_variables_[GFrac] = F[pu.phase_pos[Gas]];
}
if constexpr (Indices::enableSolvent) {
primary_variables_[SFrac] = F_solvent;
}
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
updateWellStateFromPrimaryVariables(WellState& well_state,
DeferredLogger& deferred_logger) const
{
static constexpr int Gas = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Gas;
static constexpr int Oil = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Oil;
static constexpr int Water = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Water;
const PhaseUsage& pu = baseif_.phaseUsage();
std::vector<double> F(baseif_.numPhases(), 0.0);
[[maybe_unused]] double F_solvent = 0.0;
if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ) {
const int oil_pos = pu.phase_pos[Oil];
F[oil_pos] = 1.0;
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
const int water_pos = pu.phase_pos[Water];
F[water_pos] = primary_variables_[WFrac];
F[oil_pos] -= F[water_pos];
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
const int gas_pos = pu.phase_pos[Gas];
F[gas_pos] = primary_variables_[GFrac];
F[oil_pos] -= F[gas_pos];
}
if constexpr (Indices::enableSolvent) {
F_solvent = primary_variables_[SFrac];
F[oil_pos] -= F_solvent;
}
}
else if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
const int water_pos = pu.phase_pos[Water];
F[water_pos] = 1.0;
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
const int gas_pos = pu.phase_pos[Gas];
F[gas_pos] = primary_variables_[GFrac];
F[water_pos] -= F[gas_pos];
}
}
else if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
const int gas_pos = pu.phase_pos[Gas];
F[gas_pos] = 1.0;
}
// convert the fractions to be Q_p / G_total to calculate the phase rates
for (int p = 0; p < baseif_.numPhases(); ++p) {
const double scal = baseif_.scalingFactor(p);
// for injection wells, there should only one non-zero scaling factor
if (scal > 0) {
F[p] /= scal ;
} else {
// this should only happens to injection wells
F[p] = 0.;
}
}
// F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent.
// More testing is needed to make sure this is correct for well groups and THP.
if constexpr (Indices::enableSolvent){
F_solvent /= baseif_.scalingFactor(Indices::contiSolventEqIdx);
F[pu.phase_pos[Gas]] += F_solvent;
}
auto& ws = well_state.well(baseif_.indexOfWell());
ws.bhp = primary_variables_[Bhp];
// calculate the phase rates based on the primary variables
// for producers, this is not a problem, while not sure for injectors here
if (baseif_.isProducer()) {
const double g_total = primary_variables_[WQTotal];
for (int p = 0; p < baseif_.numPhases(); ++p) {
ws.surface_rates[p] = g_total * F[p];
}
} else { // injectors
for (int p = 0; p < baseif_.numPhases(); ++p) {
ws.surface_rates[p] = 0.0;
}
switch (baseif_.wellEcl().injectorType()) {
case InjectorType::WATER:
ws.surface_rates[pu.phase_pos[Water]] = primary_variables_[WQTotal];
break;
case InjectorType::GAS:
ws.surface_rates[pu.phase_pos[Gas]] = primary_variables_[WQTotal];
break;
case InjectorType::OIL:
ws.surface_rates[pu.phase_pos[Oil]] = primary_variables_[WQTotal];
break;
case InjectorType::MULTI:
// Not supported.
deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED",
"Multi phase injectors are not supported, requested for well " + baseif_.name());
break;
}
}
this->primary_variables_.copyToWellState(well_state, deferred_logger);
WellBhpThpCalculator(baseif_).
updateThp(this->getRho(),
@ -586,81 +86,13 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
well_state, deferred_logger);
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
updateWellStateFromPrimaryVariablesPolyMW(WellState& well_state) const
{
if (baseif_.isInjector()) {
auto& ws = well_state.well(baseif_.indexOfWell());
auto& perf_data = ws.perf_data;
auto& perf_water_velocity = perf_data.water_velocity;
auto& perf_skin_pressure = perf_data.skin_pressure;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
perf_water_velocity[perf] = primary_variables_[Bhp + 1 + perf];
perf_skin_pressure[perf] = primary_variables_[Bhp + 1 + baseif_.numPerfs() + perf];
}
}
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
computeAccumWell()
{
for (int eq_idx = 0; eq_idx < numWellConservationEq; ++eq_idx) {
F0_[eq_idx] = wellSurfaceVolumeFraction(eq_idx).value();
}
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
updatePrimaryVariablesNewton(const BVectorWell& dwells,
[[maybe_unused]] const double dFLimit,
const double dBHPLimit) const
{
const std::vector<double> old_primary_variables = primary_variables_;
// for injectors, very typical one of the fractions will be one, and it is easy to get zero value
// fractions. not sure what is the best way to handle it yet, so we just use 1.0 here
[[maybe_unused]] const double relaxation_factor_fractions =
(baseif_.isProducer()) ? relaxationFactorFractionsProducer(old_primary_variables, dwells)
: 1.0;
// update the second and third well variable (The flux fractions)
if constexpr (has_wfrac_variable) {
const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit);
// primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
}
if constexpr (has_gfrac_variable) {
const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
}
if constexpr (Indices::enableSolvent) {
const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
}
processFractions();
// updating the total rates Q_t
const double relaxation_factor_rate = this->relaxationFactorRate(old_primary_variables, dwells[0][WQTotal]);
primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
// updating the bottom hole pressure
{
const int sign1 = dwells[0][Bhp] > 0 ? 1: -1;
const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]), std::abs(old_primary_variables[Bhp]) * dBHPLimit);
// 1e5 to make sure bhp will not be below 1bar
primary_variables_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5);
for (size_t eq_idx = 0; eq_idx < F0_.size(); ++eq_idx) {
F0_[eq_idx] = this->primary_variables_.surfaceVolumeFraction(eq_idx).value();
}
}
@ -676,8 +108,8 @@ getWellConvergence(const WellState& well_state,
std::vector<double>& res,
DeferredLogger& deferred_logger) const
{
res.resize(numWellEq_);
for (int eq_idx = 0; eq_idx < numWellEq_; ++eq_idx) {
res.resize(this->primary_variables_.numWellEq());
for (int eq_idx = 0; eq_idx < this->primary_variables_.numWellEq(); ++eq_idx) {
// magnitude of the residual matters
res[eq_idx] = std::abs(this->linSys_.residual()[0][eq_idx]);
}
@ -931,22 +363,21 @@ init(std::vector<double>& perf_depth,
}
// counting/updating primary variable numbers
int numWellEq = primary_variables_.numWellEq();
if (has_polymermw) {
if (baseif_.isInjector()) {
// adding a primary variable for water perforation rate per connection
numWellEq_ += baseif_.numPerfs();
numWellEq += baseif_.numPerfs();
// adding a primary variable for skin pressure per connection
numWellEq_ += baseif_.numPerfs();
numWellEq += baseif_.numPerfs();
}
}
// with the updated numWellEq_, we can initialize the primary variables and matrices now
primary_variables_.resize(numWellEq_, 0.0);
primary_variables_evaluation_.resize(numWellEq_, EvalWell{numWellEq_ + Indices::numEq, 0.0});
// with the updated numWellEq, we can initialize the primary variables and matrices now
primary_variables_.resize(numWellEq);
// setup sparsity pattern for the matrices
this->linSys_.init(num_cells, this->numWellEq_,
baseif_.numPerfs(), baseif_.cells());
this->linSys_.init(num_cells, numWellEq, baseif_.numPerfs(), baseif_.cells());
}
#define INSTANCE(...) \

109
opm/simulators/wells/StandardWellEval.hpp
View File

@ -25,6 +25,7 @@
#include <opm/simulators/wells/StandardWellEquations.hpp>
#include <opm/simulators/wells/StandardWellGeneric.hpp>
#include <opm/simulators/wells/StandardWellPrimaryVariables.hpp>
#include <opm/material/densead/DynamicEvaluation.hpp>
@ -47,59 +48,21 @@ template<class FluidSystem, class Indices, class Scalar>
class StandardWellEval : public StandardWellGeneric<Scalar>
{
protected:
// number of the conservation equations
static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents;
// number of the well control equations
static constexpr int numWellControlEq = 1;
// number of the well equations that will always be used
// based on the solution strategy, there might be other well equations be introduced
using PrimaryVariables = StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>;
static constexpr int Bhp = PrimaryVariables::Bhp;
static constexpr int WQTotal= PrimaryVariables::WQTotal;
static constexpr int numWellConservationEq = PrimaryVariables::numWellConservationEq;
static constexpr bool has_wfrac_variable = PrimaryVariables::has_wfrac_variable;
static constexpr bool has_gfrac_variable = PrimaryVariables::has_gfrac_variable;
static constexpr int WFrac = PrimaryVariables::WFrac;
static constexpr int GFrac = PrimaryVariables::GFrac;
static constexpr int SFrac = PrimaryVariables::SFrac;
public:
static constexpr int numStaticWellEq = numWellConservationEq + numWellControlEq;
protected:
// the index for Bhp in primary variables and also the index of well control equation
// they both will be the last one in their respective system.
// TODO: we should have indices for the well equations and well primary variables separately
static constexpr int Bhp = numStaticWellEq - numWellControlEq;
// the positions of the primary variables for StandardWell
// the first one is the weighted total rate (WQ_t), the second and the third ones are F_w and F_g,
// which represent the fraction of Water and Gas based on the weighted total rate, the last one is BHP.
// correspondingly, we have four well equations for blackoil model, the first three are mass
// converstation equations, and the last one is the well control equation.
// primary variables related to other components, will be before the Bhp and after F_g.
// well control equation is always the last well equation.
// TODO: in the current implementation, we use the well rate as the first primary variables for injectors,
// instead of G_t.
// 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 -1000 1 -1000
// GFrac 2 1 1 -1000 -1000
// Spres 3 2 2 2 1
static const int WQTotal = 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.
// (following implementation MultisegmentWellEval.hpp)
static const bool waterEnabled = Indices::waterEnabled;
static const bool gasEnabled = Indices::gasEnabled;
static const bool oilEnabled = Indices::oilEnabled;
static constexpr bool has_wfrac_variable = Indices::waterEnabled && Indices::oilEnabled;
static constexpr bool has_gfrac_variable = Indices::gasEnabled && Indices::numPhases > 1;
static constexpr int WFrac = has_wfrac_variable ? 1 : -1000;
static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000;
static constexpr int SFrac = !Indices::enableSolvent ? -1000 : 3;
public:
using EvalWell = DenseAd::DynamicEvaluation<Scalar, numStaticWellEq + Indices::numEq + 1>;
using EvalWell = typename PrimaryVariables::EvalWell;
using Eval = DenseAd::Evaluation<Scalar, Indices::numEq>;
using BVectorWell = typename StandardWellGeneric<Scalar>::BVectorWell;
using BVectorWell = typename StandardWellEquations<Scalar,Indices::numEq>::BVectorWell;
//! \brief Returns a const reference to equation system.
const StandardWellEquations<Scalar,Indices::numEq>& linSys() const
@ -110,28 +73,7 @@ protected:
const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif_;
void initPrimaryVariablesEvaluation() const;
const EvalWell& getBhp() const
{
return primary_variables_evaluation_[Bhp];
}
const EvalWell& getWQTotal() const
{
return primary_variables_evaluation_[WQTotal];
}
EvalWell extendEval(const Eval& in) const;
EvalWell getQs(const int compIdx) const;
EvalWell wellSurfaceVolumeFraction(const int compIdx) const;
EvalWell wellVolumeFraction(const unsigned compIdx) const;
EvalWell wellVolumeFractionScaled(const int phase) const;
// calculate a relaxation factor to avoid overshoot of the fractions for producers
// which might result in negative rates
static double relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
const BVectorWell& dwells);
// computing the accumulation term for later use in well mass equations
void computeAccumWell();
@ -160,33 +102,10 @@ protected:
const int num_cells,
const bool has_polymermw);
// handle the non reasonable fractions due to numerical overshoot
void processFractions() const;
void updatePrimaryVariables(const WellState& well_state,
DeferredLogger& deferred_logger) const;
void updatePrimaryVariablesPolyMW(const BVectorWell& dwells) const;
void updateWellStateFromPrimaryVariables(WellState& well_state,
DeferredLogger& deferred_logger) const;
void updatePrimaryVariablesNewton(const BVectorWell& dwells,
const double dFLimit,
const double dBHPLimit) const;
void updateWellStateFromPrimaryVariablesPolyMW(WellState& well_state) const;
// total number of the well equations and primary variables
// there might be extra equations be used, numWellEq will be updated during the initialization
int numWellEq_ = numStaticWellEq;
// the values for the primary varibles
// based on different solutioin strategies, the wells can have different primary variables
mutable std::vector<double> primary_variables_;
// the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
mutable std::vector<EvalWell> primary_variables_evaluation_;
mutable PrimaryVariables primary_variables_; //!< Primary variables for well
// the saturations in the well bore under surface conditions at the beginning of the time step
std::vector<double> F0_;

53
opm/simulators/wells/StandardWellGeneric.cpp
View File

@ -55,59 +55,6 @@ StandardWellGeneric(const WellInterfaceGeneric& baseif)
{
}
template<class Scalar>
double
StandardWellGeneric<Scalar>::
relaxationFactorRate(const std::vector<double>& primary_variables,
const double newton_update)
{
double relaxation_factor = 1.0;
static constexpr int WQTotal = 0;
// For injector, we only check the total rates to avoid sign change of rates
const double original_total_rate = primary_variables[WQTotal];
const double possible_update_total_rate = primary_variables[WQTotal] - newton_update;
// 0.8 here is a experimental value, which remains to be optimized
// if the original rate is zero or possible_update_total_rate is zero, relaxation_factor will
// always be 1.0, more thoughts might be needed.
if (original_total_rate * possible_update_total_rate < 0.) { // sign changed
relaxation_factor = std::abs(original_total_rate / newton_update) * 0.8;
}
assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0);
return relaxation_factor;
}
template<class Scalar>
double
StandardWellGeneric<Scalar>::
relaxationFactorFraction(const double old_value,
const double dx)
{
assert(old_value >= 0. && old_value <= 1.0);
double relaxation_factor = 1.;
// updated values without relaxation factor
const double possible_updated_value = old_value - dx;
// 0.95 is an experimental value remains to be optimized
if (possible_updated_value < 0.0) {
relaxation_factor = std::abs(old_value / dx) * 0.95;
} else if (possible_updated_value > 1.0) {
relaxation_factor = std::abs((1. - old_value) / dx) * 0.95;
}
// if possible_updated_value is between 0. and 1.0, then relaxation_factor
// remains to be one
assert(relaxation_factor >= 0. && relaxation_factor <= 1.);
return relaxation_factor;
}
template<class Scalar>
void
StandardWellGeneric<Scalar>::

24
opm/simulators/wells/StandardWellGeneric.hpp
View File

@ -49,32 +49,8 @@ template<class Scalar>
class StandardWellGeneric
{
protected:
// sparsity pattern for the matrices
//[A C^T [x = [ res
// B D ] x_well] res_well]
// the vector type for the res_well and x_well
using VectorBlockWellType = Dune::DynamicVector<Scalar>;
using BVectorWell = Dune::BlockVector<VectorBlockWellType>;
// the matrix type for the diagonal matrix D
using DiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
// the matrix type for the non-diagonal matrix B and C^T
using OffDiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
StandardWellGeneric(const WellInterfaceGeneric& baseif);
// calculate a relaxation factor to avoid overshoot of total rates
static double relaxationFactorRate(const std::vector<double>& primary_variables,
const double newton_update);
// relaxation factor considering only one fraction value
static double relaxationFactorFraction(const double old_value,
const double dx);
void computeConnectionPressureDelta();
// Base interface reference

732
opm/simulators/wells/StandardWellPrimaryVariables.cpp Normal file
View File

@ -0,0 +1,732 @@
/*
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/StandardWellPrimaryVariables.hpp>
#include <dune/common/dynvector.hh>
#include <dune/istl/bvector.hh>
#include <opm/material/densead/Evaluation.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/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/StandardWellGeneric.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <algorithm>
#include <cassert>
namespace {
//! \brief Relaxation factor considering only one fraction value.
template<class Scalar>
Scalar relaxationFactorFraction(const Scalar old_value,
const Scalar dx)
{
assert(old_value >= 0. && old_value <= 1.0);
Scalar relaxation_factor = 1.;
// updated values without relaxation factor
const Scalar possible_updated_value = old_value - dx;
// 0.95 is an experimental value remains to be optimized
if (possible_updated_value < 0.0) {
relaxation_factor = std::abs(old_value / dx) * 0.95;
} else if (possible_updated_value > 1.0) {
relaxation_factor = std::abs((1. - old_value) / dx) * 0.95;
}
// if possible_updated_value is between 0. and 1.0, then relaxation_factor
// remains to be one
assert(relaxation_factor >= 0. && relaxation_factor <= 1.);
return relaxation_factor;
}
//! \brief Calculate a relaxation factor to avoid overshoot of total rates.
template<class Scalar>
Scalar relaxationFactorRate(const Scalar old_value,
const Scalar newton_update)
{
Scalar relaxation_factor = 1.0;
// For injector, we only check the total rates to avoid sign change of rates
const Scalar original_total_rate = old_value;
const Scalar possible_update_total_rate = old_value - newton_update;
// 0.8 here is a experimental value, which remains to be optimized
// if the original rate is zero or possible_update_total_rate is zero, relaxation_factor will
// always be 1.0, more thoughts might be needed.
if (original_total_rate * possible_update_total_rate < 0.) { // sign changed
relaxation_factor = std::abs(original_total_rate / newton_update) * 0.8;
}
assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0);
return relaxation_factor;
}
}
namespace Opm {
template<class FluidSystem, class Indices, class Scalar>
void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
init()
{
for (int eqIdx = 0; eqIdx < numWellEq_; ++eqIdx) {
evaluation_[eqIdx] =
EvalWell::createVariable(numWellEq_ + Indices::numEq,
value_[eqIdx],
Indices::numEq + eqIdx);
}
}
template<class FluidSystem, class Indices, class Scalar>
void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
resize(const int numWellEq)
{
value_.resize(numWellEq, 0.0);
evaluation_.resize(numWellEq, EvalWell{numWellEq + Indices::numEq, 0.0});
numWellEq_ = numWellEq;
}
template<class FluidSystem, class Indices, class Scalar>
void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
update(const WellState& well_state, DeferredLogger& deferred_logger)
{
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
const int well_index = well_.indexOfWell();
const int np = well_.numPhases();
const auto& pu = well_.phaseUsage();
const auto& ws = well_state.well(well_index);
// the weighted total well rate
double total_well_rate = 0.0;
for (int p = 0; p < np; ++p) {
total_well_rate += well_.scalingFactor(p) * ws.surface_rates[p];
}
// Not: for the moment, the first primary variable for the injectors is not G_total. The injection rate
// under surface condition is used here
if (well_.isInjector()) {
switch (well_.wellEcl().injectorType()) {
case InjectorType::WATER:
value_[WQTotal] = ws.surface_rates[pu.phase_pos[Water]];
break;
case InjectorType::GAS:
value_[WQTotal] = ws.surface_rates[pu.phase_pos[Gas]];
break;
case InjectorType::OIL:
value_[WQTotal] = ws.surface_rates[pu.phase_pos[Oil]];
break;
case InjectorType::MULTI:
// Not supported.
deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED",
"Multi phase injectors are not supported, requested for well " + well_.name());
break;
}
} else {
value_[WQTotal] = total_well_rate;
}
if (std::abs(total_well_rate) > 0.) {
if constexpr (has_wfrac_variable) {
value_[WFrac] = well_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate;
}
if constexpr (has_gfrac_variable) {
value_[GFrac] = well_.scalingFactor(pu.phase_pos[Gas]) *
(ws.surface_rates[pu.phase_pos[Gas]] -
(Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0) ) / total_well_rate ;
}
if constexpr (Indices::enableSolvent) {
value_[SFrac] = well_.scalingFactor(pu.phase_pos[Gas]) * ws.sum_solvent_rates() / total_well_rate ;
}
} else { // total_well_rate == 0
if (well_.isInjector()) {
// only single phase injection handled
if constexpr (has_wfrac_variable) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
auto phase = well_.wellEcl().getInjectionProperties().injectorType;
if (phase == InjectorType::WATER) {
value_[WFrac] = 1.0;
} else {
value_[WFrac] = 0.0;
}
}
}
if constexpr (has_gfrac_variable) {
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
auto phase = well_.wellEcl().getInjectionProperties().injectorType;
if (phase == InjectorType::GAS) {
value_[GFrac] = (1.0 - well_.rsRvInj());
if constexpr (Indices::enableSolvent) {
value_[GFrac] = 1.0 - well_.rsRvInj() - well_.wsolvent();
value_[SFrac] = well_.wsolvent();
}
} else {
value_[GFrac] = 0.0;
}
}
}
// TODO: it is possible to leave injector as a oil well,
// when F_w and F_g both equals to zero, not sure under what kind of circumstance
// this will happen.
} else if (well_.isProducer()) { // producers
// TODO: the following are not addressed for the solvent case yet
if constexpr (has_wfrac_variable) {
value_[WFrac] = 1.0 / np;
}
if constexpr (has_gfrac_variable) {
value_[GFrac] = 1.0 / np;
}
} else {
OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
}
}
// BHP
value_[Bhp] = ws.bhp;
}
template<class FluidSystem, class Indices, class Scalar>
void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
updatePolyMW(const WellState& well_state)
{
if (well_.isInjector()) {
const auto& ws = well_state.well(well_.indexOfWell());
const auto& perf_data = ws.perf_data;
const auto& water_velocity = perf_data.water_velocity;
const auto& skin_pressure = perf_data.skin_pressure;
for (int perf = 0; perf < well_.numPerfs(); ++perf) {
value_[Bhp + 1 + perf] = water_velocity[perf];
value_[Bhp + 1 + well_.numPerfs() + perf] = skin_pressure[perf];
}
}
}
template<class FluidSystem, class Indices, class Scalar>
void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
updateNewton(const BVectorWell& dwells,
[[maybe_unused]] const double dFLimit,
const double dBHPLimit)
{
const double relaxation_factor_rate = relaxationFactorRate(value_[WQTotal],
dwells[0][WQTotal]);
// for injectors, very typical one of the fractions will be one, and it is easy to get zero value
// fractions. not sure what is the best way to handle it yet, so we just use 1.0 here
[[maybe_unused]] const double relaxation_factor_fractions =
well_.isProducer() ? this->relaxationFactorFractionsProducer(dwells) : 1.0;
// update the second and third well variable (The flux fractions)
if constexpr (has_wfrac_variable) {
const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit);
value_[WFrac] = value_[WFrac] - dx2_limited;
}
if constexpr (has_gfrac_variable) {
const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
value_[GFrac] = value_[GFrac] - dx3_limited;
}
if constexpr (Indices::enableSolvent) {
const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
value_[SFrac] = value_[SFrac] - dx4_limited;
}
this->processFractions();
// updating the total rates Q_t
value_[WQTotal] = value_[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
// updating the bottom hole pressure
const int sign1 = dwells[0][Bhp] > 0 ? 1: -1;
const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]),
std::abs(value_[Bhp]) * dBHPLimit);
// 1e5 to make sure bhp will not be below 1bar
value_[Bhp] = std::max(value_[Bhp] - dx1_limited, 1e5);
}
template<class FluidSystem, class Indices, class Scalar>
void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
updateNewtonPolyMW(const BVectorWell& dwells)
{
if (well_.isInjector()) {
for (int perf = 0; perf < well_.numPerfs(); ++perf) {
const int wat_vel_index = Bhp + 1 + perf;
const int pskin_index = Bhp + 1 + well_.numPerfs() + perf;
const double relaxation_factor = 0.9;
const double dx_wat_vel = dwells[0][wat_vel_index];
value_[wat_vel_index] -= relaxation_factor * dx_wat_vel;
const double dx_pskin = dwells[0][pskin_index];
value_[pskin_index] -= relaxation_factor * dx_pskin;
}
}
}
template<class FluidSystem, class Indices, class Scalar>
void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
copyToWellState(WellState& well_state,
DeferredLogger& deferred_logger) const
{
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> F(well_.numPhases(), 0.0);
[[maybe_unused]] double F_solvent = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int oil_pos = pu.phase_pos[Oil];
F[oil_pos] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water];
F[water_pos] = value_[WFrac];
F[oil_pos] -= F[water_pos];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
F[gas_pos] = value_[GFrac];
F[oil_pos] -= F[gas_pos];
}
if constexpr (Indices::enableSolvent) {
F_solvent = value_[SFrac];
F[oil_pos] -= F_solvent;
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water];
F[water_pos] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
F[gas_pos] = value_[GFrac];
F[water_pos] -= F[gas_pos];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
F[gas_pos] = 1.0;
}
// convert the fractions to be Q_p / G_total to calculate the phase rates
for (int p = 0; p < well_.numPhases(); ++p) {
const double scal = well_.scalingFactor(p);
// for injection wells, there should only one non-zero scaling factor
if (scal > 0) {
F[p] /= scal ;
} else {
// this should only happens to injection wells
F[p] = 0.;
}
}
// F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent.
// More testing is needed to make sure this is correct for well groups and THP.
if constexpr (Indices::enableSolvent) {
F_solvent /= well_.scalingFactor(Indices::contiSolventEqIdx);
F[pu.phase_pos[Gas]] += F_solvent;
}
auto& ws = well_state.well(well_.indexOfWell());
ws.bhp = value_[Bhp];
// calculate the phase rates based on the primary variables
// for producers, this is not a problem, while not sure for injectors here
if (well_.isProducer()) {
const double g_total = value_[WQTotal];
for (int p = 0; p < well_.numPhases(); ++p) {
ws.surface_rates[p] = g_total * F[p];
}
} else { // injectors
for (int p = 0; p < well_.numPhases(); ++p) {
ws.surface_rates[p] = 0.0;
}
switch (well_.wellEcl().injectorType()) {
case InjectorType::WATER:
ws.surface_rates[pu.phase_pos[Water]] = value_[WQTotal];
break;
case InjectorType::GAS:
ws.surface_rates[pu.phase_pos[Gas]] = value_[WQTotal];
break;
case InjectorType::OIL:
ws.surface_rates[pu.phase_pos[Oil]] = value_[WQTotal];
break;
case InjectorType::MULTI:
// Not supported.
deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED",
"Multi phase injectors are not supported, requested for well " + well_.name());
break;
}
}
}
template<class FluidSystem, class Indices, class Scalar>
void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
copyToWellStatePolyMW(WellState& well_state) const
{
if (well_.isInjector()) {
auto& ws = well_state.well(well_.indexOfWell());
auto& perf_data = ws.perf_data;
auto& perf_water_velocity = perf_data.water_velocity;
auto& perf_skin_pressure = perf_data.skin_pressure;
for (int perf = 0; perf < well_.numPerfs(); ++perf) {
perf_water_velocity[perf] = value_[Bhp + 1 + perf];
perf_skin_pressure[perf] = value_[Bhp + 1 + well_.numPerfs() + perf];
}
}
}
template<class FluidSystem, class Indices, class Scalar>
typename StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
volumeFraction(const unsigned compIdx) const
{
if (FluidSystem::numActivePhases() == 1) {
return EvalWell(numWellEq_ + Indices::numEq, 1.0);
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
return evaluation_[WFrac];
}
if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
return evaluation_[GFrac];
}
if (Indices::enableSolvent && compIdx == (unsigned)Indices::contiSolventEqIdx) {
return evaluation_[SFrac];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
return evaluation_[GFrac];
}
}
// Oil or WATER fraction
EvalWell well_fraction(numWellEq_ + Indices::numEq, 1.0);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
well_fraction -= evaluation_[WFrac];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
well_fraction -= evaluation_[GFrac];
}
if constexpr (Indices::enableSolvent) {
well_fraction -= evaluation_[SFrac];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
well_fraction -= evaluation_[GFrac];
}
return well_fraction;
}
template<class FluidSystem, class Indices, class Scalar>
typename StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
volumeFractionScaled(const int compIdx) const
{
const int legacyCompIdx = well_.ebosCompIdxToFlowCompIdx(compIdx);
const double scal = well_.scalingFactor(legacyCompIdx);
if (scal > 0)
return this->volumeFraction(compIdx) / scal;
// the scaling factor may be zero for RESV controlled wells.
return this->volumeFraction(compIdx);
}
template<class FluidSystem, class Indices, class Scalar>
typename StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
surfaceVolumeFraction(const int compIdx) const
{
EvalWell sum_volume_fraction_scaled(numWellEq_ + Indices::numEq, 0.);
for (int idx = 0; idx < well_.numComponents(); ++idx) {
sum_volume_fraction_scaled += this->volumeFractionScaled(idx);
}
assert(sum_volume_fraction_scaled.value() != 0.);
return this->volumeFractionScaled(compIdx) / sum_volume_fraction_scaled;
}
template<class FluidSystem, class Indices, class Scalar>
typename StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::EvalWell
StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
getQs(const int comp_idx) const
{
// Note: currently, the WQTotal definition is still depends on Injector/Producer.
assert(comp_idx < well_.numComponents());
if (well_.isInjector()) { // only single phase injection
double inj_frac = 0.0;
switch (well_.wellEcl().injectorType()) {
case InjectorType::WATER:
if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx))) {
inj_frac = 1.0;
}
break;
case InjectorType::GAS:
if (Indices::enableSolvent && comp_idx == Indices::contiSolventEqIdx) { // solvent
inj_frac = well_.wsolvent();
} else if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx))) {
inj_frac = 1.0 - well_.rsRvInj();
if constexpr (Indices::enableSolvent) {
inj_frac -= well_.wsolvent();
}
} else if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx))) {
inj_frac = well_.rsRvInj();
}
break;
case InjectorType::OIL:
if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx))) {
inj_frac = 1.0 - well_.rsRvInj();
} else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx))) {
inj_frac = well_.rsRvInj();
}
break;
case InjectorType::MULTI:
// Not supported.
// deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED",
// "Multi phase injectors are not supported, requested for well " + name());
break;
}
return inj_frac * evaluation_[WQTotal];
} else { // producers
return evaluation_[WQTotal] * this->volumeFractionScaled(comp_idx);
}
}
template<class FluidSystem, class Indices, class Scalar>
void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
processFractions()
{
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
const auto pu = well_.phaseUsage();
std::vector<double> F(well_.numPhases(), 0.0);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
F[pu.phase_pos[Oil]] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] = value_[WFrac];
F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] = value_[GFrac];
F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] = value_[GFrac];
F[pu.phase_pos[Water]] -= F[pu.phase_pos[Gas]];
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] = 1.0;
}
[[maybe_unused]] double F_solvent;
if constexpr (Indices::enableSolvent) {
F_solvent = value_[SFrac];
F[pu.phase_pos[Oil]] -= F_solvent;
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if (F[Water] < 0.0) {
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]);
}
if constexpr (Indices::enableSolvent) {
F_solvent /= (1.0 - F[pu.phase_pos[Water]]);
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]);
}
F[pu.phase_pos[Water]] = 0.0;
}
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (F[pu.phase_pos[Gas]] < 0.0) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]);
}
if constexpr (Indices::enableSolvent) {
F_solvent /= (1.0 - F[pu.phase_pos[Gas]]);
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]);
}
F[pu.phase_pos[Gas]] = 0.0;
}
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
if (F[pu.phase_pos[Oil]] < 0.0) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]);
}
if constexpr (Indices::enableSolvent) {
F_solvent /= (1.0 - F[pu.phase_pos[Oil]]);
}
F[pu.phase_pos[Oil]] = 0.0;
}
}
if constexpr (has_wfrac_variable) {
value_[WFrac] = F[pu.phase_pos[Water]];
}
if constexpr (has_gfrac_variable) {
value_[GFrac] = F[pu.phase_pos[Gas]];
}
if constexpr (Indices::enableSolvent) {
value_[SFrac] = F_solvent;
}
}
template<class FluidSystem, class Indices, class Scalar>
double StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
relaxationFactorFractionsProducer(const BVectorWell& dwells) const
{
// TODO: not considering solvent yet
// 0.95 is a experimental value, which remains to be optimized
double relaxation_factor = 1.0;
if (FluidSystem::numActivePhases() > 1) {
if constexpr (has_wfrac_variable) {
const double relaxation_factor_w = relaxationFactorFraction(value_[WFrac],
dwells[0][WFrac]);
relaxation_factor = std::min(relaxation_factor, relaxation_factor_w);
}
if constexpr (has_gfrac_variable) {
const double relaxation_factor_g = relaxationFactorFraction(value_[GFrac],
dwells[0][GFrac]);
relaxation_factor = std::min(relaxation_factor, relaxation_factor_g);
}
if constexpr (has_wfrac_variable && has_gfrac_variable) {
// We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will
// not be negative oil fraction later
const double original_sum = value_[WFrac] + value_[GFrac];
const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor;
const double possible_updated_sum = original_sum - relaxed_update;
// We only relax if fraction is above 1.
// The newton solver should handle the rest
const double epsilon = 0.001;
if (possible_updated_sum > 1.0 + epsilon) {
// since the orignal sum <= 1.0 the epsilon asserts that
// the relaxed_update is non trivial.
assert(relaxed_update != 0.);
const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95;
relaxation_factor *= further_relaxation_factor;
}
}
assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0);
}
return relaxation_factor;
}
template<class FluidSystem, class Indices, class Scalar>
void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
checkFinite(DeferredLogger& deferred_logger) const
{
for (const Scalar v : value_) {
if (!isfinite(v))
OPM_DEFLOG_THROW(NumericalIssue, "Infinite primary variable after update from wellState well: " << well_.name(), deferred_logger);
}
}
#define INSTANCE(...) \
template class StandardWellPrimaryVariables<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,1u,0u,false,true,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,1u,false,false,1u,0u>)
}

171
opm/simulators/wells/StandardWellPrimaryVariables.hpp Normal file
View File

@ -0,0 +1,171 @@
/*
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/>.
*/
#ifndef OPM_STANDARDWELL_PRIMARY_VARIABLES_HEADER_INCLUDED
#define OPM_STANDARDWELL_PRIMARY_VARIABLES_HEADER_INCLUDED
#include <opm/material/densead/DynamicEvaluation.hpp>
#include <opm/simulators/wells/StandardWellEquations.hpp>
#include <vector>
namespace Opm
{
class DeferredLogger;
template<class FluidSystem, class Indices, class Scalar> class WellInterfaceIndices;
class WellState;
//! \brief Class holding primary variables for StandardWell.
template<class FluidSystem, class Indices, class Scalar>
class StandardWellPrimaryVariables {
protected:
// the positions of the primary variables for StandardWell
// the first one is the weighted total rate (WQ_t), the second and the third ones are F_w and F_g,
// which represent the fraction of Water and Gas based on the weighted total rate, the last one is BHP.
// correspondingly, we have four well equations for blackoil model, the first three are mass
// converstation equations, and the last one is the well control equation.
// primary variables related to other components, will be before the Bhp and after F_g.
// well control equation is always the last well equation.
// TODO: in the current implementation, we use the well rate as the first primary variables for injectors,
// instead of G_t.
// 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 -1000 1 -1000
// GFrac 2 1 1 -1000 -1000
// Spres 3 2 2 2 1
//! \brief Number of the well control equations.
static constexpr int numWellControlEq = 1;
public:
//! \brief Number of the conservation equations.
static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents;
//! \brief Number of the well equations that will always be used.
//! \details Based on the solution strategy, there might be other well equations be introduced.
static constexpr int numStaticWellEq = numWellConservationEq + numWellControlEq;
static constexpr int WQTotal = 0; //!< The index for the weighted total rate
//! \brief The index for Bhp in primary variables and the index of well control equation.
//! \details They both will be the last one in their respective system.
//! \todo: We should have indices for the well equations and well primary variables separately.
static constexpr int Bhp = numStaticWellEq - numWellControlEq;
static constexpr bool has_wfrac_variable = Indices::waterEnabled && Indices::oilEnabled;
static constexpr bool has_gfrac_variable = Indices::gasEnabled && Indices::numPhases > 1;
static constexpr int WFrac = has_wfrac_variable ? 1 : -1000;
static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000;
static constexpr int SFrac = !Indices::enableSolvent ? -1000 : 3;
//! \brief Evaluation for the well equations.
using EvalWell = DenseAd::DynamicEvaluation<Scalar, numStaticWellEq + Indices::numEq + 1>;
using BVectorWell = typename StandardWellEquations<Scalar,Indices::numEq>::BVectorWell;
//! \brief Constructor initializes reference to well interface.
StandardWellPrimaryVariables(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well)
: well_(well)
{}
//! \brief Initialize evaluations from values.
void init();
//! \brief Resize values and evaluations.
void resize(const int numWellEq);
//! \brief Returns number of well equations.
int numWellEq() const { return numWellEq_; }
//! \brief Copy values from well state.
void update(const WellState& well_state, DeferredLogger& deferred_logger);
//! \brief Copy polymer molecular weigt values from well state.
void updatePolyMW(const WellState& well_state);
//! \brief Update values from newton update vector.
void updateNewton(const BVectorWell& dwells,
const double dFLimit,
const double dBHPLimit);
//! \brief Update polymer molecular weight values from newton update vector.
void updateNewtonPolyMW(const BVectorWell& dwells);
//! \brief Check that all values are finite.
void checkFinite(DeferredLogger& deferred_logger) const;
//! \brief Copy values to well state.
void copyToWellState(WellState& well_state, DeferredLogger& deferred_logger) const;
//! \brief Copy polymer molecular weight values to well state.
void copyToWellStatePolyMW(WellState& well_state) const;
//! \brief Returns scaled volume fraction for a component.
EvalWell volumeFractionScaled(const int compIdx) const;
//! \brief Returns surface volume fraction for a component.
EvalWell surfaceVolumeFraction(const int compIdx) const;
//! \brief Returns scaled rate for a component.
EvalWell getQs(const int compIdx) const;
//! \brief Returns a const ref to an evaluation.
Scalar value(const int idx) const
{ return value_[idx]; }
//! \brief Returns a const ref to an evaluation.
const EvalWell& eval(const int idx) const
{ return evaluation_[idx]; }
private:
//! \brief Calculate a relaxation factor for producers.
//! \details To avoid overshoot of the fractions which might result in negative rates.
double relaxationFactorFractionsProducer(const BVectorWell& dwells) const;
//! \brief Returns volume fraction for a component.
EvalWell volumeFraction(const unsigned compIdx) const;
//! \brief Handle non-reasonable fractions due to numerical overshoot.
void processFractions();
//! \brief The values for the primary variables.
//! \details Based on different solution strategies, the wells can have different primary variables.
std::vector<Scalar> value_;
//! \brief The Evaluation for the well primary variables.
//! \details Contain derivatives and are used in AD calculation
std::vector<EvalWell> evaluation_;
const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well_; //!< Reference to well interface
//! \brief Total number of the well equations and primary variables.
//! \details There might be extra equations be used, numWellEq will be updated during the initialization
int numWellEq_ = numStaticWellEq;
};
}
#endif // OPM_STANDARDWELL_PRIMARY_VARIABLES_HEADER_INCLUDED

141
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -77,7 +77,7 @@ namespace Opm
void StandardWell<TypeTag>::
initPrimaryVariablesEvaluation() const
{
this->StdWellEval::initPrimaryVariablesEvaluation();
this->primary_variables_.init();
}
@ -105,7 +105,7 @@ namespace Opm
const EvalWell rv = this->extendEval(fs.Rv());
const EvalWell rvw = this->extendEval(fs.Rvw());
std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->numWellEq_ + Indices::numEq, 0.0});
std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
@ -125,18 +125,18 @@ namespace Opm
}
}
EvalWell skin_pressure = EvalWell{this->numWellEq_ + Indices::numEq, 0.0};
EvalWell skin_pressure = EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
if (has_polymermw) {
if (this->isInjector()) {
const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
skin_pressure = this->primary_variables_evaluation_[pskin_index];
skin_pressure = this->primary_variables_.eval(pskin_index);
}
}
// surface volume fraction of fluids within wellbore
std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->numWellEq_ + Indices::numEq});
std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->primary_variables_.numWellEq() + Indices::numEq});
for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
cmix_s[componentIdx] = this->wellSurfaceVolumeFraction(componentIdx);
cmix_s[componentIdx] = this->primary_variables_.surfaceVolumeFraction(componentIdx);
}
computePerfRate(mob,
@ -199,7 +199,7 @@ namespace Opm
if (has_polymermw) {
if (this->isInjector()) {
const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
skin_pressure = getValue(this->primary_variables_evaluation_[pskin_index]);
skin_pressure = getValue(this->primary_variables_.eval(pskin_index));
}
}
@ -210,7 +210,7 @@ namespace Opm
// surface volume fraction of fluids within wellbore
std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
cmix_s[componentIdx] = getValue(this->wellSurfaceVolumeFraction(componentIdx));
cmix_s[componentIdx] = getValue(this->primary_variables_.surfaceVolumeFraction(componentIdx));
}
computePerfRate(mob,
@ -465,9 +465,9 @@ namespace Opm
auto& perf_rates = perf_data.phase_rates;
for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
// Calculate perforation quantities.
std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
EvalWell water_flux_s{this->numWellEq_ + Indices::numEq, 0.0};
EvalWell cq_s_zfrac_effective{this->numWellEq_ + Indices::numEq, 0.0};
std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
// Equation assembly for this perforation.
@ -487,7 +487,7 @@ namespace Opm
assemblePerforationEq(cq_s_effective,
componentIdx,
cell_idx,
this->numWellEq_,
this->primary_variables_.numWellEq(),
this->linSys_);
// Store the perforation phase flux for later usage.
@ -503,7 +503,7 @@ namespace Opm
StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
assembleZFracEq(cq_s_zfrac_effective,
cell_idx,
this->numWellEq_,
this->primary_variables_.numWellEq(),
this->linSys_);
}
}
@ -526,28 +526,29 @@ namespace Opm
for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
// TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume
// since all the rates are under surface condition
EvalWell resWell_loc(this->numWellEq_ + Indices::numEq, 0.0);
EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
if (FluidSystem::numActivePhases() > 1) {
assert(dt > 0);
resWell_loc += (this->wellSurfaceVolumeFraction(componentIdx) - this->F0_[componentIdx]) * volume / dt;
resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
this->F0_[componentIdx]) * volume / dt;
}
resWell_loc -= this->getQs(componentIdx) * this->well_efficiency_factor_;
resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
assembleSourceEq(resWell_loc,
componentIdx,
this->numWellEq_,
this->primary_variables_.numWellEq(),
this->linSys_);
}
const auto& summaryState = ebosSimulator.vanguard().summaryState();
const Schedule& schedule = ebosSimulator.vanguard().schedule();
std::function<EvalWell(int)> gQ = [this](int a) { return this->getQs(a); };
std::function<EvalWell(int)> gQ = [this](int a) { return this->primary_variables_.getQs(a); };
StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
assembleControlEq(well_state, group_state,
schedule, summaryState,
this->numWellEq_,
this->getWQTotal(),
this->getBhp(),
this->primary_variables_.numWellEq(),
this->primary_variables_.eval(WQTotal),
this->primary_variables_.eval(Bhp),
gQ,
this->getRho(),
Bhp,
@ -579,10 +580,10 @@ namespace Opm
DeferredLogger& deferred_logger) const
{
const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
const EvalWell& bhp = this->getBhp();
const EvalWell& bhp = this->primary_variables_.eval(Bhp);
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
getMobilityEval(ebosSimulator, perf, mob, deferred_logger);
double perf_dis_gas_rate = 0.;
@ -627,7 +628,7 @@ namespace Opm
}
// convert to reservoir conditions
EvalWell cq_r_thermal(this->numWellEq_ + Indices::numEq, 0.);
EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
if ( !both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx ) {
@ -917,7 +918,7 @@ namespace Opm
// for the cases related to polymer molecular weight, we assume fully mixing
// as a result, the polymer and water share the same viscosity
if constexpr (!Base::has_polymermw) {
std::vector<EvalWell> mob_eval(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger);
for (size_t i = 0; i < mob.size(); ++i) {
mob[i] = getValue(mob_eval[i]);
@ -956,30 +957,14 @@ namespace Opm
{
const double dFLimit = this->param_.dwell_fraction_max_;
const double dBHPLimit = this->param_.dbhp_max_rel_;
this->StdWellEval::updatePrimaryVariablesNewton(dwells, dFLimit, dBHPLimit);
this->primary_variables_.updateNewton(dwells, dFLimit, dBHPLimit);
updateExtraPrimaryVariables(dwells);
for (double v : this->primary_variables_) {
if(!isfinite(v))
OPM_DEFLOG_THROW(NumericalIssue, "Infinite primary variable after newton update well: " << this->name(), deferred_logger);
}
}
template<typename TypeTag>
void
StandardWell<TypeTag>::
updateExtraPrimaryVariables(const BVectorWell& dwells) const
{
// for the water velocity and skin pressure
if constexpr (Base::has_polymermw) {
this->updatePrimaryVariablesPolyMW(dwells);
this->primary_variables_.updateNewtonPolyMW(dwells);
}
this->primary_variables_.checkFinite(deferred_logger);
}
@ -995,7 +980,7 @@ namespace Opm
// other primary variables related to polymer injectivity study
if constexpr (Base::has_polymermw) {
this->updateWellStateFromPrimaryVariablesPolyMW(well_state);
this->primary_variables_.copyToWellStatePolyMW(well_state);
}
}
@ -1218,7 +1203,7 @@ namespace Opm
const auto& fs = intQuants.fluidState();
const double pressure = this->getPerfCellPressure(fs).value();
const double bhp = this->getBhp().value();
const double bhp = this->primary_variables_.eval(Bhp).value();
// Pressure drawdown (also used to determine direction of flow)
const double well_pressure = bhp + this->perf_pressure_diffs_[perf];
@ -1535,7 +1520,7 @@ namespace Opm
return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
};
std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
getMobilityEval(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
const auto& fs = fluidState(subsetPerfID);
@ -1684,7 +1669,7 @@ namespace Opm
// We assemble the well equations, then we check the convergence,
// which is why we do not put the assembleWellEq here.
BVectorWell dx_well(1);
dx_well[0].resize(this->numWellEq_);
dx_well[0].resize(this->primary_variables_.numWellEq());
this->linSys_.solve( dx_well);
updateWellState(dx_well, well_state, deferred_logger);
@ -1751,7 +1736,7 @@ namespace Opm
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
BVectorWell xw(1);
xw[0].resize(this->numWellEq_);
xw[0].resize(this->primary_variables_.numWellEq());
this->linSys_.recoverSolutionWell(x, xw);
updateWellState(xw, well_state, deferred_logger);
@ -2029,26 +2014,16 @@ namespace Opm
StandardWell<TypeTag>::
updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_logger) const
{
this->StdWellEval::updatePrimaryVariables(well_state, deferred_logger);
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
this->primary_variables_.update(well_state, deferred_logger);
// other primary variables related to polymer injection
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
const auto& ws = well_state.well(this->index_of_well_);
const auto& perf_data = ws.perf_data;
const auto& water_velocity = perf_data.water_velocity;
const auto& skin_pressure = perf_data.skin_pressure;
for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
this->primary_variables_[Bhp + 1 + perf] = water_velocity[perf];
this->primary_variables_[Bhp + 1 + this->number_of_perforations_ + perf] = skin_pressure[perf];
}
}
}
for (double v : this->primary_variables_) {
if(!isfinite(v))
OPM_DEFLOG_THROW(NumericalIssue, "Infinite primary variable after update from wellState well: " << this->name(), deferred_logger);
this->primary_variables_.updatePolyMW(well_state);
}
this->primary_variables_.checkFinite(deferred_logger);
}
@ -2095,9 +2070,9 @@ namespace Opm
// compute the well water velocity with out shear effects.
// TODO: do we need to turn on crossflow here?
const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
const EvalWell& bhp = this->getBhp();
const EvalWell& bhp = this->primary_variables_.eval(Bhp);
std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
double perf_vap_wat_rate = 0.;
@ -2171,9 +2146,9 @@ namespace Opm
OPM_DEFLOG_THROW(std::runtime_error, "Unused SKPRWAT table id used for well " << name(), deferred_logger);
}
const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_water(this->numWellEq_ + Indices::numEq, 0.0);
EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
pskin_water = water_table_func.eval(throughput_eval, water_velocity);
return pskin_water;
} else {
@ -2206,9 +2181,9 @@ namespace Opm
}
const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
const double reference_concentration = skprpolytable.refConcentration;
const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_poly(this->numWellEq_ + Indices::numEq, 0.0);
EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
if (poly_inj_conc == reference_concentration) {
return sign * pskin_poly;
@ -2237,8 +2212,8 @@ namespace Opm
if constexpr (Base::has_polymermw) {
const int table_id = this->well_ecl_.getPolymerProperties().m_plymwinjtable;
const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
EvalWell molecular_weight(this->numWellEq_ + Indices::numEq, 0.);
const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
if (this->wpolymer() == 0.) { // not injecting polymer
return molecular_weight;
}
@ -2264,7 +2239,7 @@ namespace Opm
auto& ws = well_state.well(this->index_of_well_);
auto& perf_water_throughput = ws.perf_data.water_throughput;
for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
const double perf_water_vel = this->primary_variables_[Bhp + 1 + perf];
const double perf_water_vel = this->primary_variables_.value(Bhp + 1 + perf);
// we do not consider the formation damage due to water flowing from reservoir into wellbore
if (perf_water_vel > 0.) {
perf_water_throughput[perf] += perf_water_vel * dt;
@ -2295,7 +2270,7 @@ namespace Opm
// water rate is update to use the form from water velocity, since water velocity is
// a primary variable now
cq_s[water_comp_idx] = area * this->primary_variables_evaluation_[wat_vel_index] * b_w;
cq_s[water_comp_idx] = area * this->primary_variables_.eval(wat_vel_index) * b_w;
}
@ -2320,7 +2295,7 @@ namespace Opm
const int wat_vel_index = Bhp + 1 + perf;
// equation for the water velocity
const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
const EvalWell eq_wat_vel = this->primary_variables_.eval(wat_vel_index) - water_velocity;
const auto& ws = well_state.well(this->index_of_well_);
const auto& perf_data = ws.perf_data;
@ -2328,12 +2303,12 @@ namespace Opm
const double throughput = perf_water_throughput[perf];
const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
EvalWell poly_conc(this->numWellEq_ + Indices::numEq, 0.0);
EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
poly_conc.setValue(this->wpolymer());
// equation for the skin pressure
const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index]
- pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
const EvalWell eq_pskin = this->primary_variables_.eval(pskin_index)
- pskin(throughput, this->primary_variables_.eval(wat_vel_index), poly_conc, deferred_logger);
StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
assembleInjectivityEq(eq_pskin,
@ -2341,7 +2316,7 @@ namespace Opm
pskin_index,
wat_vel_index,
cell_idx,
this->numWellEq_,
this->primary_variables_.numWellEq(),
this->linSys_);
}
@ -2382,7 +2357,7 @@ namespace Opm
EvalWell cq_s_polymw = cq_s_poly;
if (this->isInjector()) {
const int wat_vel_index = Bhp + 1 + perf;
const EvalWell water_velocity = this->primary_variables_evaluation_[wat_vel_index];
const EvalWell water_velocity = this->primary_variables_.eval(wat_vel_index);
if (water_velocity > 0.) { // injecting
const auto& ws = well_state.well(this->index_of_well_);
const auto& perf_water_throughput = ws.perf_data.water_throughput;
@ -2576,7 +2551,7 @@ namespace Opm
{
// Calculate the rates that follow from the current primary variables.
std::vector<double> well_q_s(this->num_components_, 0.);
const EvalWell& bhp = this->getBhp();
const EvalWell& bhp = this->primary_variables_.eval(Bhp);
const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
const int cell_idx = this->well_cells_[perf];
@ -2674,7 +2649,7 @@ namespace Opm
deferred_logger);
}
const auto zero = EvalWell { this->numWellEq_ + Indices::numEq, 0.0 };
const auto zero = EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
connII[phase_pos] = connIICalc(mt.value() * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
}