OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-11 00:41:56 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1368 from totto82/usePhaseAndCompInfoFromFluidSystem

Browse Source 
Use phase and comp info from fluid system
This commit is contained in:
Atgeirr Flø Rasmussen 2018-01-03 11:21:08 +01:00 committed by GitHub
commit bd8a31a5c0
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
10 changed files with 423 additions and 474 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

146
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -29,8 +29,6 @@
#include <opm/autodiff/BlackoilModelParameters.hpp> #include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp> #include <opm/autodiff/BlackoilWellModel.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/BlackoilDetails.hpp> #include <opm/autodiff/BlackoilDetails.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp> #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
@ -43,7 +41,6 @@
#include <opm/common/Exceptions.hpp> #include <opm/common/Exceptions.hpp>
#include <opm/common/OpmLog/OpmLog.hpp> #include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/parser/eclipse/Units/Units.hpp> #include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/core/well_controls.h>
#include <opm/simulators/timestepping/SimulatorTimer.hpp> #include <opm/simulators/timestepping/SimulatorTimer.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp> #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@ -144,7 +141,6 @@ namespace Opm {
, grid_(ebosSimulator_.gridManager().grid()) , grid_(ebosSimulator_.gridManager().grid())
, istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) ) , istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) )
, phaseUsage_(phaseUsageFromDeck(eclState())) , phaseUsage_(phaseUsageFromDeck(eclState()))
, active_(detail::activePhases(phaseUsage_))
, has_disgas_(FluidSystem::enableDissolvedGas()) , has_disgas_(FluidSystem::enableDissolvedGas())
, has_vapoil_(FluidSystem::enableVaporizedOil()) , has_vapoil_(FluidSystem::enableVaporizedOil())
, has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
@ -153,7 +149,7 @@ namespace Opm {
, well_model_ (well_model) , well_model_ (well_model)
, terminal_output_ (terminal_output) , terminal_output_ (terminal_output)
, current_relaxation_(1.0) , current_relaxation_(1.0)
, dx_old_(AutoDiffGrid::numCells(grid_)) , dx_old_(UgGridHelpers::numCells(grid_))
{ {
// compute global sum of number of cells // compute global sum of number of cells
global_nc_ = detail::countGlobalCells(grid_); global_nc_ = detail::countGlobalCells(grid_);
@ -276,7 +272,7 @@ namespace Opm {
//residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ; //residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ;
// Compute the nonlinear update. // Compute the nonlinear update.
const int nc = AutoDiffGrid::numCells(grid_); const int nc = UgGridHelpers::numCells(grid_);
BVector x(nc); BVector x(nc);
try { try {
@ -596,8 +592,6 @@ namespace Opm {
/// \param[in, out] well_state well state variables /// \param[in, out] well_state well state variables
void updateState(const BVector& dx) void updateState(const BVector& dx)
{ {
using namespace Opm::AutoDiffGrid;
const auto& ebosProblem = ebosSimulator_.problem(); const auto& ebosProblem = ebosSimulator_.problem();
unsigned numSwitched = 0; unsigned numSwitched = 0;
@ -626,8 +620,8 @@ namespace Opm {
p = std::max(p, 0.0); p = std::max(p, 0.0);
// Saturation updates. // Saturation updates.
const double dsw = active_[Water] ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0; const double dsw = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0;
const double dxvar = active_[Gas] ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0; const double dxvar = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0;
double dso = 0.0; double dso = 0.0;
double dsg = 0.0; double dsg = 0.0;
@ -669,12 +663,12 @@ namespace Opm {
satScaleFactor = dsMax()/maxVal; satScaleFactor = dsMax()/maxVal;
} }
if (active_[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
double& sw = priVars[Indices::waterSaturationIdx]; double& sw = priVars[Indices::waterSaturationIdx];
sw -= satScaleFactor * dsw; sw -= satScaleFactor * dsw;
} }
if (active_[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
double& sg = priVars[Indices::compositionSwitchIdx]; double& sg = priVars[Indices::compositionSwitchIdx];
sg -= satScaleFactor * dsg; sg -= satScaleFactor * dsg;
@ -693,7 +687,7 @@ namespace Opm {
} }
// Update rs and rv // Update rs and rv
if (active_[Gas] && active_[Oil] ) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ) {
unsigned pvtRegionIdx = ebosSimulator_.problem().pvtRegionIndex(cell_idx); unsigned pvtRegionIdx = ebosSimulator_.problem().pvtRegionIndex(cell_idx);
const double drmaxrel = drMaxRel(); const double drmaxrel = drMaxRel();
if (has_disgas_) { if (has_disgas_) {
@ -809,8 +803,7 @@ namespace Opm {
const double tol_mb = param_.tolerance_mb_; const double tol_mb = param_.tolerance_mb_;
const double tol_cnv = param_.tolerance_cnv_; const double tol_cnv = param_.tolerance_cnv_;
const int np = numPhases(); const int numComp = numEq;
const int numComp = numComponents();
Vector R_sum(numComp, 0.0 ); Vector R_sum(numComp, 0.0 );
Vector B_avg(numComp, 0.0 ); Vector B_avg(numComp, 0.0 );
@ -840,16 +833,19 @@ namespace Opm {
const double pvValue = ebosProblem.porosity(cell_idx) * ebosModel.dofTotalVolume( cell_idx ); const double pvValue = ebosProblem.porosity(cell_idx) * ebosModel.dofTotalVolume( cell_idx );
pvSumLocal += pvValue; pvSumLocal += pvValue;
for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx ) for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
{ {
const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx); if (!FluidSystem::phaseIsActive(phaseIdx)) {
const int ebosCompIdx = flowPhaseToEbosCompIdx(phaseIdx); continue;
}
B_avg[ phaseIdx ] += 1.0 / fs.invB(ebosPhaseIdx).value(); const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
const auto R2 = ebosResid[cell_idx][ebosCompIdx];
R_sum[ phaseIdx ] += R2; B_avg[ compIdx ] += 1.0 / fs.invB(phaseIdx).value();
maxCoeff[ phaseIdx ] = std::max( maxCoeff[ phaseIdx ], std::abs( R2 ) / pvValue ); const auto R2 = ebosResid[cell_idx][compIdx];
R_sum[ compIdx ] += R2;
maxCoeff[ compIdx ] = std::max( maxCoeff[ compIdx ], std::abs( R2 ) / pvValue );
} }
if ( has_solvent_ ) { if ( has_solvent_ ) {
@ -914,9 +910,15 @@ namespace Opm {
std::string msg = "Iter"; std::string msg = "Iter";
std::vector< std::string > key( numComp ); std::vector< std::string > key( numComp );
for (int phaseIdx = 0; phaseIdx < numPhases(); ++phaseIdx) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); if (!FluidSystem::phaseIsActive(phaseIdx)) {
key[ phaseIdx ] = std::toupper( phaseName.front() ); continue;
}
const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
const std::string& compName = FluidSystem::componentName(canonicalCompIdx);
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx);
key[ compIdx ] = std::toupper( compName.front() );
} }
if (has_solvent_) { if (has_solvent_) {
key[ solventSaturationIdx ] = "S"; key[ solventSaturationIdx ] = "S";
@ -949,16 +951,25 @@ namespace Opm {
OpmLog::debug(ss.str()); OpmLog::debug(ss.str());
} }
for (int phaseIdx = 0; phaseIdx < numPhases(); ++phaseIdx) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
if (std::isnan(mass_balance_residual[phaseIdx]) const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
|| std::isnan(CNV[phaseIdx])) { const std::string& compName = FluidSystem::componentName(canonicalCompIdx);
OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName); const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx);
if (std::isnan(mass_balance_residual[compIdx])
|| std::isnan(CNV[compIdx])) {
OPM_THROW(Opm::NumericalProblem, "NaN residual for " << compName << " equation");
} }
if (mass_balance_residual[phaseIdx] > maxResidualAllowed() if (mass_balance_residual[compIdx] > maxResidualAllowed()
|| CNV[phaseIdx] > maxResidualAllowed()) { || CNV[compIdx] > maxResidualAllowed()) {
OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName); OPM_THROW(Opm::NumericalProblem, "Too large residual for " << compName << " equation");
}
if (mass_balance_residual[compIdx] < 0
|| CNV[compIdx] < 0) {
OPM_THROW(Opm::NumericalProblem, "Negative residual for " << compName << " equation");
} }
} }
@ -972,21 +983,6 @@ namespace Opm {
return phaseUsage_.num_phases; return phaseUsage_.num_phases;
} }
int numComponents() const
{
if (numPhases() == 2) {
return 2;
}
int numComp = FluidSystem::numComponents;
if (has_solvent_)
numComp ++;
if (has_polymer_)
numComp ++;
return numComp;
}
/// Wrapper required due to not following generic API /// Wrapper required due to not following generic API
template<class T> template<class T>
std::vector<std::vector<double> > std::vector<std::vector<double> >
@ -1001,7 +997,6 @@ namespace Opm {
const auto& comm = grid_.comm(); const auto& comm = grid_.comm();
const auto& gridView = ebosSimulator().gridView(); const auto& gridView = ebosSimulator().gridView();
const int nc = gridView.size(/*codim=*/0); const int nc = gridView.size(/*codim=*/0);
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
int ntFip = *std::max_element(fipnum.begin(), fipnum.end()); int ntFip = *std::max_element(fipnum.begin(), fipnum.end());
ntFip = comm.max(ntFip); ntFip = comm.max(ntFip);
@ -1044,18 +1039,20 @@ namespace Opm {
ebosSimulator_.model().dofTotalVolume(cellIdx) ebosSimulator_.model().dofTotalVolume(cellIdx)
* intQuants.porosity().value(); * intQuants.porosity().value();
for (int phase = 0; phase < maxnp; ++phase) { for (unsigned phase = 0; phase < FluidSystem::numPhases; ++phase) {
const double b = fs.invB(flowPhaseToEbosPhaseIdx(phase)).value(); if (!FluidSystem::phaseIsActive(phase)) {
const double s = fs.saturation(flowPhaseToEbosPhaseIdx(phase)).value(); continue;
fip_.fip[phase][cellIdx] = b * s * pv;
if (active_[ phase ]) {
regionValues[regionIdx][phase] += fip_.fip[phase][cellIdx];
} }
const double b = fs.invB(phase).value();
const double s = fs.saturation(phase).value();
const unsigned flowCanonicalPhaseIdx = ebosPhaseToFlowCanonicalPhaseIdx(phase);
fip_.fip[flowCanonicalPhaseIdx][cellIdx] = b * s * pv;
regionValues[regionIdx][flowCanonicalPhaseIdx] += fip_.fip[flowCanonicalPhaseIdx][cellIdx];
} }
if (active_[ Oil ] && active_[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// Account for gas dissolved in oil and vaporized oil // Account for gas dissolved in oil and vaporized oil
fip_.fip[FIPDataType::FIP_DISSOLVED_GAS][cellIdx] = fs.Rs().value() * fip_.fip[FIPDataType::FIP_LIQUID][cellIdx]; fip_.fip[FIPDataType::FIP_DISSOLVED_GAS][cellIdx] = fs.Rs().value() * fip_.fip[FIPDataType::FIP_LIQUID][cellIdx];
fip_.fip[FIPDataType::FIP_VAPORIZED_OIL][cellIdx] = fs.Rv().value() * fip_.fip[FIPDataType::FIP_VAPOUR][cellIdx]; fip_.fip[FIPDataType::FIP_VAPORIZED_OIL][cellIdx] = fs.Rv().value() * fip_.fip[FIPDataType::FIP_VAPOUR][cellIdx];
@ -1461,10 +1458,6 @@ namespace Opm {
const Grid& grid_; const Grid& grid_;
const ISTLSolverType* istlSolver_; const ISTLSolverType* istlSolver_;
const PhaseUsage phaseUsage_; const PhaseUsage phaseUsage_;
// For each canonical phase -> true if active
const std::vector<bool> active_;
// Size = # active phases. Maps active -> canonical phase indices.
const std::vector<int> cells_; // All grid cells
const bool has_disgas_; const bool has_disgas_;
const bool has_vapoil_; const bool has_vapoil_;
const bool has_solvent_; const bool has_solvent_;
@ -1494,29 +1487,14 @@ namespace Opm {
const BlackoilWellModel<TypeTag>& const BlackoilWellModel<TypeTag>&
wellModel() const { return well_model_; } wellModel() const { return well_model_; }
int flowPhaseToEbosCompIdx( const int phaseIdx ) const int ebosPhaseToFlowCanonicalPhaseIdx( const int phaseIdx ) const
{ {
const auto& pu = phaseUsage_; if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::waterPhaseIdx == phaseIdx)
if (active_[Water] && pu.phase_pos[Water] == phaseIdx) return Water;
return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::oilPhaseIdx == phaseIdx)
if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) return Oil;
return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::gasPhaseIdx == phaseIdx)
if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) return Gas;
return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// for other phases return the index
return phaseIdx;
}
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
{
const auto& pu = phaseUsage_;
if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
return FluidSystem::waterPhaseIdx;
if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
return FluidSystem::oilPhaseIdx;
if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
return FluidSystem::gasPhaseIdx;
assert(phaseIdx < 3); assert(phaseIdx < 3);
// for other phases return the index // for other phases return the index

14
opm/autodiff/BlackoilWellModel.hpp
View File

@ -69,19 +69,15 @@ namespace Opm {
typedef WellStateFullyImplicitBlackoil WellState; typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilModelParameters ModelParameters; typedef BlackoilModelParameters ModelParameters;
static const int Water = WellInterface<TypeTag>::Water;
static const int Oil = WellInterface<TypeTag>::Oil;
static const int Gas = WellInterface<TypeTag>::Gas;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
static const int numEq = BlackoilIndices::numEq; static const int numEq = Indices::numEq;
static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; static const int solventSaturationIdx = Indices::solventSaturationIdx;
// TODO: where we should put these types, WellInterface or Well Model? // TODO: where we should put these types, WellInterface or Well Model?
// or there is some other strategy, like TypeTag // or there is some other strategy, like TypeTag
@ -176,7 +172,6 @@ namespace Opm {
bool has_polymer_; bool has_polymer_;
std::vector<int> pvt_region_idx_; std::vector<int> pvt_region_idx_;
PhaseUsage phase_usage_; PhaseUsage phase_usage_;
std::vector<bool> active_;
size_t global_nc_; size_t global_nc_;
// the number of the cells in the local grid // the number of the cells in the local grid
size_t number_of_cells_; size_t number_of_cells_;
@ -241,9 +236,6 @@ namespace Opm {
int numPhases() const; int numPhases() const;
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
void resetWellControlFromState() const; void resetWellControlFromState() const;
void assembleWellEq(const double dt, void assembleWellEq(const double dt,

41
opm/autodiff/BlackoilWellModel_impl.hpp
View File

@ -18,11 +18,6 @@ namespace Opm {
const auto& eclState = ebosSimulator_.gridManager().eclState(); const auto& eclState = ebosSimulator_.gridManager().eclState();
phase_usage_ = phaseUsageFromDeck(eclState); phase_usage_ = phaseUsageFromDeck(eclState);
active_.resize(phase_usage_.MaxNumPhases, false);
for (int p = 0; p < phase_usage_.MaxNumPhases; ++p) {
active_[ p ] = phase_usage_.phase_used[ p ] != 0;
}
const auto& gridView = ebosSimulator_.gridView(); const auto& gridView = ebosSimulator_.gridView();
// calculate the number of elements of the compressed sequential grid. this needs // calculate the number of elements of the compressed sequential grid. this needs
@ -119,7 +114,7 @@ namespace Opm {
// TODO: to see whether we can postpone of the intialization of the well containers to // TODO: to see whether we can postpone of the intialization of the well containers to
// optimize the usage of the following several member variables // optimize the usage of the following several member variables
for (auto& well : well_container_) { for (auto& well : well_container_) {
well->init(&phase_usage_, &active_, depth_, gravity_, number_of_cells_); well->init(&phase_usage_, depth_, gravity_, number_of_cells_);
} }
// calculate the efficiency factors for each well // calculate the efficiency factors for each well
@ -378,26 +373,6 @@ namespace Opm {
template<typename TypeTag>
int
BlackoilWellModel<TypeTag>::
flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
{
const auto& pu = phase_usage_;
if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
return FluidSystem::waterPhaseIdx;
if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
return FluidSystem::oilPhaseIdx;
if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
return FluidSystem::gasPhaseIdx;
assert(phaseIdx < 3);
// for other phases return the index
return phaseIdx;
}
template<typename TypeTag> template<typename TypeTag>
void void
BlackoilWellModel<TypeTag>:: BlackoilWellModel<TypeTag>::
@ -1004,8 +979,6 @@ namespace Opm {
BlackoilWellModel<TypeTag>:: BlackoilWellModel<TypeTag>::
computeAverageFormationFactor(std::vector<double>& B_avg) const computeAverageFormationFactor(std::vector<double>& B_avg) const
{ {
const int np = numPhases();
const auto& grid = ebosSimulator_.gridManager().grid(); const auto& grid = ebosSimulator_.gridManager().grid();
const auto& gridView = grid.leafGridView(); const auto& gridView = grid.leafGridView();
ElementContext elemCtx(ebosSimulator_); ElementContext elemCtx(ebosSimulator_);
@ -1020,12 +993,16 @@ namespace Opm {
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState(); const auto& fs = intQuants.fluidState();
for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx ) for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
{ {
auto& B = B_avg[ phaseIdx ]; if (!FluidSystem::phaseIsActive(phaseIdx)) {
const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx); continue;
}
B += 1 / fs.invB(ebosPhaseIdx).value(); const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
auto& B = B_avg[ compIdx ];
B += 1 / fs.invB(phaseIdx).value();
} }
if (has_solvent_) { if (has_solvent_) {
auto& B = B_avg[solventSaturationIdx]; auto& B = B_avg[solventSaturationIdx];

10
opm/autodiff/MultisegmentWell.hpp
View File

@ -40,7 +40,7 @@ namespace Opm
using typename Base::FluidSystem; using typename Base::FluidSystem;
using typename Base::ModelParameters; using typename Base::ModelParameters;
using typename Base::MaterialLaw; using typename Base::MaterialLaw;
using typename Base::BlackoilIndices; using typename Base::Indices;
using typename Base::RateConverterType; using typename Base::RateConverterType;
@ -61,7 +61,7 @@ namespace Opm
// TODO: the following system looks not rather flexible. Looking into all kinds of possibilities // TODO: the following system looks not rather flexible. Looking into all kinds of possibilities
// TODO: gas is always there? how about oil water case? // TODO: gas is always there? how about oil water case?
// Is it gas oil two phase case? // Is it gas oil two phase case?
static const bool gasoil = numEq == 2 && (BlackoilIndices::compositionSwitchIdx >= 0); static const bool gasoil = numEq == 2 && (Indices::compositionSwitchIdx >= 0);
static const int GTotal = 0; static const int GTotal = 0;
static const int WFrac = gasoil? -1000: 1; static const int WFrac = gasoil? -1000: 1;
static const int GFrac = gasoil? 1 : 2; static const int GFrac = gasoil? 1 : 2;
@ -106,7 +106,6 @@ namespace Opm
const int num_components); const int num_components);
virtual void init(const PhaseUsage* phase_usage_arg, virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
const std::vector<double>& depth_arg, const std::vector<double>& depth_arg,
const double gravity_arg, const double gravity_arg,
const int num_cells); const int num_cells);
@ -187,11 +186,10 @@ namespace Opm
using Base::num_components_; using Base::num_components_;
// protected functions from the Base class // protected functions from the Base class
using Base::active;
using Base::phaseUsage; using Base::phaseUsage;
using Base::name; using Base::name;
using Base::flowPhaseToEbosPhaseIdx;
using Base::flowPhaseToEbosCompIdx; using Base::flowPhaseToEbosCompIdx;
using Base::ebosCompIdxToFlowCompIdx;
using Base::getAllowCrossFlow; using Base::getAllowCrossFlow;
using Base::scalingFactor; using Base::scalingFactor;
@ -284,7 +282,7 @@ namespace Opm
// fraction value of the primary variables // fraction value of the primary variables
// should we just use member variables to store them instead of calculating them again and again // should we just use member variables to store them instead of calculating them again and again
EvalWell volumeFraction(const int seg, const int comp_idx) const; EvalWell volumeFraction(const int seg, const unsigned comp_idx) const;
// F_p / g_p, the basic usage of this value is because Q_p = G_t * F_p / G_p // F_p / g_p, the basic usage of this value is because Q_p = G_t * F_p / G_p
EvalWell volumeFractionScaled(const int seg, const int comp_idx) const; EvalWell volumeFractionScaled(const int seg, const int comp_idx) const;

269
opm/autodiff/MultisegmentWell_impl.hpp
View File

@ -105,13 +105,11 @@ namespace Opm
void void
MultisegmentWell<TypeTag>:: MultisegmentWell<TypeTag>::
init(const PhaseUsage* phase_usage_arg, init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
const std::vector<double>& depth_arg, const std::vector<double>& depth_arg,
const double gravity_arg, const double gravity_arg,
const int num_cells) const int num_cells)
{ {
Base::init(phase_usage_arg, active_arg, Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells);
depth_arg, gravity_arg, num_cells);
// TODO: for StandardWell, we need to update the perf depth here using depth_arg. // TODO: for StandardWell, we need to update the perf depth here using depth_arg.
// for MultisegmentWell, it is much more complicated. // for MultisegmentWell, it is much more complicated.
@ -274,13 +272,13 @@ namespace Opm
/* const Opm::PhaseUsage& pu = phaseUsage(); /* const Opm::PhaseUsage& pu = phaseUsage();
std::vector<double> rates(3, 0.0); std::vector<double> rates(3, 0.0);
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ]; rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ];
} }
if (active()[ Oil ]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ]; rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ];
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ]; rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ];
} */ } */
@ -433,13 +431,13 @@ namespace Opm
// TODO: the report can not handle the segment number yet. // TODO: the report can not handle the segment number yet.
if (std::isnan(flux_residual)) { if (std::isnan(flux_residual)) {
report.nan_residual_found = true; report.nan_residual_found = true;
const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx)); const auto& compName = FluidSystem::componentName(Indices::activeToCanonicalComponentIndex(eq_idx));
const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name}; const typename ConvergenceReport::ProblemWell problem_well = {name(), compName};
report.nan_residual_wells.push_back(problem_well); report.nan_residual_wells.push_back(problem_well);
} else if (flux_residual > param_.max_residual_allowed_) { } else if (flux_residual > param_.max_residual_allowed_) {
report.too_large_residual_found = true; report.too_large_residual_found = true;
const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx)); const auto& compName = FluidSystem::componentName(Indices::activeToCanonicalComponentIndex(eq_idx));
const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name}; const typename ConvergenceReport::ProblemWell problem_well = {name(), compName};
report.nan_residual_wells.push_back(problem_well); report.nan_residual_wells.push_back(problem_well);
} else { // it is a normal residual } else { // it is a normal residual
if (flux_residual > maximum_residual[eq_idx]) { if (flux_residual > maximum_residual[eq_idx]) {
@ -578,11 +576,11 @@ namespace Opm
primary_variables_[seg][GTotal] = total_seg_rate; primary_variables_[seg][GTotal] = total_seg_rate;
if (std::abs(total_seg_rate) > 0.) { if (std::abs(total_seg_rate) > 0.) {
if (active()[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water]; const int water_pos = pu.phase_pos[Water];
primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg_index + water_pos] / total_seg_rate; primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg_index + water_pos] / total_seg_rate;
} }
if (active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas]; const int gas_pos = pu.phase_pos[Gas];
primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate; primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate;
} }
@ -590,7 +588,7 @@ namespace Opm
if (well_type_ == INJECTOR) { if (well_type_ == INJECTOR) {
// only single phase injection handled // only single phase injection handled
const double* distr = well_controls_get_current_distr(well_controls_); const double* distr = well_controls_get_current_distr(well_controls_);
if (active()[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if (distr[pu.phase_pos[Water]] > 0.0) { if (distr[pu.phase_pos[Water]] > 0.0) {
primary_variables_[seg][WFrac] = 1.0; primary_variables_[seg][WFrac] = 1.0;
} else { } else {
@ -598,7 +596,7 @@ namespace Opm
} }
} }
if (active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (distr[pu.phase_pos[Gas]] > 0.0) { if (distr[pu.phase_pos[Gas]] > 0.0) {
primary_variables_[seg][GFrac] = 1.0; primary_variables_[seg][GFrac] = 1.0;
} else { } else {
@ -606,11 +604,11 @@ namespace Opm
} }
} }
} else if (well_type_ == PRODUCER) { // producers } else if (well_type_ == PRODUCER) { // producers
if (active()[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
primary_variables_[seg][WFrac] = 1.0 / number_of_phases_; primary_variables_[seg][WFrac] = 1.0 / number_of_phases_;
} }
if (active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
primary_variables_[seg][GFrac] = 1.0 / number_of_phases_; primary_variables_[seg][GFrac] = 1.0 / number_of_phases_;
} }
} }
@ -743,13 +741,13 @@ namespace Opm
const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_; const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;
for (int seg = 0; seg < numberOfSegments(); ++seg) { for (int seg = 0; seg < numberOfSegments(); ++seg) {
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int sign = dwells[seg][WFrac] > 0. ? 1 : -1; const int sign = dwells[seg][WFrac] > 0. ? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit); const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit);
primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited; primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int sign = dwells[seg][GFrac] > 0. ? 1 : -1; const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit); const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit);
primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited; primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
@ -868,25 +866,24 @@ namespace Opm
template <typename TypeTag> template <typename TypeTag>
typename MultisegmentWell<TypeTag>::EvalWell typename MultisegmentWell<TypeTag>::EvalWell
MultisegmentWell<TypeTag>:: MultisegmentWell<TypeTag>::
volumeFraction(const int seg, const int comp_idx) const volumeFraction(const int seg, const unsigned compIdx) const
{ {
const PhaseUsage& pu = phaseUsage();
if (active()[Water] && comp_idx == pu.phase_pos[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
return primary_variables_evaluation_[seg][WFrac]; return primary_variables_evaluation_[seg][WFrac];
} }
if (active()[Gas] && comp_idx == pu.phase_pos[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
return primary_variables_evaluation_[seg][GFrac]; return primary_variables_evaluation_[seg][GFrac];
} }
// Oil fraction // Oil fraction
EvalWell oil_fraction = 1.0; EvalWell oil_fraction = 1.0;
if (active()[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
oil_fraction -= primary_variables_evaluation_[seg][WFrac]; oil_fraction -= primary_variables_evaluation_[seg][WFrac];
} }
if (active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
oil_fraction -= primary_variables_evaluation_[seg][GFrac]; oil_fraction -= primary_variables_evaluation_[seg][GFrac];
} }
/* if (has_solvent) { /* if (has_solvent) {
@ -898,7 +895,6 @@ namespace Opm
template <typename TypeTag> template <typename TypeTag>
typename MultisegmentWell<TypeTag>::EvalWell typename MultisegmentWell<TypeTag>::EvalWell
MultisegmentWell<TypeTag>:: MultisegmentWell<TypeTag>::
@ -907,7 +903,7 @@ namespace Opm
// For reservoir rate control, the distr in well control is used for the // For reservoir rate control, the distr in well control is used for the
// rate conversion coefficients. For the injection well, only the distr of the injection // rate conversion coefficients. For the injection well, only the distr of the injection
// phase is not zero. // phase is not zero.
const double scale = scalingFactor(comp_idx); const double scale = scalingFactor(ebosCompIdxToFlowCompIdx(comp_idx));
if (scale > 0.) { if (scale > 0.) {
return volumeFraction(seg, comp_idx) / scale; return volumeFraction(seg, comp_idx) / scale;
} }
@ -965,9 +961,13 @@ namespace Opm
// not using number_of_phases_ because of solvent // not using number_of_phases_ because of solvent
std::vector<EvalWell> b_perfcells(num_components_, 0.0); std::vector<EvalWell> b_perfcells(num_components_, 0.0);
for (int phase = 0; phase < number_of_phases_; ++phase) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
const int phase_idx_ebos = flowPhaseToEbosPhaseIdx(phase); if (!FluidSystem::phaseIsActive(phaseIdx)) {
b_perfcells[phase] = extendEval(fs.invB(phase_idx_ebos)); continue;
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
b_perfcells[compIdx] = extendEval(fs.invB(phaseIdx));
} }
// pressure difference between the segment and the perforation // pressure difference between the segment and the perforation
@ -979,8 +979,6 @@ namespace Opm
// TODO: not 100% sure about the sign of the seg_perf_press_diff // TODO: not 100% sure about the sign of the seg_perf_press_diff
const EvalWell drawdown = (pressure_cell + cell_perf_press_diff) - (segment_pressure + perf_seg_press_diff); const EvalWell drawdown = (pressure_cell + cell_perf_press_diff) - (segment_pressure + perf_seg_press_diff);
const Opm::PhaseUsage& pu = phaseUsage();
// producing perforations // producing perforations
if ( drawdown > 0.0) { if ( drawdown > 0.0) {
// Do nothing is crossflow is not allowed // Do nothing is crossflow is not allowed
@ -994,13 +992,13 @@ namespace Opm
cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p; cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
} }
if (active()[Oil] && active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int oilpos = pu.phase_pos[Oil]; const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const int gaspos = pu.phase_pos[Gas]; const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const EvalWell cq_s_oil = cq_s[oilpos]; const EvalWell cq_s_oil = cq_s[oilCompIdx];
const EvalWell cq_s_gas = cq_s[gaspos]; const EvalWell cq_s_gas = cq_s[gasCompIdx];
cq_s[gaspos] += rs * cq_s_oil; cq_s[gasCompIdx] += rs * cq_s_oil;
cq_s[oilpos] += rv * cq_s_gas; cq_s[oilCompIdx] += rv * cq_s_gas;
} }
} else { // injecting perforations } else { // injecting perforations
// Do nothing if crossflow is not allowed // Do nothing if crossflow is not allowed
@ -1019,14 +1017,14 @@ namespace Opm
// compute volume ratio between connection and at standard conditions // compute volume ratio between connection and at standard conditions
EvalWell volume_ratio = 0.0; EvalWell volume_ratio = 0.0;
if (active()[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int watpos = pu.phase_pos[Water]; const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
volume_ratio += cmix_s[watpos] / b_perfcells[watpos]; volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
} }
if (active()[Oil] && active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int oilpos = pu.phase_pos[Oil]; const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const int gaspos = pu.phase_pos[Gas]; const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// Incorporate RS/RV factors if both oil and gas active // Incorporate RS/RV factors if both oil and gas active
// TODO: not sure we use rs rv from the perforation cells when handling injecting perforations // TODO: not sure we use rs rv from the perforation cells when handling injecting perforations
@ -1038,19 +1036,19 @@ namespace Opm
<< " with rs " << rs << " and rv " << rv); << " with rs " << rs << " and rv " << rv);
} }
const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d; const EvalWell tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
volume_ratio += tmp_oil / b_perfcells[oilpos]; volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
const EvalWell tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d; const EvalWell tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
volume_ratio += tmp_gas / b_perfcells[gaspos]; volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
} else { // not having gas and oil at the same time } else { // not having gas and oil at the same time
if (active()[Oil]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int oilpos = pu.phase_pos[Oil]; const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
volume_ratio += cmix_s[oilpos] / b_perfcells[oilpos]; volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
} }
if (active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gaspos = pu.phase_pos[Gas]; const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
volume_ratio += cmix_s[gaspos] / b_perfcells[gaspos]; volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
} }
} }
// injecting connections total volumerates at standard conditions // injecting connections total volumerates at standard conditions
@ -1110,14 +1108,17 @@ namespace Opm
pvt_region_index = fs.pvtRegionIndex(); pvt_region_index = fs.pvtRegionIndex();
} }
std::vector<double> surf_dens(number_of_phases_); std::vector<double> surf_dens(num_components_);
// Surface density. // Surface density.
// not using num_comp here is because solvent can be component for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
for (int phase = 0; phase < number_of_phases_; ++phase) { if (!FluidSystem::phaseIsActive(phaseIdx)) {
surf_dens[phase] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx(phase), pvt_region_index ); continue;
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
surf_dens[compIdx] = FluidSystem::referenceDensity( phaseIdx, pvt_region_index );
} }
const Opm::PhaseUsage& pu = phaseUsage();
for (int seg = 0; seg < numberOfSegments(); ++seg) { for (int seg = 0; seg < numberOfSegments(); ++seg) {
// the compostion of the components inside wellbore under surface condition // the compostion of the components inside wellbore under surface condition
std::vector<EvalWell> mix_s(num_components_, 0.0); std::vector<EvalWell> mix_s(num_components_, 0.0);
@ -1126,94 +1127,95 @@ namespace Opm
} }
std::vector<EvalWell> b(num_components_, 0.0); std::vector<EvalWell> b(num_components_, 0.0);
// it is the phase viscosities asked for std::vector<EvalWell> visc(num_components_, 0.0);
std::vector<EvalWell> visc(number_of_phases_, 0.0);
const EvalWell seg_pressure = getSegmentPressure(seg); const EvalWell seg_pressure = getSegmentPressure(seg);
if (pu.phase_used[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water]; const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
b[water_pos] = b[waterCompIdx] =
FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[water_pos] = visc[waterCompIdx] =
FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure); FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure);
} }
EvalWell rv(0.0); EvalWell rv(0.0);
// gas phase // gas phase
if (pu.phase_used[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gaspos = pu.phase_pos[Gas]; const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
if (pu.phase_used[Oil]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int oilpos = pu.phase_pos[Oil]; const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure); const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure);
if (mix_s[oilpos] > 0.0) { if (mix_s[oilCompIdx] > 0.0) {
if (mix_s[gaspos] > 0.0) { if (mix_s[gasCompIdx] > 0.0) {
rv = mix_s[oilpos] / mix_s[gaspos]; rv = mix_s[oilCompIdx] / mix_s[gasCompIdx];
} }
if (rv > rvmax) { if (rv > rvmax) {
rv = rvmax; rv = rvmax;
} }
b[gaspos] = b[gasCompIdx] =
FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv); FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv);
visc[gaspos] = visc[gasCompIdx] =
FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv); FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv);
} else { // no oil exists } else { // no oil exists
b[gaspos] = b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[gaspos] = visc[gasCompIdx] =
FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
} }
} else { // no Liquid phase } else { // no Liquid phase
// it is the same with zero mix_s[Oil] // it is the same with zero mix_s[Oil]
b[gaspos] = b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[gaspos] = visc[gasCompIdx] =
FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
} }
} }
EvalWell rs(0.0); EvalWell rs(0.0);
// oil phase // oil phase
if (pu.phase_used[Oil]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int oilpos = pu.phase_pos[Oil]; const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (pu.phase_used[Oil]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gaspos = pu.phase_pos[Gas]; const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure); const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure);
if (mix_s[gaspos] > 0.0) { if (mix_s[gasCompIdx] > 0.0) {
if (mix_s[oilpos] > 0.0) { if (mix_s[oilCompIdx] > 0.0) {
rs = mix_s[gaspos] / mix_s[oilpos]; rs = mix_s[gasCompIdx] / mix_s[oilCompIdx];
} }
if (rs > rsmax) { if (rs > rsmax) {
rs = rsmax; rs = rsmax;
} }
b[oilpos] = b[oilCompIdx] =
FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs); FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs);
visc[oilpos] = visc[oilCompIdx] =
FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs); FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs);
} else { // no oil exists } else { // no oil exists
b[oilpos] = b[oilCompIdx] =
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[oilpos] = visc[oilCompIdx] =
FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
} }
} else { // no Liquid phase } else { // no Liquid phase
// it is the same with zero mix_s[Oil] // it is the same with zero mix_s[Oil]
b[oilpos] = b[oilCompIdx] =
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[oilpos] = visc[oilCompIdx] =
FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
} }
} }
std::vector<EvalWell> mix(mix_s); std::vector<EvalWell> mix(mix_s);
if (pu.phase_used[Oil] && pu.phase_used[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gaspos = pu.phase_pos[Gas]; const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const int oilpos = pu.phase_pos[Oil]; const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (rs != 0.0) { // rs > 0.0? if (rs != 0.0) { // rs > 0.0?
mix[gaspos] = (mix_s[gaspos] - mix_s[oilpos] * rs) / (1. - rs * rv); mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / (1. - rs * rv);
} }
if (rv != 0.0) { // rv > 0.0? if (rv != 0.0) { // rv > 0.0?
mix[oilpos] = (mix_s[oilpos] - mix_s[gaspos] * rv) / (1. - rs * rv); mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / (1. - rs * rv);
} }
} }
@ -1224,9 +1226,9 @@ namespace Opm
segment_viscosities_[seg] = 0.; segment_viscosities_[seg] = 0.;
// calculate the average viscosity // calculate the average viscosity
for (int p = 0; p < number_of_phases_; ++p) { for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
const EvalWell phase_fraction = mix[p] / b[p] / volrat; const EvalWell comp_fraction = mix[comp_idx] / b[comp_idx] / volrat;
segment_viscosities_[seg] += visc[p] * phase_fraction; segment_viscosities_[seg] += visc[comp_idx] * comp_fraction;
} }
EvalWell density(0.0); EvalWell density(0.0);
@ -1237,9 +1239,9 @@ namespace Opm
// calculate the mass rates // calculate the mass rates
segment_mass_rates_[seg] = 0.; segment_mass_rates_[seg] = 0.;
for (int phase = 0; phase < number_of_phases_; ++phase) { for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
const EvalWell rate = getSegmentRate(seg, phase); const EvalWell rate = getSegmentRate(seg, comp_idx);
segment_mass_rates_[seg] += rate * surf_dens[phase]; segment_mass_rates_[seg] += rate * surf_dens[comp_idx];
} }
} }
} }
@ -1293,7 +1295,6 @@ namespace Opm
std::vector<EvalWell>& mob) const std::vector<EvalWell>& mob) const
{ {
// TODO: most of this function, if not the whole function, can be moved to the base class // TODO: most of this function, if not the whole function, can be moved to the base class
const int np = number_of_phases_;
const int cell_idx = well_cells_[perf]; const int cell_idx = well_cells_[perf];
assert (int(mob.size()) == num_components_); assert (int(mob.size()) == num_components_);
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
@ -1305,9 +1306,13 @@ namespace Opm
const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell
for (int phase = 0; phase < np; ++phase) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); if (!FluidSystem::phaseIsActive(phaseIdx)) {
mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); continue;
}
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
} }
// if (has_solvent) { // if (has_solvent) {
// mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); // mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
@ -1322,9 +1327,13 @@ namespace Opm
materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
// compute the mobility // compute the mobility
for (int phase = 0; phase < np; ++phase) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); if (!FluidSystem::phaseIsActive(phaseIdx)) {
mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx)); continue;
}
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
} }
} }
@ -1371,7 +1380,7 @@ namespace Opm
if (number_phases_under_control == 1) { // single phase control if (number_phases_under_control == 1) { // single phase control
for (int phase = 0; phase < number_of_phases_; ++phase) { for (int phase = 0; phase < number_of_phases_; ++phase) {
if (distr[phase] > 0.) { // under the control of this phase if (distr[phase] > 0.) { // under the control of this phase
control_eq = getSegmentGTotal(0) * volumeFraction(0, phase) - g[phase] * target_rate; control_eq = getSegmentGTotal(0) * volumeFraction(0, flowPhaseToEbosCompIdx(phase)) - g[phase] * target_rate;
break; break;
} }
} }
@ -1380,7 +1389,7 @@ namespace Opm
const EvalWell G_total = getSegmentGTotal(0); const EvalWell G_total = getSegmentGTotal(0);
for (int phase = 0; phase < number_of_phases_; ++phase) { for (int phase = 0; phase < number_of_phases_; ++phase) {
if (distr[phase] > 0.) { if (distr[phase] > 0.) {
rate_for_control += G_total * volumeFractionScaled(0, phase); rate_for_control += G_total * volumeFractionScaled(0, flowPhaseToEbosCompIdx(phase));
} }
} }
// TODO: maybe the following equation can be scaled a little bit for gas phase // TODO: maybe the following equation can be scaled a little bit for gas phase
@ -1394,7 +1403,7 @@ namespace Opm
const double* distr = well_controls_get_current_distr(well_controls_); const double* distr = well_controls_get_current_distr(well_controls_);
for (int phase = 0; phase < number_of_phases_; ++phase) { for (int phase = 0; phase < number_of_phases_; ++phase) {
if (distr[phase] > 0.) { if (distr[phase] > 0.) {
rate_for_control += getSegmentGTotal(0) * volumeFraction(0, phase); rate_for_control += getSegmentGTotal(0) * volumeFraction(0, flowPhaseToEbosCompIdx(phase));
} }
} }
const double target_rate = well_controls_get_current_target(well_controls_); const double target_rate = well_controls_get_current_target(well_controls_);
@ -1546,26 +1555,26 @@ namespace Opm
std::vector<double> fractions(number_of_phases_, 0.0); std::vector<double> fractions(number_of_phases_, 0.0);
assert( active()[Oil] ); assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
const int oil_pos = pu.phase_pos[Oil]; const int oil_pos = pu.phase_pos[Oil];
fractions[oil_pos] = 1.0; fractions[oil_pos] = 1.0;
if ( active()[Water] ) { if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
const int water_pos = pu.phase_pos[Water]; const int water_pos = pu.phase_pos[Water];
fractions[water_pos] = primary_variables_[seg][WFrac]; fractions[water_pos] = primary_variables_[seg][WFrac];
fractions[oil_pos] -= fractions[water_pos]; fractions[oil_pos] -= fractions[water_pos];
} }
if ( active()[Gas] ) { if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
const int gas_pos = pu.phase_pos[Gas]; const int gas_pos = pu.phase_pos[Gas];
fractions[gas_pos] = primary_variables_[seg][GFrac]; fractions[gas_pos] = primary_variables_[seg][GFrac];
fractions[oil_pos] -= fractions[gas_pos]; fractions[oil_pos] -= fractions[gas_pos];
} }
if (active()[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water]; const int water_pos = pu.phase_pos[Water];
if (fractions[water_pos] < 0.0) { if (fractions[water_pos] < 0.0) {
if ( active()[Gas] ) { if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]); fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]);
} }
fractions[oil_pos] /= (1.0 - fractions[water_pos]); fractions[oil_pos] /= (1.0 - fractions[water_pos]);
@ -1573,10 +1582,10 @@ namespace Opm
} }
} }
if (active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas]; const int gas_pos = pu.phase_pos[Gas];
if (fractions[gas_pos] < 0.0) { if (fractions[gas_pos] < 0.0) {
if ( active()[Water] ) { if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
fractions[pu.phase_pos[Water]] /= (1.0 - fractions[gas_pos]); fractions[pu.phase_pos[Water]] /= (1.0 - fractions[gas_pos]);
} }
fractions[oil_pos] /= (1.0 - fractions[gas_pos]); fractions[oil_pos] /= (1.0 - fractions[gas_pos]);
@ -1585,20 +1594,20 @@ namespace Opm
} }
if (fractions[oil_pos] < 0.0) { if (fractions[oil_pos] < 0.0) {
if ( active()[Water] ) { if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
fractions[pu.phase_pos[Water]] /= (1.0 - fractions[oil_pos]); fractions[pu.phase_pos[Water]] /= (1.0 - fractions[oil_pos]);
} }
if ( active()[Gas] ) { if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]); fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]);
} }
fractions[oil_pos] = 0.0; fractions[oil_pos] = 0.0;
} }
if ( active()[Water] ) { if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]]; primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]];
} }
if ( active()[Gas] ) { if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]]; primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]];
} }
} }
@ -1613,20 +1622,20 @@ namespace Opm
updateWellStateFromPrimaryVariables(WellState& well_state) const updateWellStateFromPrimaryVariables(WellState& well_state) const
{ {
const PhaseUsage& pu = phaseUsage(); const PhaseUsage& pu = phaseUsage();
assert( active()[Oil] ); assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
const int oil_pos = pu.phase_pos[Oil]; const int oil_pos = pu.phase_pos[Oil];
for (int seg = 0; seg < numberOfSegments(); ++seg) { for (int seg = 0; seg < numberOfSegments(); ++seg) {
std::vector<double> fractions(number_of_phases_, 0.0); std::vector<double> fractions(number_of_phases_, 0.0);
fractions[oil_pos] = 1.0; fractions[oil_pos] = 1.0;
if ( active()[Water] ) { if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
const int water_pos = pu.phase_pos[Water]; const int water_pos = pu.phase_pos[Water];
fractions[water_pos] = primary_variables_[seg][WFrac]; fractions[water_pos] = primary_variables_[seg][WFrac];
fractions[oil_pos] -= fractions[water_pos]; fractions[oil_pos] -= fractions[water_pos];
} }
if ( active()[Gas] ) { if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
const int gas_pos = pu.phase_pos[Gas]; const int gas_pos = pu.phase_pos[Gas];
fractions[gas_pos] = primary_variables_[seg][GFrac]; fractions[gas_pos] = primary_variables_[seg][GFrac];
fractions[oil_pos] -= fractions[gas_pos]; fractions[oil_pos] -= fractions[gas_pos];
@ -1814,9 +1823,7 @@ namespace Opm
if (!only_wells) { if (!only_wells) {
// subtract sum of component fluxes in the reservoir equation. // subtract sum of component fluxes in the reservoir equation.
// need to consider the efficiency factor // need to consider the efficiency factor
// TODO: the name of the function flowPhaseToEbosCompIdx is prolematic, since the input ebosResid[cell_idx][comp_idx] -= cq_s_effective.value();
// is a component index :D
ebosResid[cell_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.value();
} }
// subtract sum of phase fluxes in the well equations. // subtract sum of phase fluxes in the well equations.
@ -1826,7 +1833,7 @@ namespace Opm
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
if (!only_wells) { if (!only_wells) {
// also need to consider the efficiency factor when manipulating the jacobians. // also need to consider the efficiency factor when manipulating the jacobians.
duneC_[seg][cell_idx][pv_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix
} }
// the index name for the D should be eq_idx / pv_idx // the index name for the D should be eq_idx / pv_idx
duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq); duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq);
@ -1835,7 +1842,7 @@ namespace Opm
for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) { for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) {
if (!only_wells) { if (!only_wells) {
// also need to consider the efficiency factor when manipulating the jacobians. // also need to consider the efficiency factor when manipulating the jacobians.
ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(comp_idx)][pv_idx] -= cq_s_effective.derivative(pv_idx); ebosJac[cell_idx][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx);
duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx); duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx);
} }
} }

2
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -29,8 +29,6 @@
#include <opm/autodiff/BlackoilModelParameters.hpp> #include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp> #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp> #include <opm/autodiff/BlackoilWellModel.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/SimFIBODetails.hpp>
#include <opm/autodiff/moduleVersion.hpp> #include <opm/autodiff/moduleVersion.hpp>
#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp> #include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
#include <opm/core/utility/initHydroCarbonState.hpp> #include <opm/core/utility/initHydroCarbonState.hpp>

10
opm/autodiff/StandardWell.hpp
View File

@ -46,7 +46,7 @@ namespace Opm
using typename Base::FluidSystem; using typename Base::FluidSystem;
using typename Base::MaterialLaw; using typename Base::MaterialLaw;
using typename Base::ModelParameters; using typename Base::ModelParameters;
using typename Base::BlackoilIndices; using typename Base::Indices;
using typename Base::PolymerModule; using typename Base::PolymerModule;
using typename Base::RateConverterType; using typename Base::RateConverterType;
@ -56,7 +56,7 @@ namespace Opm
// there are three primary variables, the second and the third ones are F_w and F_g // there are three primary variables, the second and the third ones are F_w and F_g
// the first one can be total rate (G_t) or bhp, based on the control // the first one can be total rate (G_t) or bhp, based on the control
static const bool gasoil = numEq == 2 && (BlackoilIndices::compositionSwitchIdx >= 0); static const bool gasoil = numEq == 2 && (Indices::compositionSwitchIdx >= 0);
static const int XvarWell = 0; static const int XvarWell = 0;
static const int WFrac = gasoil? -1000: 1; static const int WFrac = gasoil? -1000: 1;
static const int GFrac = gasoil? 1: 2; static const int GFrac = gasoil? 1: 2;
@ -117,7 +117,6 @@ namespace Opm
const int num_components); const int num_components);
virtual void init(const PhaseUsage* phase_usage_arg, virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
const std::vector<double>& depth_arg, const std::vector<double>& depth_arg,
const double gravity_arg, const double gravity_arg,
const int num_cells); const int num_cells);
@ -163,9 +162,8 @@ namespace Opm
// protected functions from the Base class // protected functions from the Base class
using Base::getAllowCrossFlow; using Base::getAllowCrossFlow;
using Base::phaseUsage; using Base::phaseUsage;
using Base::active;
using Base::flowPhaseToEbosPhaseIdx;
using Base::flowPhaseToEbosCompIdx; using Base::flowPhaseToEbosCompIdx;
using Base::ebosCompIdxToFlowCompIdx;
using Base::wsolvent; using Base::wsolvent;
using Base::wpolymer; using Base::wpolymer;
using Base::wellHasTHPConstraints; using Base::wellHasTHPConstraints;
@ -233,7 +231,7 @@ namespace Opm
EvalWell wellVolumeFractionScaled(const int phase) const; EvalWell wellVolumeFractionScaled(const int phase) const;
EvalWell wellVolumeFraction(const int phase) const; EvalWell wellVolumeFraction(const unsigned compIdx) const;
EvalWell wellSurfaceVolumeFraction(const int phase) const; EvalWell wellSurfaceVolumeFraction(const int phase) const;

320
opm/autodiff/StandardWell_impl.hpp
View File

@ -49,13 +49,11 @@ namespace Opm
void void
StandardWell<TypeTag>:: StandardWell<TypeTag>::
init(const PhaseUsage* phase_usage_arg, init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
const std::vector<double>& depth_arg, const std::vector<double>& depth_arg,
const double gravity_arg, const double gravity_arg,
const int num_cells) const int num_cells)
{ {
Base::init(phase_usage_arg, active_arg, Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells);
depth_arg, gravity_arg, num_cells);
perf_depth_.resize(number_of_perforations_, 0.); perf_depth_.resize(number_of_perforations_, 0.);
for (int perf = 0; perf < number_of_perforations_; ++perf) { for (int perf = 0; perf < number_of_perforations_; ++perf) {
@ -138,13 +136,13 @@ namespace Opm
const Opm::PhaseUsage& pu = phaseUsage(); const Opm::PhaseUsage& pu = phaseUsage();
std::vector<EvalWell> rates(3, 0.0); std::vector<EvalWell> rates(3, 0.0);
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ]= getQs(pu.phase_pos[ Water]); rates[ Water ]= getQs(pu.phase_pos[ Water]);
} }
if (active()[ Oil ]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = getQs(pu.phase_pos[ Oil ]); rates[ Oil ] = getQs(pu.phase_pos[ Oil ]);
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = getQs(pu.phase_pos[ Gas ]); rates[ Gas ] = getQs(pu.phase_pos[ Gas ]);
} }
return calculateBhpFromThp(rates, control); return calculateBhpFromThp(rates, control);
@ -170,6 +168,7 @@ namespace Opm
assert(comp_idx < num_components_); assert(comp_idx < num_components_);
const auto pu = phaseUsage(); const auto pu = phaseUsage();
const int legacyCompIdx = ebosCompIdxToFlowCompIdx(comp_idx);
// TODO: the formulation for the injectors decides it only work with single phase // TODO: the formulation for the injectors decides it only work with single phase
// surface rate injection control. Improvement will be required. // surface rate injection control. Improvement will be required.
@ -180,10 +179,10 @@ namespace Opm
double comp_frac = 0.0; double comp_frac = 0.0;
if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent
comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent(); comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent();
} else if (comp_idx == pu.phase_pos[ Gas ]) { } else if (legacyCompIdx == pu.phase_pos[ Gas ]) {
comp_frac = comp_frac_[comp_idx] * (1.0 - wsolvent()); comp_frac = comp_frac_[legacyCompIdx] * (1.0 - wsolvent());
} else { } else {
comp_frac = comp_frac_[comp_idx]; comp_frac = comp_frac_[legacyCompIdx];
} }
if (comp_frac == 0.0) { if (comp_frac == 0.0) {
return qs; //zero return qs; //zero
@ -197,7 +196,7 @@ namespace Opm
return qs; return qs;
} }
const double comp_frac = comp_frac_[comp_idx]; const double comp_frac = comp_frac_[legacyCompIdx];
if (comp_frac == 0.0) { if (comp_frac == 0.0) {
return qs; return qs;
} }
@ -242,15 +241,16 @@ namespace Opm
assert(phase_under_control >= 0); assert(phase_under_control >= 0);
EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(phase_under_control); const int compIdx_under_control = flowPhaseToEbosCompIdx(phase_under_control);
EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(compIdx_under_control);
if (has_solvent && phase_under_control == Gas) { if (has_solvent && phase_under_control == Gas) {
// for GRAT controlled wells solvent is included in the target // for GRAT controlled wells solvent is included in the target
wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(contiSolventEqIdx); wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(contiSolventEqIdx);
} }
if (comp_idx == phase_under_control) { if (comp_idx == compIdx_under_control) {
if (has_solvent && phase_under_control == Gas) { if (has_solvent && compIdx_under_control == FluidSystem::gasCompIdx) {
qs.setValue(target_rate * wellVolumeFractionScaled(Gas).value() / wellVolumeFractionScaledPhaseUnderControl.value() ); qs.setValue(target_rate * wellVolumeFractionScaled(compIdx_under_control).value() / wellVolumeFractionScaledPhaseUnderControl.value() );
return qs; return qs;
} }
qs.setValue(target_rate); qs.setValue(target_rate);
@ -270,8 +270,9 @@ namespace Opm
if (num_phases_under_rate_control == 2) { if (num_phases_under_rate_control == 2) {
EvalWell combined_volume_fraction = 0.; EvalWell combined_volume_fraction = 0.;
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
const unsigned compIdxTmp = flowPhaseToEbosCompIdx(p);
if (distr[p] == 1.0) { if (distr[p] == 1.0) {
combined_volume_fraction += wellVolumeFractionScaled(p); combined_volume_fraction += wellVolumeFractionScaled(compIdxTmp);
} }
} }
return (target_rate * wellVolumeFractionScaled(comp_idx) / combined_volume_fraction); return (target_rate * wellVolumeFractionScaled(comp_idx) / combined_volume_fraction);
@ -303,7 +304,8 @@ namespace Opm
wellVolumeFractionScaled(const int compIdx) const wellVolumeFractionScaled(const int compIdx) const
{ {
const double scal = scalingFactor(compIdx); const int legacyCompIdx = ebosCompIdxToFlowCompIdx(compIdx);
const double scal = scalingFactor(legacyCompIdx);
if (scal > 0) if (scal > 0)
return wellVolumeFraction(compIdx) / scal; return wellVolumeFraction(compIdx) / scal;
@ -318,14 +320,13 @@ namespace Opm
template<typename TypeTag> template<typename TypeTag>
typename StandardWell<TypeTag>::EvalWell typename StandardWell<TypeTag>::EvalWell
StandardWell<TypeTag>:: StandardWell<TypeTag>::
wellVolumeFraction(const int compIdx) const wellVolumeFraction(const unsigned compIdx) const
{ {
const auto& pu = phaseUsage(); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
if (active()[Water] && compIdx == pu.phase_pos[Water]) {
return primary_variables_evaluation_[WFrac]; return primary_variables_evaluation_[WFrac];
} }
if (active()[Gas] && compIdx == pu.phase_pos[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
return primary_variables_evaluation_[GFrac]; return primary_variables_evaluation_[GFrac];
} }
@ -335,11 +336,11 @@ namespace Opm
// Oil fraction // Oil fraction
EvalWell well_fraction = 1.0; EvalWell well_fraction = 1.0;
if (active()[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
well_fraction -= primary_variables_evaluation_[WFrac]; well_fraction -= primary_variables_evaluation_[WFrac];
} }
if (active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
well_fraction -= primary_variables_evaluation_[GFrac]; well_fraction -= primary_variables_evaluation_[GFrac];
} }
if (has_solvent) { if (has_solvent) {
@ -396,21 +397,22 @@ namespace Opm
const double Tw, const EvalWell& bhp, const double& cdp, const double Tw, const EvalWell& bhp, const double& cdp,
const bool& allow_cf, std::vector<EvalWell>& cq_s) const const bool& allow_cf, std::vector<EvalWell>& cq_s) const
{ {
const Opm::PhaseUsage& pu = phaseUsage();
const int np = number_of_phases_;
std::vector<EvalWell> cmix_s(num_components_,0.0); std::vector<EvalWell> cmix_s(num_components_,0.0);
for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx); cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx);
} }
const auto& fs = intQuants.fluidState(); const auto& fs = intQuants.fluidState();
const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
const EvalWell rs = extendEval(fs.Rs()); const EvalWell rs = extendEval(fs.Rs());
const EvalWell rv = extendEval(fs.Rv()); const EvalWell rv = extendEval(fs.Rv());
std::vector<EvalWell> b_perfcells_dense(num_components_, 0.0); std::vector<EvalWell> b_perfcells_dense(num_components_, 0.0);
for (int phase = 0; phase < np; ++phase) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); if (!FluidSystem::phaseIsActive(phaseIdx)) {
b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); continue;
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
b_perfcells_dense[compIdx] = extendEval(fs.invB(phaseIdx));
} }
if (has_solvent) { if (has_solvent) {
b_perfcells_dense[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor()); b_perfcells_dense[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
@ -433,13 +435,13 @@ namespace Opm
cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p; cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
} }
if (active()[Oil] && active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int oilpos = pu.phase_pos[Oil]; const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const int gaspos = pu.phase_pos[Gas]; const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const EvalWell cq_sOil = cq_s[oilpos]; const EvalWell cq_sOil = cq_s[oilCompIdx];
const EvalWell cq_sGas = cq_s[gaspos]; const EvalWell cq_sGas = cq_s[gasCompIdx];
cq_s[gaspos] += rs * cq_sOil; cq_s[gasCompIdx] += rs * cq_sOil;
cq_s[oilpos] += rv * cq_sGas; cq_s[oilCompIdx] += rv * cq_sGas;
} }
} else { } else {
@ -459,19 +461,18 @@ namespace Opm
// compute volume ratio between connection at standard conditions // compute volume ratio between connection at standard conditions
EvalWell volumeRatio = 0.0; EvalWell volumeRatio = 0.0;
if (active()[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int watpos = pu.phase_pos[Water]; const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
} }
if (has_solvent) { if (has_solvent) {
volumeRatio += cmix_s[contiSolventEqIdx] / b_perfcells_dense[contiSolventEqIdx]; volumeRatio += cmix_s[contiSolventEqIdx] / b_perfcells_dense[contiSolventEqIdx];
} }
if (active()[Oil] && active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int oilpos = pu.phase_pos[Oil]; const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const int gaspos = pu.phase_pos[Gas]; const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// Incorporate RS/RV factors if both oil and gas active // Incorporate RS/RV factors if both oil and gas active
const EvalWell d = 1.0 - rv * rs; const EvalWell d = 1.0 - rv * rs;
@ -480,22 +481,22 @@ namespace Opm
<< " with rs " << rs << " and rv " << rv); << " with rs " << rs << " and rv " << rv);
} }
const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d; const EvalWell tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
//std::cout << "tmp_oil " <<tmp_oil << std::endl; //std::cout << "tmp_oil " <<tmp_oil << std::endl;
volumeRatio += tmp_oil / b_perfcells_dense[oilpos]; volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
const EvalWell tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d; const EvalWell tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
//std::cout << "tmp_gas " <<tmp_gas << std::endl; //std::cout << "tmp_gas " <<tmp_gas << std::endl;
volumeRatio += tmp_gas / b_perfcells_dense[gaspos]; volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
} }
else { else {
if (active()[Oil]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int oilpos = pu.phase_pos[Oil]; const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
volumeRatio += cmix_s[oilpos] / b_perfcells_dense[oilpos]; volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
} }
if (active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gaspos = pu.phase_pos[Gas]; const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
volumeRatio += cmix_s[gaspos] / b_perfcells_dense[gaspos]; volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
} }
} }
@ -556,7 +557,7 @@ namespace Opm
if (!only_wells) { if (!only_wells) {
// subtract sum of component fluxes in the reservoir equation. // subtract sum of component fluxes in the reservoir equation.
// need to consider the efficiency factor // need to consider the efficiency factor
ebosResid[cell_idx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.value(); ebosResid[cell_idx][componentIdx] -= cq_s_effective.value();
} }
// subtract sum of phase fluxes in the well equations. // subtract sum of phase fluxes in the well equations.
@ -566,7 +567,7 @@ namespace Opm
for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) { for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
if (!only_wells) { if (!only_wells) {
// also need to consider the efficiency factor when manipulating the jacobians. // also need to consider the efficiency factor when manipulating the jacobians.
duneC_[0][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix
} }
invDuneD_[0][0][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx+numEq); invDuneD_[0][0][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx+numEq);
} }
@ -574,22 +575,23 @@ namespace Opm
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
if (!only_wells) { if (!only_wells) {
// also need to consider the efficiency factor when manipulating the jacobians. // also need to consider the efficiency factor when manipulating the jacobians.
ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][pvIdx] -= cq_s_effective.derivative(pvIdx); ebosJac[cell_idx][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx);
duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx); duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx);
} }
} }
// Store the perforation phase flux for later usage. // Store the perforation phase flux for later usage.
if (has_solvent && componentIdx == contiSolventEqIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent) if (has_solvent && componentIdx == contiSolventEqIdx) {
well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value(); well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value();
} else { } else {
well_state.perfPhaseRates()[(first_perf_ + perf) * np + componentIdx] = cq_s[componentIdx].value(); well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
} }
} }
if (has_polymer) { if (has_polymer) {
// TODO: the application of well efficiency factor has not been tested with an example yet // TODO: the application of well efficiency factor has not been tested with an example yet
EvalWell cq_s_poly = cq_s[Water] * well_efficiency_factor_; const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
EvalWell cq_s_poly = cq_s[waterCompIdx] * well_efficiency_factor_;
if (well_type_ == INJECTOR) { if (well_type_ == INJECTOR) {
cq_s_poly *= wpolymer(); cq_s_poly *= wpolymer();
} else { } else {
@ -671,7 +673,6 @@ namespace Opm
const int perf, const int perf,
std::vector<EvalWell>& mob) const std::vector<EvalWell>& mob) const
{ {
const int np = number_of_phases_;
const int cell_idx = well_cells_[perf]; const int cell_idx = well_cells_[perf];
assert (int(mob.size()) == num_components_); assert (int(mob.size()) == num_components_);
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
@ -683,9 +684,13 @@ namespace Opm
const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell
for (int phase = 0; phase < np; ++phase) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); if (!FluidSystem::phaseIsActive(phaseIdx)) {
mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); continue;
}
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
} }
if (has_solvent) { if (has_solvent) {
mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
@ -700,9 +705,13 @@ namespace Opm
materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
// compute the mobility // compute the mobility
for (int phase = 0; phase < np; ++phase) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); if (!FluidSystem::phaseIsActive(phaseIdx)) {
mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx)); continue;
}
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
} }
// this may not work if viscosity and relperms has been modified? // this may not work if viscosity and relperms has been modified?
@ -713,7 +722,7 @@ namespace Opm
// modify the water mobility if polymer is present // modify the water mobility if polymer is present
if (has_polymer) { if (has_polymer) {
if (!active()[Water]) { if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
OPM_THROW(std::runtime_error, "Water is required when polymer is active"); OPM_THROW(std::runtime_error, "Water is required when polymer is active");
} }
@ -740,13 +749,13 @@ namespace Opm
// update the second and third well variable (The flux fractions) // update the second and third well variable (The flux fractions)
std::vector<double> F(np,0.0); std::vector<double> F(np,0.0);
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int sign2 = dwells[0][WFrac] > 0 ? 1: -1; const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac]),dFLimit); const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac]),dFLimit);
primary_variables_[WFrac] = xvar_well_old[WFrac] - dx2_limited; primary_variables_[WFrac] = xvar_well_old[WFrac] - dx2_limited;
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int sign3 = dwells[0][GFrac] > 0 ? 1: -1; const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac]),dFLimit); const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac]),dFLimit);
primary_variables_[GFrac] = xvar_well_old[GFrac] - dx3_limited; primary_variables_[GFrac] = xvar_well_old[GFrac] - dx3_limited;
@ -758,15 +767,15 @@ namespace Opm
primary_variables_[SFrac] = xvar_well_old[SFrac] - dx4_limited; primary_variables_[SFrac] = xvar_well_old[SFrac] - dx4_limited;
} }
assert(active()[ Oil ]); assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
F[pu.phase_pos[Oil]] = 1.0; F[pu.phase_pos[Oil]] = 1.0;
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] = primary_variables_[WFrac]; F[pu.phase_pos[Water]] = primary_variables_[WFrac];
F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]]; F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]];
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] = primary_variables_[GFrac]; F[pu.phase_pos[Gas]] = primary_variables_[GFrac];
F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]]; F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]];
} }
@ -777,9 +786,9 @@ namespace Opm
F[pu.phase_pos[Oil]] -= F_solvent; F[pu.phase_pos[Oil]] -= F_solvent;
} }
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if (F[Water] < 0.0) { if (F[Water] < 0.0) {
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]); F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]);
} }
if (has_solvent) { if (has_solvent) {
@ -790,9 +799,9 @@ namespace Opm
} }
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (F[pu.phase_pos[Gas]] < 0.0) { if (F[pu.phase_pos[Gas]] < 0.0) {
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]); F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]);
} }
if (has_solvent) { if (has_solvent) {
@ -804,10 +813,10 @@ namespace Opm
} }
if (F[pu.phase_pos[Oil]] < 0.0) { if (F[pu.phase_pos[Oil]] < 0.0) {
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]); F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]);
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]); F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]);
} }
if (has_solvent) { if (has_solvent) {
@ -816,10 +825,10 @@ namespace Opm
F[pu.phase_pos[Oil]] = 0.0; F[pu.phase_pos[Oil]] = 0.0;
} }
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
primary_variables_[WFrac] = F[pu.phase_pos[Water]]; primary_variables_[WFrac] = F[pu.phase_pos[Water]];
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
primary_variables_[GFrac] = F[pu.phase_pos[Gas]]; primary_variables_[GFrac] = F[pu.phase_pos[Gas]];
} }
if(has_solvent) { if(has_solvent) {
@ -876,13 +885,13 @@ namespace Opm
std::vector<double> rates(3, 0.0); // the vfp related only supports three phases for the moment std::vector<double> rates(3, 0.0); // the vfp related only supports three phases for the moment
const Opm::PhaseUsage& pu = phaseUsage(); const Opm::PhaseUsage& pu = phaseUsage();
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ]; rates[ Water ] = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ];
} }
if (active()[ Oil ]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ]; rates[ Oil ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ]; rates[ Gas ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ];
} }
@ -943,13 +952,13 @@ namespace Opm
const Opm::PhaseUsage& pu = phaseUsage(); const Opm::PhaseUsage& pu = phaseUsage();
std::vector<double> rates(3, 0.0); std::vector<double> rates(3, 0.0);
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Water ] ]; rates[ Water ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Water ] ];
} }
if (active()[ Oil ]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Oil ] ]; rates[ Oil ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Oil ] ];
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Gas ] ]; rates[ Gas ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Gas ] ];
} }
@ -999,13 +1008,13 @@ namespace Opm
const Opm::PhaseUsage& pu = phaseUsage(); const Opm::PhaseUsage& pu = phaseUsage();
std::vector<double> rates(3, 0.0); std::vector<double> rates(3, 0.0);
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Water ] ]; rates[ Water ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
} }
if (active()[ Oil ]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Oil ] ]; rates[ Oil ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Gas ] ]; rates[ Gas ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
} }
@ -1096,8 +1105,12 @@ namespace Opm
surf_dens_perf.resize(nperf * num_components_); surf_dens_perf.resize(nperf * num_components_);
const int w = index_of_well_; const int w = index_of_well_;
const bool waterPresent = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
const bool oilPresent = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
const bool gasPresent = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
//rs and rv are only used if both oil and gas is present //rs and rv are only used if both oil and gas is present
if (pu.phase_used[Gas] && pu.phase_used[Oil]) { if (oilPresent && gasPresent) {
rsmax_perf.resize(nperf); rsmax_perf.resize(nperf);
rvmax_perf.resize(nperf); rvmax_perf.resize(nperf);
} }
@ -1114,16 +1127,18 @@ namespace Opm
const double p_avg = (well_state.perfPress()[first_perf_ + perf] + p_above)/2; const double p_avg = (well_state.perfPress()[first_perf_ + perf] + p_above)/2;
const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
if (pu.phase_used[Water]) { if (waterPresent) {
b_perf[ pu.phase_pos[Water] + perf * num_components_] = const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
b_perf[ waterCompIdx + perf * num_components_] =
FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
} }
if (pu.phase_used[Gas]) { if (gasPresent) {
const int gaspos = pu.phase_pos[Gas] + perf * num_components_; const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const int gaspos = gasCompIdx + perf * num_components_;
const int gaspos_well = pu.phase_pos[Gas] + w * pu.num_phases; const int gaspos_well = pu.phase_pos[Gas] + w * pu.num_phases;
if (pu.phase_used[Oil]) { if (oilPresent) {
const int oilpos_well = pu.phase_pos[Oil] + w * pu.num_phases; const int oilpos_well = pu.phase_pos[Oil] + w * pu.num_phases;
const double oilrate = std::abs(well_state.wellRates()[oilpos_well]); //in order to handle negative rates in producers const double oilrate = std::abs(well_state.wellRates()[oilpos_well]); //in order to handle negative rates in producers
rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
@ -1146,10 +1161,11 @@ namespace Opm
} }
} }
if (pu.phase_used[Oil]) { if (oilPresent) {
const int oilpos = pu.phase_pos[Oil] + perf * num_components_; const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const int oilpos = oilCompIdx + perf * num_components_;
const int oilpos_well = pu.phase_pos[Oil] + w * pu.num_phases; const int oilpos_well = pu.phase_pos[Oil] + w * pu.num_phases;
if (pu.phase_used[Gas]) { if (gasPresent) {
rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
const int gaspos_well = pu.phase_pos[Gas] + w * pu.num_phases; const int gaspos_well = pu.phase_pos[Gas] + w * pu.num_phases;
const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - well_state.solventWellRate(w); const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - well_state.solventWellRate(w);
@ -1170,8 +1186,13 @@ namespace Opm
} }
// Surface density. // Surface density.
for (int p = 0; p < pu.num_phases; ++p) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
surf_dens_perf[num_components_ * perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
surf_dens_perf[num_components_ * perf + compIdx] = FluidSystem::referenceDensity( phaseIdx, fs.pvtRegionIndex() );
} }
// We use cell values for solvent injector // We use cell values for solvent injector
@ -1199,7 +1220,6 @@ namespace Opm
const int np = number_of_phases_; const int np = number_of_phases_;
const int nperf = number_of_perforations_; const int nperf = number_of_perforations_;
const int num_comp = num_components_; const int num_comp = num_components_;
const PhaseUsage& phase_usage = phaseUsage();
// 1. Compute the flow (in surface volume units for each // 1. Compute the flow (in surface volume units for each
// component) exiting up the wellbore from each perforation, // component) exiting up the wellbore from each perforation,
@ -1227,12 +1247,9 @@ namespace Opm
// absolute values of the surface rates divided by their sum. // absolute values of the surface rates divided by their sum.
// Then compute volume ratios (formation factors) for each perforation. // Then compute volume ratios (formation factors) for each perforation.
// Finally compute densities for the segments associated with each perforation. // Finally compute densities for the segments associated with each perforation.
const int gaspos = phase_usage.phase_pos[Gas];
const int oilpos = phase_usage.phase_pos[Oil];
std::vector<double> mix(num_comp,0.0); std::vector<double> mix(num_comp,0.0);
std::vector<double> x(num_comp); std::vector<double> x(num_comp);
std::vector<double> surf_dens(num_comp); std::vector<double> surf_dens(num_comp);
std::vector<double> dens(nperf);
for (int perf = 0; perf < nperf; ++perf) { for (int perf = 0; perf < nperf; ++perf) {
// Find component mix. // Find component mix.
@ -1244,28 +1261,36 @@ namespace Opm
} }
} else { } else {
// No flow => use well specified fractions for mix. // No flow => use well specified fractions for mix.
for (int phase = 0; phase < np; ++phase) { for (int component = 0; component < num_comp; ++component) {
mix[phase] = comp_frac_[phase]; if (component < np) {
mix[component] = comp_frac_[ ebosCompIdxToFlowCompIdx(component)];
}
} }
// intialize 0.0 for comIdx >= np; // intialize 0.0 for comIdx >= np;
} }
// Compute volume ratio. // Compute volume ratio.
x = mix; x = mix;
double rs = 0.0;
double rv = 0.0; // Subtract dissolved gas from oil phase and vapporized oil from gas phase
if (!rsmax_perf.empty() && mix[oilpos] > 0.0) { if (FluidSystem::phaseIsActive(FluidSystem::gasCompIdx) && FluidSystem::phaseIsActive(FluidSystem::oilCompIdx)) {
rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]); const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
} const unsigned oilpos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (!rvmax_perf.empty() && mix[gaspos] > 0.0) { double rs = 0.0;
rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); double rv = 0.0;
} if (!rsmax_perf.empty() && mix[oilpos] > 0.0) {
if (rs != 0.0) { rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]);
// Subtract gas in oil from gas mixture }
x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv); if (!rvmax_perf.empty() && mix[gaspos] > 0.0) {
} rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]);
if (rv != 0.0) { }
// Subtract oil in gas from oil mixture if (rs != 0.0) {
x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);; // Subtract gas in oil from gas mixture
x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv);
}
if (rv != 0.0) {
// Subtract oil in gas from oil mixture
x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);;
}
} }
double volrat = 0.0; double volrat = 0.0;
for (int component = 0; component < num_comp; ++component) { for (int component = 0; component < num_comp; ++component) {
@ -1329,8 +1354,6 @@ namespace Opm
StandardWell<TypeTag>:: StandardWell<TypeTag>::
getWellConvergence(const std::vector<double>& B_avg) const getWellConvergence(const std::vector<double>& B_avg) const
{ {
const int np = number_of_phases_;
// the following implementation assume that the polymer is always after the w-o-g phases // the following implementation assume that the polymer is always after the w-o-g phases
// For the polymer case, there is one more mass balance equations of reservoir than wells // For the polymer case, there is one more mass balance equations of reservoir than wells
assert((int(B_avg.size()) == num_components_) || has_polymer); assert((int(B_avg.size()) == num_components_) || has_polymer);
@ -1356,18 +1379,23 @@ namespace Opm
ConvergenceReport report; ConvergenceReport report;
// checking if any NaN or too large residuals found // checking if any NaN or too large residuals found
// TODO: not understand why phase here while component in other places. for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) {
const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); continue;
}
if (std::isnan(well_flux_residual[phaseIdx])) { const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
const std::string& compName = FluidSystem::componentName(canonicalCompIdx);
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx);
if (std::isnan(well_flux_residual[compIdx])) {
report.nan_residual_found = true; report.nan_residual_found = true;
const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName}; const typename ConvergenceReport::ProblemWell problem_well = {name(), compName};
report.nan_residual_wells.push_back(problem_well); report.nan_residual_wells.push_back(problem_well);
} else { } else {
if (well_flux_residual[phaseIdx] > maxResidualAllowed) { if (well_flux_residual[compIdx] > maxResidualAllowed) {
report.too_large_residual_found = true; report.too_large_residual_found = true;
const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName}; const typename ConvergenceReport::ProblemWell problem_well = {name(), compName};
report.too_large_residual_wells.push_back(problem_well); report.too_large_residual_wells.push_back(problem_well);
} }
} }
@ -1405,8 +1433,8 @@ namespace Opm
std::vector<double> perfRates(b_perf.size(),0.0); std::vector<double> perfRates(b_perf.size(),0.0);
for (int perf = 0; perf < nperf; ++perf) { for (int perf = 0; perf < nperf; ++perf) {
for (int phase = 0; phase < np; ++phase) { for (int comp = 0; comp < np; ++comp) {
perfRates[perf * num_components_ + phase] = well_state.perfPhaseRates()[(first_perf_ + perf) * np + phase]; perfRates[perf * num_components_ + comp] = well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(comp)];
} }
if(has_solvent) { if(has_solvent) {
perfRates[perf * num_components_ + contiSolventEqIdx] = well_state.perfRateSolvent()[first_perf_ + perf]; perfRates[perf * num_components_ + contiSolventEqIdx] = well_state.perfRateSolvent()[first_perf_ + perf];
@ -1584,7 +1612,7 @@ namespace Opm
computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s); computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
for(int p = 0; p < np; ++p) { for(int p = 0; p < np; ++p) {
well_flux[p] += cq_s[p].value(); well_flux[ebosCompIdxToFlowCompIdx(p)] += cq_s[p].value();
} }
} }
} }
@ -1633,13 +1661,13 @@ namespace Opm
const Opm::PhaseUsage& pu = phaseUsage(); const Opm::PhaseUsage& pu = phaseUsage();
std::vector<double> rates(3, 0.0); std::vector<double> rates(3, 0.0);
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = potentials[pu.phase_pos[ Water ] ]; rates[ Water ] = potentials[pu.phase_pos[ Water ] ];
} }
if (active()[ Oil ]) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = potentials[pu.phase_pos[ Oil ] ]; rates[ Oil ] = potentials[pu.phase_pos[ Oil ] ];
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = potentials[pu.phase_pos[ Gas ] ]; rates[ Gas ] = potentials[pu.phase_pos[ Gas ] ];
} }
@ -1786,10 +1814,10 @@ namespace Opm
tot_well_rate += scalingFactor(p) * well_state.wellRates()[np*well_index + p]; tot_well_rate += scalingFactor(p) * well_state.wellRates()[np*well_index + p];
} }
if(std::abs(tot_well_rate) > 0) { if(std::abs(tot_well_rate) > 0) {
if (active()[ Water ]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / tot_well_rate; primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / tot_well_rate;
} }
if (active()[ Gas ]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / tot_well_rate ; primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / tot_well_rate ;
} }
if (has_solvent) { if (has_solvent) {
@ -1798,7 +1826,7 @@ namespace Opm
} else { // tot_well_rate == 0 } else { // tot_well_rate == 0
if (well_type_ == INJECTOR) { if (well_type_ == INJECTOR) {
// only single phase injection handled // only single phase injection handled
if (active()[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if (distr[Water] > 0.0) { if (distr[Water] > 0.0) {
primary_variables_[WFrac] = 1.0; primary_variables_[WFrac] = 1.0;
} else { } else {
@ -1806,7 +1834,7 @@ namespace Opm
} }
} }
if (active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (distr[pu.phase_pos[Gas]] > 0.0) { if (distr[pu.phase_pos[Gas]] > 0.0) {
primary_variables_[GFrac] = 1.0 - wsolvent(); primary_variables_[GFrac] = 1.0 - wsolvent();
if (has_solvent) { if (has_solvent) {
@ -1822,10 +1850,10 @@ namespace Opm
// this will happen. // this will happen.
} else if (well_type_ == PRODUCER) { // producers } else if (well_type_ == PRODUCER) { // producers
// TODO: the following are not addressed for the solvent case yet // TODO: the following are not addressed for the solvent case yet
if (active()[Water]) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
primary_variables_[WFrac] = 1.0 / np; primary_variables_[WFrac] = 1.0 / np;
} }
if (active()[Gas]) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
primary_variables_[GFrac] = 1.0 / np; primary_variables_[GFrac] = 1.0 / np;
} }
} else { } else {
@ -1953,7 +1981,8 @@ namespace Opm
if (well_type_ == INJECTOR) { if (well_type_ == INJECTOR) {
// assume fully mixing within injecting wellbore // assume fully mixing within injecting wellbore
const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex()); const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
mob[Water] /= (extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration, /*extrapolate=*/true) ); const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
mob[waterCompIdx] /= (extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration, /*extrapolate=*/true) );
} }
if (PolymerModule::hasPlyshlog()) { if (PolymerModule::hasPlyshlog()) {
@ -1974,10 +2003,11 @@ namespace Opm
material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx); material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
const double swcr = scaled_drainage_info.Swcr; const double swcr = scaled_drainage_info.Swcr;
const EvalWell poro = extendEval(int_quant.porosity()); const EvalWell poro = extendEval(int_quant.porosity());
const EvalWell sw = extendEval(int_quant.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water))); const EvalWell sw = extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
// guard against zero porosity and no water // guard against zero porosity and no water
const EvalWell denom = Opm::max( (area * poro * (sw - swcr)), 1e-12); const EvalWell denom = Opm::max( (area * poro * (sw - swcr)), 1e-12);
EvalWell water_velocity = cq_s[Water] / denom * extendEval(int_quant.fluidState().invB(flowPhaseToEbosPhaseIdx(Water))); const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
EvalWell water_velocity = cq_s[waterCompIdx] / denom * extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
if (PolymerModule::hasShrate()) { if (PolymerModule::hasShrate()) {
// the equation for the water velocity conversion for the wells and reservoir are from different version // the equation for the water velocity conversion for the wells and reservoir are from different version
@ -1988,7 +2018,7 @@ namespace Opm
int_quant.pvtRegionIndex(), int_quant.pvtRegionIndex(),
water_velocity); water_velocity);
// modify the mobility with the shear factor. // modify the mobility with the shear factor.
mob[Water] /= shear_factor; mob[waterCompIdx] /= shear_factor;
} }
} }
} }

