diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 7d80c5e37..567e2db9c 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -29,8 +29,6 @@ #include #include -#include -#include #include #include @@ -43,7 +41,6 @@ #include #include #include -#include #include #include #include @@ -144,7 +141,6 @@ namespace Opm { , grid_(ebosSimulator_.gridManager().grid()) , istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) ) , phaseUsage_(phaseUsageFromDeck(eclState())) - , active_(detail::activePhases(phaseUsage_)) , has_disgas_(FluidSystem::enableDissolvedGas()) , has_vapoil_(FluidSystem::enableVaporizedOil()) , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) @@ -153,7 +149,7 @@ namespace Opm { , well_model_ (well_model) , terminal_output_ (terminal_output) , current_relaxation_(1.0) - , dx_old_(AutoDiffGrid::numCells(grid_)) + , dx_old_(UgGridHelpers::numCells(grid_)) { // compute global sum of number of cells global_nc_ = detail::countGlobalCells(grid_); @@ -276,7 +272,7 @@ namespace Opm { //residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ; // Compute the nonlinear update. - const int nc = AutoDiffGrid::numCells(grid_); + const int nc = UgGridHelpers::numCells(grid_); BVector x(nc); try { @@ -596,8 +592,6 @@ namespace Opm { /// \param[in, out] well_state well state variables void updateState(const BVector& dx) { - using namespace Opm::AutoDiffGrid; - const auto& ebosProblem = ebosSimulator_.problem(); unsigned numSwitched = 0; @@ -626,8 +620,8 @@ namespace Opm { p = std::max(p, 0.0); // Saturation updates. - const double dsw = active_[Water] ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0; - const double dxvar = active_[Gas] ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0; + const double dsw = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0; + const double dxvar = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0; double dso = 0.0; double dsg = 0.0; @@ -669,12 +663,12 @@ namespace Opm { satScaleFactor = dsMax()/maxVal; } - if (active_[Water]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { double& sw = priVars[Indices::waterSaturationIdx]; sw -= satScaleFactor * dsw; } - if (active_[Gas]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { double& sg = priVars[Indices::compositionSwitchIdx]; sg -= satScaleFactor * dsg; @@ -693,7 +687,7 @@ namespace Opm { } // 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); const double drmaxrel = drMaxRel(); if (has_disgas_) { @@ -809,8 +803,7 @@ namespace Opm { const double tol_mb = param_.tolerance_mb_; const double tol_cnv = param_.tolerance_cnv_; - const int np = numPhases(); - const int numComp = numComponents(); + const int numComp = numEq; Vector R_sum(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 ); pvSumLocal += pvValue; - for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx ) + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx); - const int ebosCompIdx = flowPhaseToEbosCompIdx(phaseIdx); + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } - B_avg[ phaseIdx ] += 1.0 / fs.invB(ebosPhaseIdx).value(); - const auto R2 = ebosResid[cell_idx][ebosCompIdx]; + const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - R_sum[ phaseIdx ] += R2; - maxCoeff[ phaseIdx ] = std::max( maxCoeff[ phaseIdx ], std::abs( R2 ) / pvValue ); + B_avg[ compIdx ] += 1.0 / fs.invB(phaseIdx).value(); + const auto R2 = ebosResid[cell_idx][compIdx]; + + R_sum[ compIdx ] += R2; + maxCoeff[ compIdx ] = std::max( maxCoeff[ compIdx ], std::abs( R2 ) / pvValue ); } if ( has_solvent_ ) { @@ -914,9 +910,15 @@ namespace Opm { std::string msg = "Iter"; std::vector< std::string > key( numComp ); - for (int phaseIdx = 0; phaseIdx < numPhases(); ++phaseIdx) { - const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); - key[ phaseIdx ] = std::toupper( phaseName.front() ); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + 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_) { key[ solventSaturationIdx ] = "S"; @@ -949,16 +951,25 @@ namespace Opm { OpmLog::debug(ss.str()); } - for (int phaseIdx = 0; phaseIdx < numPhases(); ++phaseIdx) { - const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; - if (std::isnan(mass_balance_residual[phaseIdx]) - || std::isnan(CNV[phaseIdx])) { - OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName); + const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx); + const std::string& compName = FluidSystem::componentName(canonicalCompIdx); + 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() - || CNV[phaseIdx] > maxResidualAllowed()) { - OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName); + if (mass_balance_residual[compIdx] > maxResidualAllowed() + || CNV[compIdx] > maxResidualAllowed()) { + 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; } - 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 template std::vector > @@ -1001,7 +997,6 @@ namespace Opm { const auto& comm = grid_.comm(); const auto& gridView = ebosSimulator().gridView(); const int nc = gridView.size(/*codim=*/0); - const int maxnp = Opm::BlackoilPhases::MaxNumPhases; int ntFip = *std::max_element(fipnum.begin(), fipnum.end()); ntFip = comm.max(ntFip); @@ -1044,18 +1039,20 @@ namespace Opm { ebosSimulator_.model().dofTotalVolume(cellIdx) * intQuants.porosity().value(); - for (int phase = 0; phase < maxnp; ++phase) { - const double b = fs.invB(flowPhaseToEbosPhaseIdx(phase)).value(); - const double s = fs.saturation(flowPhaseToEbosPhaseIdx(phase)).value(); - - fip_.fip[phase][cellIdx] = b * s * pv; - - if (active_[ phase ]) { - regionValues[regionIdx][phase] += fip_.fip[phase][cellIdx]; + for (unsigned phase = 0; phase < FluidSystem::numPhases; ++phase) { + if (!FluidSystem::phaseIsActive(phase)) { + continue; } + + 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 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]; @@ -1461,10 +1458,6 @@ namespace Opm { const Grid& grid_; const ISTLSolverType* istlSolver_; const PhaseUsage phaseUsage_; - // For each canonical phase -> true if active - const std::vector active_; - // Size = # active phases. Maps active -> canonical phase indices. - const std::vector cells_; // All grid cells const bool has_disgas_; const bool has_vapoil_; const bool has_solvent_; @@ -1494,29 +1487,14 @@ namespace Opm { const BlackoilWellModel& wellModel() const { return well_model_; } - int flowPhaseToEbosCompIdx( const int phaseIdx ) const + int ebosPhaseToFlowCanonicalPhaseIdx( const int phaseIdx ) const { - const auto& pu = phaseUsage_; - if (active_[Water] && pu.phase_pos[Water] == phaseIdx) - return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) - return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) - 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; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::waterPhaseIdx == phaseIdx) + return Water; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::oilPhaseIdx == phaseIdx) + return Oil; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::gasPhaseIdx == phaseIdx) + return Gas; assert(phaseIdx < 3); // for other phases return the index diff --git a/opm/autodiff/BlackoilWellModel.hpp b/opm/autodiff/BlackoilWellModel.hpp index f81fb6294..b59342a5f 100644 --- a/opm/autodiff/BlackoilWellModel.hpp +++ b/opm/autodiff/BlackoilWellModel.hpp @@ -69,19 +69,15 @@ namespace Opm { typedef WellStateFullyImplicitBlackoil WellState; typedef BlackoilModelParameters ModelParameters; - static const int Water = WellInterface::Water; - static const int Oil = WellInterface::Oil; - static const int Gas = WellInterface::Gas; - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; 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, Scalar) Scalar; - static const int numEq = BlackoilIndices::numEq; - static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; + static const int numEq = Indices::numEq; + static const int solventSaturationIdx = Indices::solventSaturationIdx; // TODO: where we should put these types, WellInterface or Well Model? // or there is some other strategy, like TypeTag @@ -176,7 +172,6 @@ namespace Opm { bool has_polymer_; std::vector pvt_region_idx_; PhaseUsage phase_usage_; - std::vector active_; size_t global_nc_; // the number of the cells in the local grid size_t number_of_cells_; @@ -241,9 +236,6 @@ namespace Opm { int numPhases() const; - - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; - void resetWellControlFromState() const; void assembleWellEq(const double dt, diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index 5f350cff8..87b0a989c 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -18,11 +18,6 @@ namespace Opm { const auto& eclState = ebosSimulator_.gridManager().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(); // 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 // optimize the usage of the following several member variables 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 @@ -378,26 +373,6 @@ namespace Opm { - - template - int - BlackoilWellModel:: - 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 void BlackoilWellModel:: @@ -1004,8 +979,6 @@ namespace Opm { BlackoilWellModel:: computeAverageFormationFactor(std::vector& B_avg) const { - const int np = numPhases(); - const auto& grid = ebosSimulator_.gridManager().grid(); const auto& gridView = grid.leafGridView(); ElementContext elemCtx(ebosSimulator_); @@ -1020,12 +993,16 @@ namespace Opm { const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); 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 ]; - const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx); + if (!FluidSystem::phaseIsActive(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_) { auto& B = B_avg[solventSaturationIdx]; diff --git a/opm/autodiff/MultisegmentWell.hpp b/opm/autodiff/MultisegmentWell.hpp index e14626884..ac55dd610 100644 --- a/opm/autodiff/MultisegmentWell.hpp +++ b/opm/autodiff/MultisegmentWell.hpp @@ -40,7 +40,7 @@ namespace Opm using typename Base::FluidSystem; using typename Base::ModelParameters; using typename Base::MaterialLaw; - using typename Base::BlackoilIndices; + using typename Base::Indices; 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: gas is always there? how about oil water 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 WFrac = gasoil? -1000: 1; static const int GFrac = gasoil? 1 : 2; @@ -106,7 +106,6 @@ namespace Opm const int num_components); virtual void init(const PhaseUsage* phase_usage_arg, - const std::vector* active_arg, const std::vector& depth_arg, const double gravity_arg, const int num_cells); @@ -187,11 +186,10 @@ namespace Opm using Base::num_components_; // protected functions from the Base class - using Base::active; using Base::phaseUsage; using Base::name; - using Base::flowPhaseToEbosPhaseIdx; using Base::flowPhaseToEbosCompIdx; + using Base::ebosCompIdxToFlowCompIdx; using Base::getAllowCrossFlow; using Base::scalingFactor; @@ -284,7 +282,7 @@ namespace Opm // fraction value of the primary variables // should we just use member variables to store them instead of calculating them again and again - EvalWell volumeFraction(const int seg, const 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 EvalWell volumeFractionScaled(const int seg, const int comp_idx) const; diff --git a/opm/autodiff/MultisegmentWell_impl.hpp b/opm/autodiff/MultisegmentWell_impl.hpp index d78eed4a8..d89cc12d8 100644 --- a/opm/autodiff/MultisegmentWell_impl.hpp +++ b/opm/autodiff/MultisegmentWell_impl.hpp @@ -105,13 +105,11 @@ namespace Opm void MultisegmentWell:: init(const PhaseUsage* phase_usage_arg, - const std::vector* active_arg, const std::vector& depth_arg, const double gravity_arg, const int num_cells) { - Base::init(phase_usage_arg, active_arg, - depth_arg, gravity_arg, num_cells); + Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells); // TODO: for StandardWell, we need to update the perf depth here using depth_arg. // for MultisegmentWell, it is much more complicated. @@ -274,13 +272,13 @@ namespace Opm /* const Opm::PhaseUsage& pu = phaseUsage(); std::vector 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 ] ]; } - if (active()[ Oil ]) { + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { 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 ] ]; } */ @@ -433,13 +431,13 @@ namespace Opm // TODO: the report can not handle the segment number yet. if (std::isnan(flux_residual)) { report.nan_residual_found = true; - const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx)); - const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name}; + const auto& compName = FluidSystem::componentName(Indices::activeToCanonicalComponentIndex(eq_idx)); + const typename ConvergenceReport::ProblemWell problem_well = {name(), compName}; report.nan_residual_wells.push_back(problem_well); } else if (flux_residual > param_.max_residual_allowed_) { report.too_large_residual_found = true; - const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx)); - const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name}; + const auto& compName = FluidSystem::componentName(Indices::activeToCanonicalComponentIndex(eq_idx)); + const typename ConvergenceReport::ProblemWell problem_well = {name(), compName}; report.nan_residual_wells.push_back(problem_well); } else { // it is a normal residual if (flux_residual > maximum_residual[eq_idx]) { @@ -578,11 +576,11 @@ namespace Opm primary_variables_[seg][GTotal] = total_seg_rate; if (std::abs(total_seg_rate) > 0.) { - if (active()[Water]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { 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; } - if (active()[Gas]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { 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; } @@ -590,7 +588,7 @@ namespace Opm if (well_type_ == INJECTOR) { // only single phase injection handled 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) { primary_variables_[seg][WFrac] = 1.0; } else { @@ -598,7 +596,7 @@ namespace Opm } } - if (active()[Gas]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (distr[pu.phase_pos[Gas]] > 0.0) { primary_variables_[seg][GFrac] = 1.0; } else { @@ -606,11 +604,11 @@ namespace Opm } } } else if (well_type_ == PRODUCER) { // producers - if (active()[Water]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { 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_; } } @@ -743,13 +741,13 @@ namespace Opm const std::vector > old_primary_variables = primary_variables_; 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 double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit); primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited; } - if (active()[ Gas ]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { const int sign = dwells[seg][GFrac] > 0. ? 1 : -1; const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit); primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited; @@ -868,25 +866,24 @@ namespace Opm template typename MultisegmentWell::EvalWell MultisegmentWell:: - 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]; } - 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]; } // Oil fraction EvalWell oil_fraction = 1.0; - if (active()[Water]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { oil_fraction -= primary_variables_evaluation_[seg][WFrac]; } - if (active()[Gas]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { oil_fraction -= primary_variables_evaluation_[seg][GFrac]; } /* if (has_solvent) { @@ -898,7 +895,6 @@ namespace Opm - template typename MultisegmentWell::EvalWell MultisegmentWell:: @@ -907,7 +903,7 @@ namespace Opm // For reservoir rate control, the distr in well control is used for the // rate conversion coefficients. For the injection well, only the distr of the injection // phase is not zero. - const double scale = scalingFactor(comp_idx); + const double scale = scalingFactor(ebosCompIdxToFlowCompIdx(comp_idx)); if (scale > 0.) { return volumeFraction(seg, comp_idx) / scale; } @@ -965,9 +961,13 @@ namespace Opm // not using number_of_phases_ because of solvent std::vector b_perfcells(num_components_, 0.0); - for (int phase = 0; phase < number_of_phases_; ++phase) { - const int phase_idx_ebos = flowPhaseToEbosPhaseIdx(phase); - b_perfcells[phase] = extendEval(fs.invB(phase_idx_ebos)); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + b_perfcells[compIdx] = extendEval(fs.invB(phaseIdx)); } // 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 const EvalWell drawdown = (pressure_cell + cell_perf_press_diff) - (segment_pressure + perf_seg_press_diff); - const Opm::PhaseUsage& pu = phaseUsage(); - // producing perforations if ( drawdown > 0.0) { // Do nothing is crossflow is not allowed @@ -994,13 +992,13 @@ namespace Opm cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p; } - if (active()[Oil] && active()[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const EvalWell cq_s_oil = cq_s[oilpos]; - const EvalWell cq_s_gas = cq_s[gaspos]; - cq_s[gaspos] += rs * cq_s_oil; - cq_s[oilpos] += rv * cq_s_gas; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const EvalWell cq_s_oil = cq_s[oilCompIdx]; + const EvalWell cq_s_gas = cq_s[gasCompIdx]; + cq_s[gasCompIdx] += rs * cq_s_oil; + cq_s[oilCompIdx] += rv * cq_s_gas; } } else { // injecting perforations // Do nothing if crossflow is not allowed @@ -1019,14 +1017,14 @@ namespace Opm // compute volume ratio between connection and at standard conditions EvalWell volume_ratio = 0.0; - if (active()[Water]) { - const int watpos = pu.phase_pos[Water]; - volume_ratio += cmix_s[watpos] / b_perfcells[watpos]; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx]; } - if (active()[Oil] && active()[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); // 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 @@ -1038,19 +1036,19 @@ namespace Opm << " with rs " << rs << " and rv " << rv); } - const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d; - volume_ratio += tmp_oil / b_perfcells[oilpos]; + const EvalWell tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d; + volume_ratio += tmp_oil / b_perfcells[oilCompIdx]; - const EvalWell tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d; - volume_ratio += tmp_gas / b_perfcells[gaspos]; + const EvalWell tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d; + volume_ratio += tmp_gas / b_perfcells[gasCompIdx]; } else { // not having gas and oil at the same time - if (active()[Oil]) { - const int oilpos = pu.phase_pos[Oil]; - volume_ratio += cmix_s[oilpos] / b_perfcells[oilpos]; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx]; } - if (active()[Gas]) { - const int gaspos = pu.phase_pos[Gas]; - volume_ratio += cmix_s[gaspos] / b_perfcells[gaspos]; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx]; } } // injecting connections total volumerates at standard conditions @@ -1110,14 +1108,17 @@ namespace Opm pvt_region_index = fs.pvtRegionIndex(); } - std::vector surf_dens(number_of_phases_); + std::vector surf_dens(num_components_); // Surface density. - // not using num_comp here is because solvent can be component - for (int phase = 0; phase < number_of_phases_; ++phase) { - surf_dens[phase] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx(phase), pvt_region_index ); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + 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) { // the compostion of the components inside wellbore under surface condition std::vector mix_s(num_components_, 0.0); @@ -1126,94 +1127,95 @@ namespace Opm } std::vector b(num_components_, 0.0); - // it is the phase viscosities asked for - std::vector visc(number_of_phases_, 0.0); + std::vector visc(num_components_, 0.0); + const EvalWell seg_pressure = getSegmentPressure(seg); - if (pu.phase_used[Water]) { - const int water_pos = pu.phase_pos[Water]; - b[water_pos] = + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + b[waterCompIdx] = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - visc[water_pos] = + visc[waterCompIdx] = FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure); } EvalWell rv(0.0); // gas phase - if (pu.phase_used[Gas]) { - const int gaspos = pu.phase_pos[Gas]; - if (pu.phase_used[Oil]) { - const int oilpos = pu.phase_pos[Oil]; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure); - if (mix_s[oilpos] > 0.0) { - if (mix_s[gaspos] > 0.0) { - rv = mix_s[oilpos] / mix_s[gaspos]; + if (mix_s[oilCompIdx] > 0.0) { + if (mix_s[gasCompIdx] > 0.0) { + rv = mix_s[oilCompIdx] / mix_s[gasCompIdx]; } if (rv > rvmax) { rv = rvmax; } - b[gaspos] = + b[gasCompIdx] = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv); - visc[gaspos] = + visc[gasCompIdx] = FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv); } else { // no oil exists - b[gaspos] = + b[gasCompIdx] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - visc[gaspos] = + visc[gasCompIdx] = FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); } } else { // no Liquid phase // it is the same with zero mix_s[Oil] - b[gaspos] = + b[gasCompIdx] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - visc[gaspos] = + visc[gasCompIdx] = FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); } } EvalWell rs(0.0); // oil phase - if (pu.phase_used[Oil]) { - const int oilpos = pu.phase_pos[Oil]; - if (pu.phase_used[Oil]) { - const int gaspos = pu.phase_pos[Gas]; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure); - if (mix_s[gaspos] > 0.0) { - if (mix_s[oilpos] > 0.0) { - rs = mix_s[gaspos] / mix_s[oilpos]; + if (mix_s[gasCompIdx] > 0.0) { + if (mix_s[oilCompIdx] > 0.0) { + rs = mix_s[gasCompIdx] / mix_s[oilCompIdx]; } if (rs > rsmax) { rs = rsmax; } - b[oilpos] = + b[oilCompIdx] = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs); - visc[oilpos] = + visc[oilCompIdx] = FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs); } else { // no oil exists - b[oilpos] = + b[oilCompIdx] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - visc[oilpos] = + visc[oilCompIdx] = FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); } } else { // no Liquid phase // it is the same with zero mix_s[Oil] - b[oilpos] = + b[oilCompIdx] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - visc[oilpos] = + visc[oilCompIdx] = FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); } } std::vector mix(mix_s); - if (pu.phase_used[Oil] && pu.phase_used[Gas]) { - const int gaspos = pu.phase_pos[Gas]; - const int oilpos = pu.phase_pos[Oil]; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + 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? - 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.; // calculate the average viscosity - for (int p = 0; p < number_of_phases_; ++p) { - const EvalWell phase_fraction = mix[p] / b[p] / volrat; - segment_viscosities_[seg] += visc[p] * phase_fraction; + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { + const EvalWell comp_fraction = mix[comp_idx] / b[comp_idx] / volrat; + segment_viscosities_[seg] += visc[comp_idx] * comp_fraction; } EvalWell density(0.0); @@ -1237,9 +1239,9 @@ namespace Opm // calculate the mass rates segment_mass_rates_[seg] = 0.; - for (int phase = 0; phase < number_of_phases_; ++phase) { - const EvalWell rate = getSegmentRate(seg, phase); - segment_mass_rates_[seg] += rate * surf_dens[phase]; + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { + const EvalWell rate = getSegmentRate(seg, comp_idx); + segment_mass_rates_[seg] += rate * surf_dens[comp_idx]; } } } @@ -1293,7 +1295,6 @@ namespace Opm std::vector& mob) const { // 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]; assert (int(mob.size()) == num_components_); const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); @@ -1305,9 +1306,13 @@ namespace Opm 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 - for (int phase = 0; phase < np; ++phase) { - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx)); } // if (has_solvent) { // mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); @@ -1322,9 +1327,13 @@ namespace Opm materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); // compute the mobility - for (int phase = 0; phase < np; ++phase) { - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx)); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + 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 for (int phase = 0; phase < number_of_phases_; ++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; } } @@ -1380,7 +1389,7 @@ namespace Opm const EvalWell G_total = getSegmentGTotal(0); for (int phase = 0; phase < number_of_phases_; ++phase) { 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 @@ -1394,7 +1403,7 @@ namespace Opm const double* distr = well_controls_get_current_distr(well_controls_); for (int phase = 0; phase < number_of_phases_; ++phase) { 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_); @@ -1546,26 +1555,26 @@ namespace Opm std::vector fractions(number_of_phases_, 0.0); - assert( active()[Oil] ); + assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ); const int oil_pos = pu.phase_pos[Oil]; fractions[oil_pos] = 1.0; - if ( active()[Water] ) { + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { const int water_pos = pu.phase_pos[Water]; fractions[water_pos] = primary_variables_[seg][WFrac]; fractions[oil_pos] -= fractions[water_pos]; } - if ( active()[Gas] ) { + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { const int gas_pos = pu.phase_pos[Gas]; fractions[gas_pos] = primary_variables_[seg][GFrac]; fractions[oil_pos] -= fractions[gas_pos]; } - if (active()[Water]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { const int water_pos = pu.phase_pos[Water]; 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[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]; 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[oil_pos] /= (1.0 - fractions[gas_pos]); @@ -1585,20 +1594,20 @@ namespace Opm } if (fractions[oil_pos] < 0.0) { - if ( active()[Water] ) { + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { 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[oil_pos] = 0.0; } - if ( active()[Water] ) { + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { 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]]; } } @@ -1613,20 +1622,20 @@ namespace Opm updateWellStateFromPrimaryVariables(WellState& well_state) const { const PhaseUsage& pu = phaseUsage(); - assert( active()[Oil] ); + assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ); const int oil_pos = pu.phase_pos[Oil]; for (int seg = 0; seg < numberOfSegments(); ++seg) { std::vector fractions(number_of_phases_, 0.0); fractions[oil_pos] = 1.0; - if ( active()[Water] ) { + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { const int water_pos = pu.phase_pos[Water]; fractions[water_pos] = primary_variables_[seg][WFrac]; fractions[oil_pos] -= fractions[water_pos]; } - if ( active()[Gas] ) { + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { const int gas_pos = pu.phase_pos[Gas]; fractions[gas_pos] = primary_variables_[seg][GFrac]; fractions[oil_pos] -= fractions[gas_pos]; @@ -1814,9 +1823,7 @@ namespace Opm if (!only_wells) { // subtract sum of component fluxes in the reservoir equation. // need to consider the efficiency factor - // TODO: the name of the function flowPhaseToEbosCompIdx is prolematic, since the input - // is a component index :D - ebosResid[cell_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.value(); + ebosResid[cell_idx][comp_idx] -= cq_s_effective.value(); } // 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) { if (!only_wells) { // 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 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) { if (!only_wells) { // 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); } } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index ce24dd942..9964bea19 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -29,8 +29,6 @@ #include #include #include -#include -#include #include #include #include diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 617346fbc..22fa6bfa8 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -46,7 +46,7 @@ namespace Opm using typename Base::FluidSystem; using typename Base::MaterialLaw; using typename Base::ModelParameters; - using typename Base::BlackoilIndices; + using typename Base::Indices; using typename Base::PolymerModule; 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 // 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 WFrac = gasoil? -1000: 1; static const int GFrac = gasoil? 1: 2; @@ -117,7 +117,6 @@ namespace Opm const int num_components); virtual void init(const PhaseUsage* phase_usage_arg, - const std::vector* active_arg, const std::vector& depth_arg, const double gravity_arg, const int num_cells); @@ -163,9 +162,8 @@ namespace Opm // protected functions from the Base class using Base::getAllowCrossFlow; using Base::phaseUsage; - using Base::active; - using Base::flowPhaseToEbosPhaseIdx; using Base::flowPhaseToEbosCompIdx; + using Base::ebosCompIdxToFlowCompIdx; using Base::wsolvent; using Base::wpolymer; using Base::wellHasTHPConstraints; @@ -233,7 +231,7 @@ namespace Opm EvalWell wellVolumeFractionScaled(const int phase) const; - EvalWell wellVolumeFraction(const int phase) const; + EvalWell wellVolumeFraction(const unsigned compIdx) const; EvalWell wellSurfaceVolumeFraction(const int phase) const; diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 00b5d1be3..7c9406993 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -49,13 +49,11 @@ namespace Opm void StandardWell:: init(const PhaseUsage* phase_usage_arg, - const std::vector* active_arg, const std::vector& depth_arg, const double gravity_arg, const int num_cells) { - Base::init(phase_usage_arg, active_arg, - depth_arg, gravity_arg, num_cells); + Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells); perf_depth_.resize(number_of_perforations_, 0.); for (int perf = 0; perf < number_of_perforations_; ++perf) { @@ -138,13 +136,13 @@ namespace Opm const Opm::PhaseUsage& pu = phaseUsage(); std::vector rates(3, 0.0); - if (active()[ Water ]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { rates[ Water ]= getQs(pu.phase_pos[ Water]); } - if (active()[ Oil ]) { + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { rates[ Oil ] = getQs(pu.phase_pos[ Oil ]); } - if (active()[ Gas ]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { rates[ Gas ] = getQs(pu.phase_pos[ Gas ]); } return calculateBhpFromThp(rates, control); @@ -170,6 +168,7 @@ namespace Opm assert(comp_idx < num_components_); const auto pu = phaseUsage(); + const int legacyCompIdx = ebosCompIdxToFlowCompIdx(comp_idx); // TODO: the formulation for the injectors decides it only work with single phase // surface rate injection control. Improvement will be required. @@ -180,10 +179,10 @@ namespace Opm double comp_frac = 0.0; if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent(); - } else if (comp_idx == pu.phase_pos[ Gas ]) { - comp_frac = comp_frac_[comp_idx] * (1.0 - wsolvent()); + } else if (legacyCompIdx == pu.phase_pos[ Gas ]) { + comp_frac = comp_frac_[legacyCompIdx] * (1.0 - wsolvent()); } else { - comp_frac = comp_frac_[comp_idx]; + comp_frac = comp_frac_[legacyCompIdx]; } if (comp_frac == 0.0) { return qs; //zero @@ -197,7 +196,7 @@ namespace Opm return qs; } - const double comp_frac = comp_frac_[comp_idx]; + const double comp_frac = comp_frac_[legacyCompIdx]; if (comp_frac == 0.0) { return qs; } @@ -242,15 +241,16 @@ namespace Opm 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) { // for GRAT controlled wells solvent is included in the target wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(contiSolventEqIdx); } - if (comp_idx == phase_under_control) { - if (has_solvent && phase_under_control == Gas) { - qs.setValue(target_rate * wellVolumeFractionScaled(Gas).value() / wellVolumeFractionScaledPhaseUnderControl.value() ); + if (comp_idx == compIdx_under_control) { + if (has_solvent && compIdx_under_control == FluidSystem::gasCompIdx) { + qs.setValue(target_rate * wellVolumeFractionScaled(compIdx_under_control).value() / wellVolumeFractionScaledPhaseUnderControl.value() ); return qs; } qs.setValue(target_rate); @@ -270,8 +270,9 @@ namespace Opm if (num_phases_under_rate_control == 2) { EvalWell combined_volume_fraction = 0.; for (int p = 0; p < np; ++p) { + const unsigned compIdxTmp = flowPhaseToEbosCompIdx(p); if (distr[p] == 1.0) { - combined_volume_fraction += wellVolumeFractionScaled(p); + combined_volume_fraction += wellVolumeFractionScaled(compIdxTmp); } } return (target_rate * wellVolumeFractionScaled(comp_idx) / combined_volume_fraction); @@ -303,7 +304,8 @@ namespace Opm wellVolumeFractionScaled(const int compIdx) const { - const double scal = scalingFactor(compIdx); + const int legacyCompIdx = ebosCompIdxToFlowCompIdx(compIdx); + const double scal = scalingFactor(legacyCompIdx); if (scal > 0) return wellVolumeFraction(compIdx) / scal; @@ -318,14 +320,13 @@ namespace Opm template typename StandardWell::EvalWell StandardWell:: - wellVolumeFraction(const int compIdx) const + wellVolumeFraction(const unsigned compIdx) const { - const auto& pu = phaseUsage(); - if (active()[Water] && compIdx == pu.phase_pos[Water]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { 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]; } @@ -335,11 +336,11 @@ namespace Opm // Oil fraction EvalWell well_fraction = 1.0; - if (active()[Water]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { well_fraction -= primary_variables_evaluation_[WFrac]; } - if (active()[Gas]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { well_fraction -= primary_variables_evaluation_[GFrac]; } if (has_solvent) { @@ -396,21 +397,22 @@ namespace Opm const double Tw, const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector& cq_s) const { - const Opm::PhaseUsage& pu = phaseUsage(); - const int np = number_of_phases_; std::vector cmix_s(num_components_,0.0); for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx); } const auto& fs = intQuants.fluidState(); - const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); const EvalWell rs = extendEval(fs.Rs()); const EvalWell rv = extendEval(fs.Rv()); std::vector b_perfcells_dense(num_components_, 0.0); - for (int phase = 0; phase < np; ++phase) { - const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + b_perfcells_dense[compIdx] = extendEval(fs.invB(phaseIdx)); } if (has_solvent) { b_perfcells_dense[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor()); @@ -433,13 +435,13 @@ namespace Opm cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p; } - if (active()[Oil] && active()[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const EvalWell cq_sOil = cq_s[oilpos]; - const EvalWell cq_sGas = cq_s[gaspos]; - cq_s[gaspos] += rs * cq_sOil; - cq_s[oilpos] += rv * cq_sGas; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const EvalWell cq_sOil = cq_s[oilCompIdx]; + const EvalWell cq_sGas = cq_s[gasCompIdx]; + cq_s[gasCompIdx] += rs * cq_sOil; + cq_s[oilCompIdx] += rv * cq_sGas; } } else { @@ -459,19 +461,18 @@ namespace Opm // compute volume ratio between connection at standard conditions EvalWell volumeRatio = 0.0; - if (active()[Water]) { - const int watpos = pu.phase_pos[Water]; - volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx]; } if (has_solvent) { volumeRatio += cmix_s[contiSolventEqIdx] / b_perfcells_dense[contiSolventEqIdx]; } - if (active()[Oil] && active()[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); // Incorporate RS/RV factors if both oil and gas active const EvalWell d = 1.0 - rv * rs; @@ -480,22 +481,22 @@ namespace Opm << " 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 " <& mob) const { - const int np = number_of_phases_; const int cell_idx = well_cells_[perf]; assert (int(mob.size()) == num_components_); const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); @@ -683,9 +684,13 @@ namespace Opm 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 - for (int phase = 0; phase < np; ++phase) { - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx)); } if (has_solvent) { mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); @@ -700,9 +705,13 @@ namespace Opm materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); // compute the mobility - for (int phase = 0; phase < np; ++phase) { - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx)); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + 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? @@ -713,7 +722,7 @@ namespace Opm // modify the water mobility if polymer is present if (has_polymer) { - if (!active()[Water]) { + if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { 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) std::vector F(np,0.0); - if (active()[ Water ]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { const int sign2 = dwells[0][WFrac] > 0 ? 1: -1; const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac]),dFLimit); 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 double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac]),dFLimit); primary_variables_[GFrac] = xvar_well_old[GFrac] - dx3_limited; @@ -758,15 +767,15 @@ namespace Opm primary_variables_[SFrac] = xvar_well_old[SFrac] - dx4_limited; } - assert(active()[ Oil ]); + assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)); 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[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[Oil]] -= F[pu.phase_pos[Gas]]; } @@ -777,9 +786,9 @@ namespace Opm F[pu.phase_pos[Oil]] -= F_solvent; } - if (active()[ Water ]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { 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]]); } 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 (active()[ Water ]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]); } if (has_solvent) { @@ -804,10 +813,10 @@ namespace Opm } 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]]); } - if (active()[ Gas ]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]); } if (has_solvent) { @@ -816,10 +825,10 @@ namespace Opm F[pu.phase_pos[Oil]] = 0.0; } - if (active()[ Water ]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { primary_variables_[WFrac] = F[pu.phase_pos[Water]]; } - if (active()[ Gas ]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { primary_variables_[GFrac] = F[pu.phase_pos[Gas]]; } if(has_solvent) { @@ -876,13 +885,13 @@ namespace Opm std::vector rates(3, 0.0); // the vfp related only supports three phases for the moment 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 ] ]; } - if (active()[ Oil ]) { + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { 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 ] ]; } @@ -943,13 +952,13 @@ namespace Opm const Opm::PhaseUsage& pu = phaseUsage(); std::vector 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 ] ]; } - if (active()[ Oil ]) { + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { 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 ] ]; } @@ -999,13 +1008,13 @@ namespace Opm const Opm::PhaseUsage& pu = phaseUsage(); std::vector rates(3, 0.0); - if (active()[ Water ]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { 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 ] ]; } - if (active()[ Gas ]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { 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_); 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 - if (pu.phase_used[Gas] && pu.phase_used[Oil]) { + if (oilPresent && gasPresent) { rsmax_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 temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); - if (pu.phase_used[Water]) { - b_perf[ pu.phase_pos[Water] + perf * num_components_] = + if (waterPresent) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + b_perf[ waterCompIdx + perf * num_components_] = FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); } - if (pu.phase_used[Gas]) { - const int gaspos = pu.phase_pos[Gas] + perf * num_components_; + if (gasPresent) { + 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; - if (pu.phase_used[Oil]) { + if (oilPresent) { 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 rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); @@ -1146,10 +1161,11 @@ namespace Opm } } - if (pu.phase_used[Oil]) { - const int oilpos = pu.phase_pos[Oil] + perf * num_components_; + if (oilPresent) { + 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; - if (pu.phase_used[Gas]) { + if (gasPresent) { rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); 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); @@ -1170,8 +1186,13 @@ namespace Opm } // Surface density. - for (int p = 0; p < pu.num_phases; ++p) { - surf_dens_perf[num_components_ * perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + 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 @@ -1199,7 +1220,6 @@ namespace Opm const int np = number_of_phases_; const int nperf = number_of_perforations_; const int num_comp = num_components_; - const PhaseUsage& phase_usage = phaseUsage(); // 1. Compute the flow (in surface volume units for each // component) exiting up the wellbore from each perforation, @@ -1227,12 +1247,9 @@ namespace Opm // absolute values of the surface rates divided by their sum. // Then compute volume ratios (formation factors) for 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 mix(num_comp,0.0); std::vector x(num_comp); std::vector surf_dens(num_comp); - std::vector dens(nperf); for (int perf = 0; perf < nperf; ++perf) { // Find component mix. @@ -1244,28 +1261,36 @@ namespace Opm } } else { // No flow => use well specified fractions for mix. - for (int phase = 0; phase < np; ++phase) { - mix[phase] = comp_frac_[phase]; + for (int component = 0; component < num_comp; ++component) { + if (component < np) { + mix[component] = comp_frac_[ ebosCompIdxToFlowCompIdx(component)]; + } } // intialize 0.0 for comIdx >= np; } // Compute volume ratio. x = mix; - double rs = 0.0; - double rv = 0.0; - if (!rsmax_perf.empty() && mix[oilpos] > 0.0) { - rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]); - } - if (!rvmax_perf.empty() && mix[gaspos] > 0.0) { - rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); - } - if (rs != 0.0) { - // 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);; + + // Subtract dissolved gas from oil phase and vapporized oil from gas phase + if (FluidSystem::phaseIsActive(FluidSystem::gasCompIdx) && FluidSystem::phaseIsActive(FluidSystem::oilCompIdx)) { + const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const unsigned oilpos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + double rs = 0.0; + double rv = 0.0; + if (!rsmax_perf.empty() && mix[oilpos] > 0.0) { + rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]); + } + if (!rvmax_perf.empty() && mix[gaspos] > 0.0) { + rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); + } + if (rs != 0.0) { + // 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; for (int component = 0; component < num_comp; ++component) { @@ -1329,8 +1354,6 @@ namespace Opm StandardWell:: getWellConvergence(const std::vector& B_avg) const { - const int np = number_of_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 assert((int(B_avg.size()) == num_components_) || has_polymer); @@ -1356,18 +1379,23 @@ namespace Opm ConvergenceReport report; // checking if any NaN or too large residuals found - // TODO: not understand why phase here while component in other places. - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(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; - const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName}; + const typename ConvergenceReport::ProblemWell problem_well = {name(), compName}; report.nan_residual_wells.push_back(problem_well); } else { - if (well_flux_residual[phaseIdx] > maxResidualAllowed) { + if (well_flux_residual[compIdx] > maxResidualAllowed) { 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); } } @@ -1405,8 +1433,8 @@ namespace Opm std::vector perfRates(b_perf.size(),0.0); for (int perf = 0; perf < nperf; ++perf) { - for (int phase = 0; phase < np; ++phase) { - perfRates[perf * num_components_ + phase] = well_state.perfPhaseRates()[(first_perf_ + perf) * np + phase]; + for (int comp = 0; comp < np; ++comp) { + perfRates[perf * num_components_ + comp] = well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(comp)]; } if(has_solvent) { 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); 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(); std::vector rates(3, 0.0); - if (active()[ Water ]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { rates[ Water ] = potentials[pu.phase_pos[ Water ] ]; } - if (active()[ Oil ]) { + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { rates[ Oil ] = potentials[pu.phase_pos[ Oil ] ]; } - if (active()[ Gas ]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { 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]; } 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; } - 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 ; } if (has_solvent) { @@ -1798,7 +1826,7 @@ namespace Opm } else { // tot_well_rate == 0 if (well_type_ == INJECTOR) { // only single phase injection handled - if (active()[Water]) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { if (distr[Water] > 0.0) { primary_variables_[WFrac] = 1.0; } else { @@ -1806,7 +1834,7 @@ namespace Opm } } - if (active()[Gas]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (distr[pu.phase_pos[Gas]] > 0.0) { primary_variables_[GFrac] = 1.0 - wsolvent(); if (has_solvent) { @@ -1822,10 +1850,10 @@ namespace Opm // this will happen. } else if (well_type_ == PRODUCER) { // producers // 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; } - if (active()[Gas]) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { primary_variables_[GFrac] = 1.0 / np; } } else { @@ -1953,7 +1981,8 @@ namespace Opm if (well_type_ == INJECTOR) { // assume fully mixing within injecting wellbore 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()) { @@ -1974,10 +2003,11 @@ namespace Opm material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx); const double swcr = scaled_drainage_info.Swcr; 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 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()) { // 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(), water_velocity); // modify the mobility with the shear factor. - mob[Water] /= shear_factor; + mob[waterCompIdx] /= shear_factor; } } } diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index 979fdf115..28fe69afa 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -74,11 +74,11 @@ namespace Opm typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; 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, MaterialLaw) MaterialLaw; - static const int numEq = BlackoilIndices::numEq; + static const int numEq = Indices::numEq; typedef double Scalar; typedef Dune::FieldVector VectorBlockType; @@ -91,8 +91,8 @@ namespace Opm static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); - static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; - static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; + static const int contiSolventEqIdx = Indices::contiSolventEqIdx; + static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx; // For the conversion between the surface volume rate and resrevoir voidage rate using RateConverterType = RateConverter:: @@ -123,7 +123,6 @@ namespace Opm void setVFPProperties(const VFPProperties* vfp_properties_arg); virtual void init(const PhaseUsage* phase_usage_arg, - const std::vector* active_arg, const std::vector& depth_arg, const double gravity_arg, const int num_cells); @@ -269,8 +268,6 @@ namespace Opm bool getAllowCrossFlow() const; - const std::vector* active_; - const VFPProperties* vfp_properties_; double gravity_; @@ -284,13 +281,11 @@ namespace Opm const int num_components_; - const std::vector& active() const; - const PhaseUsage& phaseUsage() const; int flowPhaseToEbosCompIdx( const int phaseIdx ) const; - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; + int ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const; double wsolvent() const; diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index bcc62dc22..2abe15155 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -113,13 +113,11 @@ namespace Opm void WellInterface:: init(const PhaseUsage* phase_usage_arg, - const std::vector* active_arg, const std::vector& /* depth_arg */, const double gravity_arg, const int /* num_cells */) { phase_usage_ = phase_usage_arg; - active_ = active_arg; gravity_ = gravity_arg; } @@ -186,21 +184,6 @@ namespace Opm - - template - const std::vector& - WellInterface:: - active() const - { - assert(active_); - - return *active_; - } - - - - - template void WellInterface:: @@ -233,39 +216,32 @@ namespace Opm flowPhaseToEbosCompIdx( const int phaseIdx ) const { const auto& pu = phaseUsage(); - if (active()[Water] && pu.phase_pos[Water] == phaseIdx) - return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) - return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) - return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && pu.phase_pos[Water] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && pu.phase_pos[Oil] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && pu.phase_pos[Gas] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); // for other phases return the index return phaseIdx; } - - - template int WellInterface:: - flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + ebosCompIdxToFlowCompIdx( const unsigned compIdx ) 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; - } + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == compIdx) + return pu.phase_pos[Water]; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == compIdx) + return pu.phase_pos[Oil]; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == compIdx) + return pu.phase_pos[Gas]; - assert(phaseIdx < 3); // for other phases return the index - return phaseIdx; + return compIdx; } @@ -461,7 +437,7 @@ namespace Opm const int np = number_of_phases_; 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 min_oil_rate = econ_production_limits.minOilRate(); if (std::abs(oil_rate) < min_oil_rate) { @@ -470,7 +446,7 @@ namespace Opm } 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 min_gas_rate = econ_production_limits.minGasRate(); if (std::abs(gas_rate) < min_gas_rate) { @@ -479,8 +455,8 @@ namespace Opm } if (econ_production_limits.onMinLiquidRate() ) { - assert(active()[Oil]); - assert(active()[Water]); + assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)); + assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)); 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 liquid_rate = oil_rate + water_rate; @@ -517,8 +493,8 @@ namespace Opm const Opm::PhaseUsage& pu = phaseUsage(); const int well_number = index_of_well_; - assert(active()[Oil]); - assert(active()[Water]); + assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)); + assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)); 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 ] ]; @@ -854,11 +830,11 @@ namespace Opm return distr[phaseIdx]; } 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; - if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && pu.phase_pos[Oil] == phaseIdx) return 1.0; - if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && pu.phase_pos[Gas] == phaseIdx) return 0.01; if (has_solvent && phaseIdx == contiSolventEqIdx ) return 0.01;