diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 06c285e5f..63c826111 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -3,7 +3,7 @@ Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services Copyright 2014, 2015 Statoil ASA. Copyright 2015 NTNU - Copyright 2015 IRIS AS + Copyright 2015, 2016, 2017 IRIS AS This file is part of the Open Porous Media project (OPM). @@ -33,7 +33,6 @@ #include #include #include -#include #include #include @@ -126,10 +125,6 @@ namespace Opm { typedef ISTLSolver< MatrixBlockType, VectorBlockType, Indices::pressureSwitchIdx > ISTLSolverType; //typedef typename SolutionVector :: value_type PrimaryVariables ; - // For the conversion between the surface volume rate and resrevoir voidage rate - using RateConverterType = RateConverter:: - SurfaceToReservoirVoidage >; - typedef Opm::FIPData FIPDataType; // --------- Public methods --------- diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index 86566e1a5..25feab433 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -221,8 +221,11 @@ namespace Opm { const Well* well_ecl = wells_ecl_[index_well]; + const int well_cell_top = wells()->well_cells[wells()->well_connpos[w]]; + const int pvtreg = pvt_region_idx_[well_cell_top]; + if ( !well_ecl->isMultiSegment(time_step) || !param_.use_multisegment_well_) { - well_container.emplace_back(new StandardWell(well_ecl, time_step, wells(), param_) ); + well_container.emplace_back(new StandardWell(well_ecl, time_step, wells(), param_, *rateConverter_, pvtreg ) ); } else { well_container.emplace_back(new MultisegmentWell(well_ecl, time_step, wells(), param_) ); } diff --git a/opm/autodiff/RateConverter.hpp b/opm/autodiff/RateConverter.hpp index d27dce076..b04a7a12c 100644 --- a/opm/autodiff/RateConverter.hpp +++ b/opm/autodiff/RateConverter.hpp @@ -1,6 +1,7 @@ /* Copyright 2014, 2015 SINTEF ICT, Applied Mathematics. Copyright 2014, 2015 Statoil ASA. + Copyright 2017, IRIS This file is part of the Open Porous Media Project (OPM). @@ -638,6 +639,31 @@ namespace Opm { } } + + /** + * Compute coefficients for surface-to-reservoir voidage + * conversion for solvent. + * + * + * \param[in] r Fluid-in-place region of the well + * \param[in] pvtRegionIdx PVT region of the well + * + * + * \param[out] double Surface-to-reservoir conversion + * coefficients for solvent. + */ + template + void + calcCoeffSolvent(const RegionId r, const int pvtRegionIdx, double& coeff) const + { + const auto& ra = attr_.attributes(r); + const double p = ra.pressure; + const double T = ra.temperature; + const auto& solventPvt = SolventModule::solventPvt(); + const double bs = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p); + coeff = 1.0 / bs; + } + private: /** * Fluid property object. diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 63d86dd33..ac763aeb2 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -67,7 +67,6 @@ public: typedef BlackoilModelParameters ModelParameters; typedef NonlinearSolver Solver; typedef BlackoilWellModel WellModel; - typedef RateConverter::SurfaceToReservoirVoidage > RateConverterType; /// Initialise from parameters and objects to observe. diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index ecf6fb645..ae05d849c 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -26,6 +26,7 @@ #include #include +#include namespace Opm { @@ -104,6 +105,10 @@ namespace Opm typedef DenseAd::Evaluation EvalWell; + // For the conversion between the surface volume rate and resrevoir voidage rate + using RateConverterType = RateConverter:: + SurfaceToReservoirVoidage >; + // TODO: should these go to WellInterface? static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; @@ -111,7 +116,10 @@ namespace Opm static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx; - StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param); + StandardWell(const Well* well, const int time_step, const Wells* wells, + const ModelParameters& param, + const RateConverterType& rate_converter, + const int pvtRegionIdx); virtual void init(const PhaseUsage* phase_usage_arg, const std::vector* active_arg, @@ -220,6 +228,9 @@ namespace Opm // the saturations in the well bore under surface conditions at the beginning of the time step std::vector F0_; + const RateConverterType& rateConverter_; + int pvtRegionIdx_; + // TODO: this function should be moved to the base class. // while it faces chanllenges for MSWell later, since the calculation of bhp // based on THP is never implemented for MSWell yet. diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 983db1775..b998ec203 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -24,13 +24,18 @@ namespace Opm { template StandardWell:: - StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param) + StandardWell(const Well* well, const int time_step, const Wells* wells, + const ModelParameters& param, + const RateConverterType& rate_converter, + const int pvtRegionIdx) : Base(well, time_step, wells, param) , perf_densities_(number_of_perforations_) , perf_pressure_diffs_(number_of_perforations_) , primary_variables_(numWellEq, 0.0) , primary_variables_evaluation_(numWellEq) // the number of the primary variables , F0_(numWellEq) + , rateConverter_(rate_converter) + , pvtRegionIdx_(pvtRegionIdx) { duneB_.setBuildMode( OffDiagMatWell::row_wise ); duneC_.setBuildMode( OffDiagMatWell::row_wise ); @@ -861,13 +866,7 @@ namespace Opm primary_variables_[SFrac] = 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[pu.phase_pos[Gas]] += F_solvent; - } - - // The interpretation of the first well variable depends on the well control + // The interpretation of the first well variable depends on the well control const WellControls* wc = well_controls_; // TODO: we should only maintain one current control either from the well_state or from well_controls struct. @@ -884,6 +883,18 @@ namespace Opm } } + // 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){ + const double scal = scalingFactor(contiSolventEqIdx); + if (scal > 0) { + F_solvent /= scal ; + } else { + F_solvent = 0.; + } + F[pu.phase_pos[Gas]] += F_solvent; + } + switch (well_controls_iget_type(wc, current)) { case THP: // The BHP and THP both uses the total rate as first well variable. case BHP: @@ -1981,9 +1992,13 @@ namespace Opm const double* distr = well_controls_get_current_distr(wc); if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { - if (has_solvent && phaseIdx == contiSolventEqIdx ) - OPM_THROW(std::runtime_error, "RESERVOIR_RATE control in combination with solvent is not implemented"); - + if (has_solvent && phaseIdx == contiSolventEqIdx ) { + typedef Ewoms::BlackOilSolventModule SolventModule; + double coeff = 0; + rateConverter_.template calcCoeffSolvent(0, pvtRegionIdx_, coeff); + return coeff; + } + // TODO: use the rateConverter here as well. return distr[phaseIdx]; } const auto& pu = phaseUsage(); diff --git a/tests/test_wellmodel.cpp b/tests/test_wellmodel.cpp index d1f633b55..81af57adf 100644 --- a/tests/test_wellmodel.cpp +++ b/tests/test_wellmodel.cpp @@ -132,9 +132,22 @@ BOOST_AUTO_TEST_CASE(TestStandardWellInput) { BOOST_CHECK_EQUAL( wells_ecl.size(), 2); const Opm::Well* well = wells_ecl[1]; const Opm::BlackoilModelParameters param; - BOOST_CHECK_THROW( StandardWell( well, -1, wells, param), std::invalid_argument); - BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells, param), std::invalid_argument); - BOOST_CHECK_THROW( StandardWell( well, 4, nullptr, param), std::invalid_argument); + + // For the conversion between the surface volume rate and resrevoir voidage rate + typedef Opm::FluidSystems::BlackOil FluidSystem; + using RateConverterType = Opm::RateConverter:: + SurfaceToReservoirVoidage >; + // Compute reservoir volumes for RESV controls. + Opm::PhaseUsage phaseUsage; + std::unique_ptr rateConverter; + // Compute reservoir volumes for RESV controls. + rateConverter.reset(new RateConverterType (phaseUsage, + std::vector(10, 0))); + + const int pvtIdx = 0; + BOOST_CHECK_THROW( StandardWell( well, -1, wells, param, *rateConverter, pvtIdx), std::invalid_argument); + BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells, param , *rateConverter, pvtIdx), std::invalid_argument); + BOOST_CHECK_THROW( StandardWell( well, 4, nullptr, param , *rateConverter, pvtIdx), std::invalid_argument); } @@ -160,8 +173,20 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) { } // we should always be able to find the well in wells_ecl BOOST_CHECK(index_well != wells_ecl.size()); + // For the conversion between the surface volume rate and resrevoir voidage rate + typedef Opm::FluidSystems::BlackOil FluidSystem; + using RateConverterType = Opm::RateConverter:: + SurfaceToReservoirVoidage >; + // Compute reservoir volumes for RESV controls. + Opm::PhaseUsage phaseUsage; + std::unique_ptr rateConverter; + // Compute reservoir volumes for RESV controls. + rateConverter.reset(new RateConverterType (phaseUsage, + std::vector(10, 0))); - wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct, param) ); + const int pvtIdx = 0; + + wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct, param, *rateConverter, pvtIdx) ); } }