15
opm/autodiff/WellInterface.hpp
View File

@ -74,11 +74,11 @@ namespace Opm
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
static const int numEq = BlackoilIndices::numEq; static const int numEq = Indices::numEq;
typedef double Scalar; typedef double Scalar;
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType; typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
@ -91,8 +91,8 @@ namespace Opm
static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
// For the conversion between the surface volume rate and resrevoir voidage rate // For the conversion between the surface volume rate and resrevoir voidage rate
using RateConverterType = RateConverter:: using RateConverterType = RateConverter::
@ -123,7 +123,6 @@ namespace Opm
void setVFPProperties(const VFPProperties* vfp_properties_arg); void setVFPProperties(const VFPProperties* vfp_properties_arg);
virtual void init(const PhaseUsage* phase_usage_arg, virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
const std::vector<double>& depth_arg, const std::vector<double>& depth_arg,
const double gravity_arg, const double gravity_arg,
const int num_cells); const int num_cells);
@ -269,8 +268,6 @@ namespace Opm
bool getAllowCrossFlow() const; bool getAllowCrossFlow() const;
const std::vector<bool>* active_;
const VFPProperties* vfp_properties_; const VFPProperties* vfp_properties_;
double gravity_; double gravity_;
@ -284,13 +281,11 @@ namespace Opm
const int num_components_; const int num_components_;
const std::vector<bool>& active() const;
const PhaseUsage& phaseUsage() const; const PhaseUsage& phaseUsage() const;
int flowPhaseToEbosCompIdx( const int phaseIdx ) const; int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; int ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const;
double wsolvent() const; double wsolvent() const;

