diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 2392fa88d..fcabf8bb6 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -48,7 +48,6 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/VFPInjProperties.cpp opm/autodiff/WellMultiSegment.cpp opm/autodiff/MultisegmentWells.cpp - opm/autodiff/BlackoilSolventState.cpp opm/autodiff/MissingFeatures.cpp opm/polymer/PolymerState.cpp opm/polymer/PolymerBlackoilState.cpp @@ -169,7 +168,6 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/BlackoilSequentialModel.hpp opm/autodiff/BlackoilSolventModel.hpp opm/autodiff/BlackoilSolventModel_impl.hpp - opm/autodiff/BlackoilSolventState.hpp opm/autodiff/BlackoilMultiSegmentModel.hpp opm/autodiff/BlackoilMultiSegmentModel_impl.hpp opm/autodiff/BlackoilTransportModel.hpp diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 81a1dc5ab..797d67f81 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -84,6 +84,7 @@ NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem)); SET_BOOL_PROP(EclFlowProblem, DisableWells, true); SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false); SET_BOOL_PROP(EclFlowProblem, ExportGlobalTransmissibility, true); +SET_BOOL_PROP(EclFlowProblem, EnableSolvent, false); // SWATINIT is done by the flow part of flow_ebos. this can be removed once the legacy // code for fluid and satfunc handling gets fully retired. @@ -119,6 +120,7 @@ namespace Opm { typedef double Scalar; static const int numEq = BlackoilIndices::numEq; + static const int solventCompIdx = 3; //TODO get this from ebos typedef Dune::FieldVector VectorBlockType; typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; @@ -162,6 +164,7 @@ namespace Opm { , active_(detail::activePhases(fluid.phaseUsage())) , has_disgas_(FluidSystem::enableDissolvedGas()) , has_vapoil_(FluidSystem::enableVaporizedOil()) + , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) , param_( param ) , well_model_ (well_model) , terminal_output_ (terminal_output) @@ -624,8 +627,14 @@ namespace Opm { } dso -= dsg; + // solvent + const double dss = has_solvent_ ? dx[cell_idx][BlackoilIndices::solventSaturationIdx] : 0.0; + dso -= dss; + // Appleyard chop process. maxVal = std::max(std::abs(dsg),maxVal); + maxVal = std::max(std::abs(dss),maxVal); + double step = dsMax()/maxVal; step = std::min(step, 1.0); @@ -638,6 +647,12 @@ namespace Opm { double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; sg -= step * dsg; } + + if (has_solvent_) { + double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx]; + ss -= step * dss; + } + double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]]; so -= step * dso; @@ -682,12 +697,20 @@ namespace Opm { if (sg <= 0.0 && has_disgas_) { reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs sg = 0; - so = 1.0 - sw - sg; + so = 1.0 - sw; + if (has_solvent_) { + double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx]; + so -= ss; + } rs *= (1-epsilon); } else if (so <= 0.0 && has_vapoil_) { reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasOnly; // sg --> rv so = 0; - sg = 1.0 - sw - so; + sg = 1.0 - sw; + if (has_solvent_) { + double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx]; + sg -= ss; + } rv *= (1-epsilon); } break; @@ -862,6 +885,7 @@ namespace Opm { const double tol_wells = param_.tolerance_wells_; const int nc = Opm::AutoDiffGrid::numCells(grid_); + const int np = numPhases(); const int numComp = numComponents(); Vector R_sum(numComp); @@ -895,16 +919,24 @@ namespace Opm { const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& fs = intQuants.fluidState(); - for ( int compIdx = 0; compIdx < numComp; ++compIdx ) + for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx ) { - Vector& R2_idx = R2[ compIdx ]; - Vector& B_idx = B[ compIdx ]; - const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(compIdx); - const int ebosCompIdx = flowPhaseToEbosCompIdx(compIdx); + Vector& R2_idx = R2[ phaseIdx ]; + Vector& B_idx = B[ phaseIdx ]; + const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx); + const int ebosCompIdx = flowPhaseToEbosCompIdx(phaseIdx); B_idx [cell_idx] = 1.0 / fs.invB(ebosPhaseIdx).value(); R2_idx[cell_idx] = ebosResid[cell_idx][ebosCompIdx]; } + + if (has_solvent_ ) { + Vector& R2_idx = R2[ solventCompIdx ]; + Vector& B_idx = B[ solventCompIdx ]; + B_idx [cell_idx] = 1.0 / intQuants.solventInverseFormationVolumeFactor().value(); + R2_idx[cell_idx] = ebosResid[cell_idx][solventCompIdx]; + } + } Vector pv_vector; @@ -969,6 +1001,9 @@ namespace Opm { const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); key[ phaseIdx ] = std::toupper( phaseName.front() ); } + if (has_solvent_) { + key[ solventCompIdx ] = "S"; + } for (int compIdx = 0; compIdx < numComp; ++compIdx) { msg += " MB(" + key[ compIdx ] + ") "; @@ -1029,7 +1064,11 @@ namespace Opm { if (numPhases() == 2) { return 2; } - return FluidSystem::numComponents; + int numComp = FluidSystem::numComponents; + if (has_solvent_) + numComp ++; + + return numComp; } /// Wrapper required due to not following generic API @@ -1273,6 +1312,11 @@ namespace Opm { VectorType& pcSwMdc_ow = simData.getCellData( "PCSWMDC_OW" ); VectorType& krnSwMdc_ow = simData.getCellData( "KRNSWMDC_OW" ); + if (has_solvent_) { + simData.registerCellData( "SSOL", 1 ); + } + VectorType& ssol = has_solvent_ ? simData.getCellData( "SSOL" ) : zero; + std::vector failed_cells_pb; std::vector failed_cells_pd; const auto& gridView = ebosSimulator().gridView(); @@ -1356,6 +1400,11 @@ namespace Opm { krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value(); } + if (has_solvent_) + { + ssol[cellIdx] = intQuants.solventSaturation().value(); + } + // hack to make the intial output of rs and rv Ecl compatible. // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite @@ -1464,6 +1513,7 @@ namespace Opm { const std::vector cells_; // All grid cells const bool has_disgas_; const bool has_vapoil_; + const bool has_solvent_; ModelParameters param_; SimulatorReport failureReport_; @@ -1524,6 +1574,10 @@ namespace Opm { // set water saturation cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]]; + if (has_solvent_) { + cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx]; + } + // set switching variable and interpretation if (active_[Gas] ) { if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ ) @@ -1599,19 +1653,18 @@ namespace Opm { int flowToEbosPvIdx( const int flowPv ) const { - assert(flowPv < 3); - const int flowToEbos[ 3 ] = { + const int flowToEbos[ 4 ] = { BlackoilIndices::pressureSwitchIdx, BlackoilIndices::waterSaturationIdx, - BlackoilIndices::compositionSwitchIdx + BlackoilIndices::compositionSwitchIdx, + BlackoilIndices::solventSaturationIdx }; return flowToEbos[ flowPv ]; } int flowPhaseToEbosCompIdx( const int phaseIdx ) const { - assert(phaseIdx < 3); - const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; + const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx }; return phaseToComp[ phaseIdx ]; } @@ -1644,6 +1697,13 @@ namespace Opm { cellRes[ flowPhaseToEbosCompIdx( flowPhaseIdx ) ] /= refDens; cellRes[ flowPhaseToEbosCompIdx( flowPhaseIdx ) ] *= cellVolume; } + if (has_solvent_) { + // no need to store refDens for all cells? + const auto& intQuants = ebosSimulator_.model().cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0); + const auto& refDens = intQuants->solventRefDensity(); + cellRes[ solventCompIdx ] /= refDens; + cellRes[ solventCompIdx ] *= cellVolume; + } } for( auto row = ebosJac.begin(); row != endrow; ++row ) @@ -1669,13 +1729,24 @@ namespace Opm { (*col)[ebosCompIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume; } } + if (has_solvent_) { + // TODO store refDens pr pvtRegion? + const auto& intQuants = ebosSimulator_.model().cachedIntensiveQuantities(rowIdx, /*timeIdx=*/0); + const auto& refDens = intQuants->solventRefDensity(); + for( int pvIdx=0; pvIdx < numEq; ++pvIdx ) + { + (*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] /= refDens; + (*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume; + } + } } } } int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const { - const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; + assert(phaseIdx < 3); + const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}; return flowToEbos[ phaseIdx ]; } diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index 2c56f4fbe..87606e24a 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -22,7 +22,6 @@ #include #include -#include #include #include #include @@ -255,7 +254,7 @@ namespace Opm { template struct ModelTraits< BlackoilSolventModel > { - typedef BlackoilSolventState ReservoirState; + typedef BlackoilState ReservoirState; typedef WellStateFullyImplicitBlackoilSolvent WellState; typedef BlackoilModelParameters ModelParameters; typedef BlackoilSolventSolutionState SolutionState; diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index 807bb6fbf..422c67eda 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -135,7 +135,7 @@ namespace Opm { // Initial solvent concentration. if (has_solvent_) { - const auto& solvent_saturation = x.getCellData( BlackoilSolventState::SSOL ); + const auto& solvent_saturation = x.getCellData( BlackoilState::SSOL ); const int nc = solvent_saturation.size(); const V ss = Eigen::Map(solvent_saturation.data() , nc); diff --git a/opm/autodiff/BlackoilSolventState.cpp b/opm/autodiff/BlackoilSolventState.cpp deleted file mode 100644 index 90e8bccd7..000000000 --- a/opm/autodiff/BlackoilSolventState.cpp +++ /dev/null @@ -1,34 +0,0 @@ -/* - Copyright 2015 IRIS AS, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - - -#include - -namespace Opm -{ - const std::string BlackoilSolventState::SSOL = "SSOL"; - - BlackoilSolventState::BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases) - : BlackoilState( number_of_cells , number_of_faces , number_of_phases) - { - registerCellData( SSOL , 1 ); - } - -} // namespace Opm - diff --git a/opm/autodiff/BlackoilSolventState.hpp b/opm/autodiff/BlackoilSolventState.hpp deleted file mode 100644 index 0103e0904..000000000 --- a/opm/autodiff/BlackoilSolventState.hpp +++ /dev/null @@ -1,46 +0,0 @@ -/* - Copyright 2015 IRIS AS, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED -#define OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED - -#include - - #include - - -namespace Opm -{ - - /// Simulator state for blackoil simulator with solvent. - /// We use the Blackoil state parameters. - class BlackoilSolventState : public BlackoilState - { - public: - static const std::string SSOL; - - BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases); - }; - -} // namespace Opm - - - - -#endif // OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED diff --git a/opm/autodiff/Compat.cpp b/opm/autodiff/Compat.cpp index 6a5735394..7dbc4746d 100644 --- a/opm/autodiff/Compat.cpp +++ b/opm/autodiff/Compat.cpp @@ -29,7 +29,6 @@ #include #include #include -#include #include #include #include @@ -111,7 +110,7 @@ data::Solution simToSolution( const SimulationDataContainer& reservoir, } if (phases.has_solvent) { - sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilSolventState::SSOL ) , data::TargetType::RESTART_SOLUTION ); + sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilState::SSOL ) , data::TargetType::RESTART_SOLUTION ); } return sol; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 387a055b6..277d3cac9 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -252,7 +252,9 @@ public: // Run a multiple steps of the solver depending on the time step control. solver_timer.start(); - WellModel well_model(wells, &(wells_manager.wellCollection()), model_param_, terminal_output_); + const auto& wells_ecl = eclState().getSchedule().getWells(timer.currentStepNum()); + WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_, terminal_output_, timer.currentStepNum()); + auto solver = createSolver(well_model); std::vector> currentFluidInPlace; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp index 4bc7ee446..45c1ac6cc 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp @@ -22,7 +22,6 @@ #include #include -#include #include #include @@ -85,7 +84,7 @@ namespace Opm struct SimulatorTraits > { typedef WellStateFullyImplicitBlackoilSolvent WellState; - typedef BlackoilSolventState ReservoirState; + typedef BlackoilState ReservoirState; typedef BlackoilOutputWriter OutputWriter; typedef GridT Grid; typedef BlackoilSolventModel Model; diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 5c2906a49..8ba09e20e 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -62,7 +62,8 @@ namespace Opm { enum WellVariablePositions { XvarWell = 0, WFrac = 1, - GFrac = 2 + GFrac = 2, + SFrac = 3 }; @@ -85,6 +86,7 @@ enum WellVariablePositions { typedef double Scalar; static const int numEq = BlackoilIndices::numEq; static const int numWellEq = numEq; //number of wellEq is the same as numEq in the model + static const int solventCompIdx = 3; //TODO get this from ebos typedef Dune::FieldVector VectorBlockType; typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; @@ -102,8 +104,10 @@ enum WellVariablePositions { // --------- Public methods --------- StandardWellsDense(const Wells* wells_arg, WellCollection* well_collection, + const std::vector< const Well* >& wells_ecl, const ModelParameters& param, - const bool terminal_output); + const bool terminal_output, + const int current_index); void init(const PhaseUsage phase_usage_arg, const std::vector& active_arg, @@ -118,10 +122,15 @@ enum WellVariablePositions { /// The number of components in the model. int numComponents() const { - if (phase_usage_.num_phases == 2) { + if (numPhases() == 2) { return 2; } - return FluidSystem::numComponents; + int numComp = FluidSystem::numComponents; + if (has_solvent_) { + numComp ++; + } + + return numComp; } @@ -279,12 +288,15 @@ enum WellVariablePositions { protected: bool wells_active_; const Wells* wells_; + const std::vector< const Well* > wells_ecl_; // Well collection is used to enforce the group control WellCollection* well_collection_; ModelParameters param_; bool terminal_output_; + bool has_solvent_; + int current_timeIdx_; PhaseUsage phase_usage_; std::vector active_; @@ -385,6 +397,8 @@ enum WellVariablePositions { const int well_index, const double initial_bhp, // bhp from BHP constraints const std::vector& initial_potential) const; + + double wsolvent(const int well_index) const; }; diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index d1af3d5dd..11781c98b 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -7,13 +7,18 @@ namespace Opm { StandardWellsDense:: StandardWellsDense(const Wells* wells_arg, WellCollection* well_collection, + const std::vector< const Well* >& wells_ecl, const ModelParameters& param, - const bool terminal_output) + const bool terminal_output, + const int current_timeIdx) : wells_active_(wells_arg!=nullptr) , wells_(wells_arg) + , wells_ecl_(wells_ecl) , well_collection_(well_collection) , param_(param) , terminal_output_(terminal_output) + , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) + , current_timeIdx_(current_timeIdx) , well_perforation_efficiency_factors_((wells_!=nullptr ? wells_->well_connpos[wells_->number_of_wells] : 0), 1.0) , well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) , well_perforation_pressure_diffs_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) @@ -168,6 +173,7 @@ namespace Opm { { const int nw = wells().number_of_wells; const int numComp = numComponents(); + const int np = numPhases(); // clear all entries duneB_ = 0.0; @@ -223,7 +229,11 @@ namespace Opm { } // Store the perforation phase flux for later usage. - well_state.perfPhaseRates()[perf*numComp + componentIdx] = cq_s[componentIdx].value(); + if (componentIdx == solventCompIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent) + well_state.perfRateSolvent()[perf] = cq_s[componentIdx].value(); + } else { + well_state.perfPhaseRates()[perf*np + componentIdx] = cq_s[componentIdx].value(); + } } // Store the perforation pressure for later usage. @@ -267,6 +277,9 @@ namespace Opm { int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); } + if (has_solvent_) { + mob[solventCompIdx] = extendEval(intQuants.solventMobility()); + } } else { const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); @@ -281,6 +294,11 @@ namespace Opm { int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx)); } + + // this may not work if viscosity and relperms has been modified? + if (has_solvent_) { + OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent"); + } } } @@ -443,8 +461,7 @@ namespace Opm { StandardWellsDense:: flowPhaseToEbosCompIdx( const int phaseIdx ) const { - assert(phaseIdx < 3); - const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; + const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx }; return phaseToComp[ phaseIdx ]; } @@ -457,11 +474,11 @@ namespace Opm { StandardWellsDense:: flowToEbosPvIdx( const int flowPv ) const { - assert(flowPv < 3); - const int flowToEbos[ 3 ] = { + const int flowToEbos[ 4 ] = { BlackoilIndices::pressureSwitchIdx, BlackoilIndices::waterSaturationIdx, - BlackoilIndices::compositionSwitchIdx + BlackoilIndices::compositionSwitchIdx, + BlackoilIndices::solventSaturationIdx }; return flowToEbos[ flowPv ]; } @@ -754,6 +771,9 @@ namespace Opm { int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); } + if (has_solvent_) { + b_perfcells_dense[solventCompIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor()); + } // Pressure drawdown (also used to determine direction of flow) EvalWell well_pressure = bhp + cdp; @@ -803,6 +823,10 @@ namespace Opm { volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; } + if (has_solvent_) { + volumeRatio += cmix_s[solventCompIdx] / b_perfcells_dense[solventCompIdx]; + } + if (active_[Oil] && active_[Gas]) { const int oilpos = pu.phase_pos[Oil]; const int gaspos = pu.phase_pos[Gas]; @@ -989,6 +1013,10 @@ namespace Opm { B += 1 / fs.invB(ebosPhaseIdx).value(); } + if (has_solvent_) { + auto& B = B_avg[ solventCompIdx ]; + B += 1 / intQuants.solventInverseFormationVolumeFactor().value(); + } } // compute global average @@ -1129,14 +1157,14 @@ namespace Opm { if (pu.phase_used[BlackoilPhases::Vapour]) { const int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * numComp; - const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * numComp; + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; if (pu.phase_used[BlackoilPhases::Liquid]) { - const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * numComp; + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); if (oilrate > 0) { - const double gasrate = std::abs(xw.wellRates()[gaspos_well]); + const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); double rv = 0.0; if (gasrate > 0) { rv = oilrate / gasrate; @@ -1156,11 +1184,11 @@ namespace Opm { if (pu.phase_used[BlackoilPhases::Liquid]) { const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * numComp; - const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * numComp; + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; if (pu.phase_used[BlackoilPhases::Vapour]) { rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); - const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * numComp; - const double gasrate = std::abs(xw.wellRates()[gaspos_well]); + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); if (gasrate > 0) { const double oilrate = std::abs(xw.wellRates()[oilpos_well]); double rs = 0.0; @@ -1181,6 +1209,12 @@ namespace Opm { for (int p = 0; p < pu.num_phases; ++p) { surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); } + + // We use cell values for solvent injector + if (has_solvent_) { + b_perf[numComp*perf + solventCompIdx] = intQuants.solventInverseFormationVolumeFactor().value(); + surf_dens_perf[numComp*perf + solventCompIdx] = intQuants.solventRefDensity(); + } } } } @@ -1220,6 +1254,12 @@ namespace Opm { well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited; } + if (has_solvent_) { + const int sign4 = dwells[w][flowPhaseToEbosCompIdx(SFrac)] > 0 ? 1: -1; + const double dx4_limited = sign4 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(SFrac)]),dFLimit); + well_state.wellSolutions()[SFrac*nw + w] = xvar_well_old[SFrac*nw + w] - dx4_limited; + } + assert(active_[ Oil ]); F[Oil] = 1.0; if (active_[ Water ]) { @@ -1232,11 +1272,20 @@ namespace Opm { F[Oil] -= F[Gas]; } + double F_solvent = 0.0; + if (has_solvent_) { + F_solvent = well_state.wellSolutions()[SFrac*nw + w]; + F[Oil] -= F_solvent; + } + if (active_[ Water ]) { if (F[Water] < 0.0) { if (active_[ Gas ]) { F[Gas] /= (1.0 - F[Water]); } + if (has_solvent_) { + F_solvent /= (1.0 - F[Water]); + } F[Oil] /= (1.0 - F[Water]); F[Water] = 0.0; } @@ -1246,6 +1295,9 @@ namespace Opm { if (active_[ Water ]) { F[Water] /= (1.0 - F[Gas]); } + if (has_solvent_) { + F_solvent /= (1.0 - F[Gas]); + } F[Oil] /= (1.0 - F[Gas]); F[Gas] = 0.0; } @@ -1257,6 +1309,9 @@ namespace Opm { if (active_[ Gas ]) { F[Gas] /= (1.0 - F[Oil]); } + if (has_solvent_) { + F_solvent /= (1.0 - F[Oil]); + } F[Oil] = 0.0; } @@ -1266,6 +1321,15 @@ namespace Opm { if (active_[ Gas ]) { well_state.wellSolutions()[GFrac*nw + w] = F[Gas]; } + if(has_solvent_) { + well_state.wellSolutions()[SFrac*nw + w] = F_solvent; + } + + // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent. + // More testing is needed to make sure this is correct for well groups and THP. + if (has_solvent_){ + F[Gas] += F_solvent; + } // The interpretation of the first well variable depends on the well control const WellControls* wc = wells().ctrls[w]; @@ -1699,9 +1763,21 @@ namespace Opm { const double grav) { // Compute densities + const int nperf = depth_perf.size(); + const int numComponent = b_perf.size() / nperf; + const int np = wells().number_of_phases; + std::vector perfRates(b_perf.size(),0.0); + for (int perf = 0; perf < nperf; ++perf) { + for (int phase = 0; phase < np; ++phase) { + perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[perf*np + phase]; + } + if(has_solvent_) { + perfRates[perf*numComponent + solventCompIdx] = xw.perfRateSolvent()[perf]; + } + } well_perforation_densities_ = WellDensitySegmented::computeConnectionDensities( - wells(), phase_usage_, xw.perfPhaseRates(), + wells(), phase_usage_, perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); // Compute pressure deltas @@ -2065,13 +2141,34 @@ namespace Opm { EvalWell qs = 0.0; const WellControls* wc = wells().ctrls[wellIdx]; const int np = wells().number_of_phases; - assert(compIdx < np); + assert(compIdx < numComponents()); const int nw = wells().number_of_wells; + const auto pu = phase_usage_; const double target_rate = well_controls_get_current_target(wc); // TODO: the formulation for the injectors decides it only work with single phase // surface rate injection control. Improvement will be required. if (wells().type[wellIdx] == INJECTOR) { + if (has_solvent_ ) { + double comp_frac = 0.0; + if (compIdx == solventCompIdx) { // solvent + comp_frac = wells().comp_frac[np*wellIdx + pu.phase_pos[ Gas ]] * wsolvent(wellIdx); + } else if (compIdx == pu.phase_pos[ Gas ]) { + comp_frac = wells().comp_frac[np*wellIdx + compIdx] * (1.0 - wsolvent(wellIdx)); + } else { + comp_frac = wells().comp_frac[np*wellIdx + compIdx]; + } + if (comp_frac == 0.0) { + return qs; //zero + } + + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { + return comp_frac * wellVariables_[nw*XvarWell + wellIdx]; + } + + qs.setValue(comp_frac * target_rate); + return qs; + } const double comp_frac = wells().comp_frac[np*wellIdx + compIdx]; if (comp_frac == 0.0) { return qs; @@ -2117,17 +2214,27 @@ namespace Opm { assert(phase_under_control >= 0); + EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(wellIdx, phase_under_control); + if (has_solvent_ && phase_under_control == Gas) { + // for GRAT controlled wells solvent is included in the target + wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(wellIdx, solventCompIdx); + } + if (compIdx == phase_under_control) { + if (has_solvent_ && phase_under_control == Gas) { + qs.setValue(target_rate * wellVolumeFractionScaled(wellIdx, Gas).value() / wellVolumeFractionScaledPhaseUnderControl.value() ); + return qs; + } qs.setValue(target_rate); return qs; } // TODO: not sure why the single phase under control will have near zero fraction const double eps = 1e-6; - if (wellVolumeFractionScaled(wellIdx, phase_under_control) < eps) { + if (wellVolumeFractionScaledPhaseUnderControl < eps) { return qs; } - return (target_rate * wellVolumeFractionScaled(wellIdx, compIdx) / wellVolumeFractionScaled(wellIdx, phase_under_control)); + return (target_rate * wellVolumeFractionScaled(wellIdx, compIdx) / wellVolumeFractionScaledPhaseUnderControl); } // when it is a combined two phase rate limit, such like LRAT @@ -2175,6 +2282,10 @@ namespace Opm { return wellVariables_[GFrac * nw + wellIdx]; } + if (compIdx == solventCompIdx) { + return wellVariables_[SFrac * nw + wellIdx]; + } + // Oil fraction EvalWell well_fraction = 1.0; if (active_[Water]) { @@ -2184,6 +2295,9 @@ namespace Opm { if (active_[Gas]) { well_fraction -= wellVariables_[GFrac * nw + wellIdx]; } + if (has_solvent_) { + well_fraction -= wellVariables_[SFrac * nw + wellIdx]; + } return well_fraction; } @@ -2198,6 +2312,10 @@ namespace Opm { { const WellControls* wc = wells().ctrls[wellIdx]; if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { + + if (has_solvent_ && compIdx == solventCompIdx) { + return wellVolumeFraction(wellIdx, compIdx); + } const double* distr = well_controls_get_current_distr(wc); assert(compIdx < 3); if (distr[compIdx] > 0.) { @@ -2208,8 +2326,7 @@ namespace Opm { return wellVolumeFraction(wellIdx, compIdx); } } - assert(compIdx < 3); - std::vector g = {1,1,0.01}; + std::vector g = {1,1,0.01,0.01}; return (wellVolumeFraction(wellIdx, compIdx) / g[compIdx]); } @@ -2601,7 +2718,10 @@ namespace Opm { xw.wellSolutions()[WFrac*nw + well_index] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate; } if (active_[ Gas ]) { - xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * xw.wellRates()[np*well_index + Gas] / tot_well_rate ; + xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * (1.0 - wsolvent(well_index)) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ; + } + if (has_solvent_) { + xw.wellSolutions()[SFrac*nw + well_index] = g[Gas] * wsolvent(well_index) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ; } } else { const WellType& well_type = wells().type[well_index]; @@ -2617,7 +2737,10 @@ namespace Opm { if (active_[Gas]) { if (distr[Gas] > 0.0) { - xw.wellSolutions()[GFrac * nw + well_index] = 1.0; + xw.wellSolutions()[GFrac * nw + well_index] = 1.0 - wsolvent(well_index); + if (has_solvent_) { + xw.wellSolutions()[SFrac * nw + well_index] = wsolvent(well_index); + } } else { xw.wellSolutions()[GFrac * nw + well_index] = 0.0; } @@ -2672,6 +2795,7 @@ namespace Opm { std::vector& well_flux) const { const int np = wells().number_of_phases; + const int numComp = numComponents(); well_flux.resize(np, 0.0); const bool allow_cf = allow_cross_flow(well_index, ebosSimulator); @@ -2679,8 +2803,8 @@ namespace Opm { const int cell_index = wells().well_cells[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_index, /*timeIdx=*/ 0)); // flux for each perforation - std::vector cq_s(np, 0.0); - std::vector mob(np, 0.0); + std::vector cq_s(numComp, 0.0); + std::vector mob(numComp, 0.0); getMobility(ebosSimulator, perf, cell_index, mob); computeWellFlux(well_index, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); @@ -2879,4 +3003,38 @@ namespace Opm { return potentials; } + template + double + StandardWellsDense:: + wsolvent(const int well_index) const { + + if (!has_solvent_) { + return 0.0; + } + + // loop over all wells until we find the well with the matching name + for (const auto& well : wells_ecl_) { + if (well->getStatus( current_timeIdx_ ) == WellCommon::SHUT) { + continue; + } + + WellInjectionProperties injection = well->getInjectionProperties(current_timeIdx_); + if (injection.injectorType == WellInjector::GAS) { + + double solventFraction = well->getSolventFraction(current_timeIdx_); + + // Look until we find the correct well + if (well->name() == wells().name[well_index]) { + return solventFraction; + } + } + } + // we didn't find it return 0; + assert(false); + return 0.0; + } + + + + } // namespace Opm diff --git a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp index de4b945f2..840c7ed33 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp @@ -67,6 +67,48 @@ namespace Opm // call init on base class BaseType :: init(wells, state, prevState); + + const int nw = wells->number_of_wells; + if (nw == 0) { + return; + } + const int nperf = wells->well_connpos[nw]; + perfRateSolvent_.clear(); + perfRateSolvent_.resize(nperf, 0.0); + + if (pu.has_solvent) { + + // intialize wells that have been there before + // order may change so the mapping is based on the well name + if( ! prevState.wellMap().empty() ) + { + typedef typename WellMapType :: const_iterator const_iterator; + const_iterator end = prevState.wellMap().end(); + for (int w = 0; w < nw; ++w) { + std::string name( wells->name[ w ] ); + const_iterator it = prevState.wellMap().find( name ); + if( it != end ) + { + const int newIndex = w; + + // perfSolventRates + int oldPerf_idx = (*it).second[ 1 ]; + const int num_perf_old_well = (*it).second[ 2 ]; + const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex]; + if( num_perf_old_well == num_perf_this_well ) + { + for (int perf = wells->well_connpos[ newIndex ]; + perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx ) + { + perfRateSolvent()[ perf ] = prevState.perfRateSolvent()[ oldPerf_idx ]; + } + } + } + } + } + } + + // TODO: the reason to keep this is to avoid getting defaulted value BHP // some facilities needed from opm-parser or opm-core // It is a little tricky, since sometimes before applying group control, the only @@ -88,9 +130,12 @@ namespace Opm if (nw == 0) { return; } + + const int np = wells_->number_of_phases; + const int numComp = pu.has_solvent? np+1:np; well_solutions_.clear(); - well_solutions_.resize(nw * np, 0.0); + well_solutions_.resize(nw * numComp, 0.0); std::vector g = {1.0,1.0,0.01}; for (int w = 0; w < nw; ++w) { WellControls* wc = wells_->ctrls[w]; @@ -136,8 +181,12 @@ namespace Opm wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + waterpos] / total_rates; } if( pu.phase_used[Gas] ) { - wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + gaspos] / total_rates ; + wellSolutions()[2*nw + w] = g[Gas] * (wellRates()[np*w + gaspos] - solventWellRate(w)) / total_rates ; } + if( pu.has_solvent) { + wellSolutions()[3*nw + w] = g[Gas] * solventWellRate(w) / total_rates; + } + } else { @@ -147,6 +196,10 @@ namespace Opm if( pu.phase_used[Gas] ) { wellSolutions()[2*nw + w] = wells_->comp_frac[np*w + gaspos]; } + if (pu.has_solvent) { + wellSolutions()[3*nw + w] = 0; + } + } } } @@ -163,14 +216,44 @@ namespace Opm std::vector& wellSolutions() { return well_solutions_; } const std::vector& wellSolutions() const { return well_solutions_; } + /// One rate pr well connection. + std::vector& perfRateSolvent() { return perfRateSolvent_; } + const std::vector& perfRateSolvent() const { return perfRateSolvent_; } + + /// One rate pr well + const double solventWellRate(const int w) const { + double solvent_well_rate = 0.0; + for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) { + solvent_well_rate += perfRateSolvent_[perf]; + } + return solvent_well_rate; + } + + data::Wells report(const PhaseUsage& pu) const override { data::Wells res = BaseType::report(pu); + const int nw = WellState::numWells(); + // If there are now wells numPhases throws a floating point + // exception. + if (nw == 0) { + return res; + } + if (pu.has_solvent) { + // add solvent component + for( int w = 0; w < nw; ++w ) { + using rt = data::Rates::opt; + res.at( wells_->name[ w ]).rates.set( rt::solvent, solventWellRate(w) ); + } + } + return res; } private: std::vector well_solutions_; + std::vector perfRateSolvent_; + }; } // namespace Opm