OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-25 16:51:00 -06:00
Code Issues Packages Projects Releases Wiki Activity
f4e0a996b9
opm-simulators/opm/simulators/wells/StandardWellPrimaryVariables.cpp
Kai Bao 2c5a4398c9 make sure zero production rates are obtained for the following wells 
1. stopped production wells
2. production wells under zero rate control

We guarantee the objective through enforce zero values for the WQTotal
primary variable during the initialization and update process during the
Newton solution.
2023-03-27 13:24:08 +02:00

746 lines
29 KiB
C++
Raw Blame History

/*
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 <opm/common/Exceptions.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/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 an 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,
const bool zero_rate_target,
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 (zero_rate_target) {
value_[WQTotal] = 0.;
}
}
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,
const bool zero_rate_target,
[[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;
if (zero_rate_target) {
value_[WQTotal] = 0.;
}
// TODO: here, we make sure it is zero for zero rated wells
// 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(NumericalProblem,
"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>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,0u,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>)
}
Reference in New Issue View Git Blame Copy Permalink