diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index f26132841..caf9ab29b 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -163,7 +163,7 @@ namespace Opm { const ModelParameters& param, const BlackoilPropsAdFromDeck& fluid, const DerivedGeology& geo , - const StandardWellsDense& well_model, + const StandardWellsDense& well_model, const NewtonIterationBlackoilInterface& linsolver, const bool terminal_output) : ebosSimulator_(ebosSimulator) @@ -1428,7 +1428,7 @@ namespace Opm { ModelParameters param_; // Well Model - StandardWellsDense well_model_; + StandardWellsDense well_model_; /// \brief Whether we print something to std::cout bool terminal_output_; @@ -1446,9 +1446,9 @@ namespace Opm { public: /// return the StandardWells object - StandardWellsDense& + StandardWellsDense& wellModel() { return well_model_; } - const StandardWellsDense& + const StandardWellsDense& wellModel() const { return well_model_; } /// return the Well struct in the StandardWells diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index a1c225009..733d90c67 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -57,6 +57,8 @@ public: 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, MaterialLaw) MaterialLaw; + typedef WellStateFullyImplicitBlackoilDense WellState; typedef BlackoilState ReservoirState; @@ -64,7 +66,7 @@ public: typedef BlackoilModelEbos Model; typedef BlackoilModelParameters ModelParameters; typedef NonlinearSolver Solver; - typedef StandardWellsDense WellModel; + typedef StandardWellsDense WellModel; /// Initialise from parameters and objects to observe. diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 3c13344e2..c4a185ca2 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -67,7 +67,7 @@ enum WellVariablePositions { /// Class for handling the standard well model. - template + template class StandardWellsDense { public: // --------- Types --------- @@ -114,6 +114,13 @@ enum WellVariablePositions { WellState& well_state, bool only_wells); + template + void + getMobility(const Simulator& ebosSimulator, + const int perf, + const int cell_idx, + std::vector& mob) const; + template bool allow_cross_flow(const int w, Simulator& ebosSimulator) const; @@ -180,9 +187,9 @@ enum WellVariablePositions { void computeAccumWells(); - template + template void - computeWellFlux(const int& w, const double& Tw, const intensiveQuants& intQuants, + computeWellFlux(const int& w, const double& Tw, const FluidState& fs, const std::vector& mob_perfcells_dense, const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector& cq_s) const; template diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 7d1d811c7..92739cd19 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -3,8 +3,8 @@ namespace Opm { - template - StandardWellsDense:: + template + StandardWellsDense:: StandardWellsDense(const Wells* wells_arg, WellCollection* well_collection, const ModelParameters& param, @@ -32,9 +32,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: init(const PhaseUsage phase_usage_arg, const std::vector& active_arg, const VFPProperties* vfp_properties_arg, @@ -120,10 +120,10 @@ namespace Opm { - template + template template SimulatorReport - StandardWellsDense:: + StandardWellsDense:: assemble(Simulator& ebosSimulator, const int iterationIdx, const double dt, @@ -162,10 +162,10 @@ namespace Opm { - template + template template void - StandardWellsDense:: + StandardWellsDense:: assembleWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state, @@ -192,7 +192,10 @@ namespace Opm { const int cell_idx = wells().well_cells[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); std::vector cq_s(np,0.0); - computeWellFlux(w, wells().WI[perf], intQuants, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); + + std::vector mob(np, 0.0); + getMobility(ebosSimulator, perf, cell_idx, mob); + computeWellFlux(w, wells().WI[perf], intQuants.fluidState(), mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); for (int p1 = 0; p1 < np; ++p1) { @@ -249,13 +252,51 @@ namespace Opm { } + template + template + void + StandardWellsDense:: + getMobility(const Simulator& ebosSimulator, const int perf, const int cell_idx, std::vector& mob) const + { + + const int np = wells().number_of_phases; + assert (mob.size() == np); + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); + + // either use mobility of the perforation cell or calcualte its own + // based on passing the saturation table index + const int satid = wells().sat_table_id[perf] - 1; + 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)); + } + } else { + + const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); + Eval relativePerms[3]; + MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); + + // reset the satnumvalue back to original + 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)); + } + } + } - template + template template bool - StandardWellsDense:: + StandardWellsDense:: allow_cross_flow(const int w, Simulator& ebosSimulator) const { if (wells().allow_cf[w]) { @@ -290,9 +331,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: localInvert(Mat& istlA) const { for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { @@ -307,9 +348,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: print(Mat& istlA) const { for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { @@ -323,9 +364,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: apply( BVector& r) const { if ( ! localWellsActive() ) { @@ -342,9 +383,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: apply(const BVector& x, BVector& Ax) { if ( ! localWellsActive() ) { @@ -365,9 +406,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) { if ( ! localWellsActive() ) { @@ -387,9 +428,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: recoverVariable(const BVector& x, BVector& xw) const { if ( ! localWellsActive() ) { @@ -404,9 +445,9 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: flowPhaseToEbosCompIdx( const int phaseIdx ) const { const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; @@ -417,9 +458,9 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: flowToEbosPvIdx( const int flowPv ) const { const int flowToEbos[ 3 ] = { @@ -434,9 +475,9 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: flowPhaseToEbosPhaseIdx( const int phaseIdx ) const { const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; @@ -447,9 +488,9 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: ebosCompToFlowPhaseIdx( const int compIdx ) const { const int compToPhase[ 3 ] = { Oil, Water, Gas }; @@ -460,9 +501,9 @@ namespace Opm { - template + template std::vector - StandardWellsDense:: + StandardWellsDense:: extractPerfData(const std::vector& in) const { const int nw = wells().number_of_wells; @@ -481,9 +522,9 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: numPhases() const { return wells().number_of_phases; @@ -493,9 +534,9 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: numCells() const { return pv_.size(); @@ -505,9 +546,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: resetWellControlFromState(WellState xw) const { const int nw = wells_->number_of_wells; @@ -521,9 +562,9 @@ namespace Opm { - template + template const Wells& - StandardWellsDense:: + StandardWellsDense:: wells() const { assert(wells_ != 0); @@ -534,9 +575,9 @@ namespace Opm { - template + template const Wells* - StandardWellsDense:: + StandardWellsDense:: wellsPointer() const { return wells_; @@ -546,9 +587,9 @@ namespace Opm { - template + template bool - StandardWellsDense:: + StandardWellsDense:: wellsActive() const { return wells_active_; @@ -558,9 +599,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: setWellsActive(const bool wells_active) { wells_active_ = wells_active; @@ -570,9 +611,9 @@ namespace Opm { - template + template bool - StandardWellsDense:: + StandardWellsDense:: localWellsActive() const { return wells_ ? (wells_->number_of_wells > 0 ) : false; @@ -582,9 +623,9 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: numWellVars() const { if ( !localWellsActive() ) { @@ -600,9 +641,9 @@ namespace Opm { - template + template const std::vector& - StandardWellsDense:: + StandardWellsDense:: wellPerforationDensities() const { return well_perforation_densities_; @@ -612,9 +653,9 @@ namespace Opm { - template + template const std::vector& - StandardWellsDense:: + StandardWellsDense:: wellPerforationPressureDiffs() const { return well_perforation_pressure_diffs_; @@ -624,9 +665,9 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: extendEval(Eval in) const { EvalWell out = 0.0; out.setValue(in.value()); @@ -640,9 +681,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: setWellVariables(const WellState& xw) { const int np = wells().number_of_phases; @@ -660,9 +701,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: print(EvalWell in) const { std::cout << in.value() << std::endl; @@ -675,9 +716,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: computeAccumWells() { const int np = wells().number_of_phases; @@ -693,12 +734,13 @@ namespace Opm { - template - template + template + template void - StandardWellsDense:: + StandardWellsDense:: computeWellFlux(const int& w, const double& Tw, - const intensiveQuants& intQuants, + const FluidState& fs, + const std::vector& mob_perfcells_dense, const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector& cq_s) const { @@ -710,16 +752,13 @@ namespace Opm { cmix_s[phase] = wellSurfaceVolumeFraction(w, phase); } - const auto& fs = intQuants.fluidState(); EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); EvalWell rs = extendEval(fs.Rs()); EvalWell rv = extendEval(fs.Rv()); std::vector b_perfcells_dense(np, 0.0); - std::vector mob_perfcells_dense(np, 0.0); for (int phase = 0; phase < np; ++phase) { int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); - mob_perfcells_dense[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); } // Pressure drawdown (also used to determine direction of flow) @@ -823,10 +862,10 @@ namespace Opm { - template + template template SimulatorReport - StandardWellsDense:: + StandardWellsDense:: solveWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state) @@ -875,9 +914,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: printIf(const int c, const double x, const double y, const double eps, const std::string type) const { if (std::abs(x-y) > eps) { @@ -889,9 +928,9 @@ namespace Opm { - template + template std::vector - StandardWellsDense:: + StandardWellsDense:: residual() const { if( ! wellsActive() ) @@ -916,10 +955,10 @@ namespace Opm { - template + template template bool - StandardWellsDense:: + StandardWellsDense:: getWellConvergence(Simulator& ebosSimulator, const int iteration) const { @@ -1029,10 +1068,10 @@ namespace Opm { - template + template template void - StandardWellsDense:: + StandardWellsDense:: computeWellConnectionPressures(const Simulator& ebosSimulator, const WellState& xw) { @@ -1053,10 +1092,10 @@ namespace Opm { - template + template template void - StandardWellsDense:: + StandardWellsDense:: computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator, const WellState& xw, std::vector& b_perf, @@ -1156,9 +1195,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: updateWellState(const BVector& dwells, WellState& well_state) const { @@ -1428,9 +1467,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: updateWellControls(WellState& xw) const { if( !localWellsActive() ) return ; @@ -1524,9 +1563,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: updateListEconLimited(const Schedule& schedule, const int current_step, const Wells* wells_struct, @@ -1632,9 +1671,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: computeWellConnectionDensitesPressures(const WellState& xw, const std::vector& b_perf, const std::vector& rsmax_perf, @@ -1659,10 +1698,10 @@ namespace Opm { - template + template template void - StandardWellsDense:: + StandardWellsDense:: computeWellPotentials(const Simulator& ebosSimulator, WellState& well_state) const { @@ -1781,7 +1820,9 @@ namespace Opm { const int cell_index = wells().well_cells[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_index, /*timeIdx=*/ 0)); std::vector well_potentials(np, 0.0); - computeWellFlux(w, wells().WI[perf], intQuants, bhp, wellPerforationPressureDiffs()[perf], allow_cf, well_potentials); + std::vector mob(np, 0.0); + getMobility(ebosSimulator, perf, cell_index, mob); + computeWellFlux(w, wells().WI[perf], intQuants.fluidState(), mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, well_potentials); for(int p = 0; p < np; ++p) { well_state.wellPotentials()[perf * np + p] = well_potentials[p].value(); } @@ -1793,9 +1834,9 @@ namespace Opm { - template + template WellCollection* - StandardWellsDense:: + StandardWellsDense:: wellCollection() const { return well_collection_; @@ -1804,9 +1845,9 @@ namespace Opm { - template + template const std::vector& - StandardWellsDense:: + StandardWellsDense:: wellPerfEfficiencyFactors() const { return well_perforation_efficiency_factors_; @@ -1816,9 +1857,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: calculateEfficiencyFactors() { if ( !localWellsActive() ) { @@ -1844,9 +1885,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: computeWellVoidageRates(const WellState& well_state, std::vector& well_voidage_rates, std::vector& voidage_conversion_coeffs) const @@ -1904,9 +1945,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: applyVREPGroupControl(WellState& well_state) const { if ( wellCollection()->havingVREPGroups() ) { @@ -1929,9 +1970,9 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: getBhp(const int wellIdx) const { const WellControls* wc = wells().ctrls[wellIdx]; if (well_controls_get_current_type(wc) == BHP) { @@ -1985,9 +2026,9 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: getQs(const int wellIdx, const int phaseIdx) const { EvalWell qs = 0.0; @@ -2088,9 +2129,9 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: wellVolumeFraction(const int wellIdx, const int phaseIdx) const { const int nw = wells().number_of_wells; @@ -2118,9 +2159,9 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const { const WellControls* wc = wells().ctrls[wellIdx]; @@ -2136,9 +2177,9 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: wellSurfaceVolumeFraction(const int well_index, const int phase) const { EvalWell sum_volume_fraction_scaled = 0.; @@ -2156,9 +2197,9 @@ namespace Opm { - template + template bool - StandardWellsDense:: + StandardWellsDense:: checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, const WellState& well_state, const int well_number) const @@ -2207,9 +2248,9 @@ namespace Opm { - template - typename StandardWellsDense::RatioCheckTuple - StandardWellsDense:: + template + typename StandardWellsDense::RatioCheckTuple + StandardWellsDense:: checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, const WellState& well_state, const WellMapEntryType& map_entry) const @@ -2265,9 +2306,9 @@ namespace Opm { - template - typename StandardWellsDense::RatioCheckTuple - StandardWellsDense:: + template + typename StandardWellsDense::RatioCheckTuple + StandardWellsDense:: checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, const WellState& well_state, const WellMapEntryType& map_entry) const @@ -2345,9 +2386,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: updateWellStateWithTarget(const WellControls* wc, const int current, const int well_index, diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index 7b17d1e66..c9949e9f2 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -448,7 +448,7 @@ BOOST_AUTO_TEST_CASE(GetTable) std::shared_ptr wells(create_wells(nphases, nwells, nperfs), destroy_wells); const int cells[] = {5}; - add_well(INJECTOR, 100, 1, NULL, cells, NULL, NULL, true, wells.get()); + add_well(INJECTOR, 100, 1, NULL, cells, NULL, 0, NULL, true, wells.get()); //Create interpolation points double aqua_d = -0.15; @@ -785,7 +785,7 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs) std::stringstream ss; ss << "WELL_" << i; const bool ok = add_well(INJECTOR, 0.0, 1, NULL, &cells, - NULL, ss.str().c_str(), true, wells.get()); + NULL, 0, ss.str().c_str(), true, wells.get()); BOOST_REQUIRE(ok); } diff --git a/tests/test_welldensitysegmented.cpp b/tests/test_welldensitysegmented.cpp index 5dd928de6..470aa9e48 100644 --- a/tests/test_welldensitysegmented.cpp +++ b/tests/test_welldensitysegmented.cpp @@ -52,9 +52,9 @@ BOOST_AUTO_TEST_CASE(TestPressureDeltas) const bool allow_crossflow = true; std::shared_ptr wells(create_wells(np, 2, nperf), destroy_wells); BOOST_REQUIRE(wells); - int ok = add_well(INJECTOR, ref_depth, nperf/2, comp_frac_w, cells, WI, "INJ", allow_crossflow, wells.get()); + int ok = add_well(INJECTOR, ref_depth, nperf/2, comp_frac_w, cells, WI, 0, "INJ", allow_crossflow, wells.get()); BOOST_REQUIRE(ok); - ok = add_well(PRODUCER, ref_depth, nperf/2, comp_frac_o, cells, WI, "PROD", allow_crossflow, wells.get()); + ok = add_well(PRODUCER, ref_depth, nperf/2, comp_frac_o, cells, WI, 0, "PROD", allow_crossflow, wells.get()); BOOST_REQUIRE(ok); std::vector rates = { 1.0, 0.0, 0.0, 1.0, 0.0, 0.0,