70
opm/autodiff/WellInterface_impl.hpp
View File

@ -113,13 +113,11 @@ namespace Opm
void void
WellInterface<TypeTag>:: WellInterface<TypeTag>::
init(const PhaseUsage* phase_usage_arg, init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
const std::vector<double>& /* depth_arg */, const std::vector<double>& /* depth_arg */,
const double gravity_arg, const double gravity_arg,
const int /* num_cells */) const int /* num_cells */)
{ {
phase_usage_ = phase_usage_arg; phase_usage_ = phase_usage_arg;
active_ = active_arg;
gravity_ = gravity_arg; gravity_ = gravity_arg;
} }
@ -186,21 +184,6 @@ namespace Opm
template<typename TypeTag>
const std::vector<bool>&
WellInterface<TypeTag>::
active() const
{
assert(active_);
return *active_;
}
template<typename TypeTag> template<typename TypeTag>
void void
WellInterface<TypeTag>:: WellInterface<TypeTag>::
@ -233,39 +216,32 @@ namespace Opm
flowPhaseToEbosCompIdx( const int phaseIdx ) const flowPhaseToEbosCompIdx( const int phaseIdx ) const
{ {
const auto& pu = phaseUsage(); const auto& pu = phaseUsage();
if (active()[Water] && pu.phase_pos[Water] == phaseIdx) if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && pu.phase_pos[Water] == phaseIdx)
return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && pu.phase_pos[Oil] == phaseIdx)
return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && pu.phase_pos[Gas] == phaseIdx)
return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// for other phases return the index // for other phases return the index
return phaseIdx; return phaseIdx;
} }
template<typename TypeTag> template<typename TypeTag>
int int
WellInterface<TypeTag>:: WellInterface<TypeTag>::
flowPhaseToEbosPhaseIdx( const int phaseIdx ) const ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const
{ {
const auto& pu = phaseUsage(); const auto& pu = phaseUsage();
if (active()[Water] && pu.phase_pos[Water] == phaseIdx) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == compIdx)
return FluidSystem::waterPhaseIdx; return pu.phase_pos[Water];
} if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == compIdx)
if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) { return pu.phase_pos[Oil];
return FluidSystem::oilPhaseIdx; if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == compIdx)
} return pu.phase_pos[Gas];
if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) {
return FluidSystem::gasPhaseIdx;
}
assert(phaseIdx < 3);
// for other phases return the index // for other phases return the index
return phaseIdx; return compIdx;
} }
@ -461,7 +437,7 @@ namespace Opm
const int np = number_of_phases_; const int np = number_of_phases_;
if (econ_production_limits.onMinOilRate()) { if (econ_production_limits.onMinOilRate()) {
assert(active()[Oil]); assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ]; const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
const double min_oil_rate = econ_production_limits.minOilRate(); const double min_oil_rate = econ_production_limits.minOilRate();
if (std::abs(oil_rate) < min_oil_rate) { if (std::abs(oil_rate) < min_oil_rate) {
@ -470,7 +446,7 @@ namespace Opm
} }
if (econ_production_limits.onMinGasRate() ) { if (econ_production_limits.onMinGasRate() ) {
assert(active()[Gas]); assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
const double gas_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ]; const double gas_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ];
const double min_gas_rate = econ_production_limits.minGasRate(); const double min_gas_rate = econ_production_limits.minGasRate();
if (std::abs(gas_rate) < min_gas_rate) { if (std::abs(gas_rate) < min_gas_rate) {
@ -479,8 +455,8 @@ namespace Opm
} }
if (econ_production_limits.onMinLiquidRate() ) { if (econ_production_limits.onMinLiquidRate() ) {
assert(active()[Oil]); assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
assert(active()[Water]); assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ]; const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
const double water_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ]; const double water_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ];
const double liquid_rate = oil_rate + water_rate; const double liquid_rate = oil_rate + water_rate;
@ -517,8 +493,8 @@ namespace Opm
const Opm::PhaseUsage& pu = phaseUsage(); const Opm::PhaseUsage& pu = phaseUsage();
const int well_number = index_of_well_; const int well_number = index_of_well_;
assert(active()[Oil]); assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
assert(active()[Water]); assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
@ -854,11 +830,11 @@ namespace Opm
return distr[phaseIdx]; return distr[phaseIdx];
} }
const auto& pu = phaseUsage(); const auto& pu = phaseUsage();
if (active()[Water] && pu.phase_pos[Water] == phaseIdx) if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && pu.phase_pos[Water] == phaseIdx)
return 1.0; return 1.0;
if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && pu.phase_pos[Oil] == phaseIdx)
return 1.0; return 1.0;
if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && pu.phase_pos[Gas] == phaseIdx)
return 0.01; return 0.01;
if (has_solvent && phaseIdx == contiSolventEqIdx ) if (has_solvent && phaseIdx == contiSolventEqIdx )
return 0.01; return 0.01;