diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 2551bd7fb..8b388ac10 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -387,11 +387,19 @@ namespace Opm { void computeWellConnectionPressures(const SolutionState& state, const WellState& xw); + void computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf); + void assembleMassBalanceEq(const SolutionState& state); void - extractWellPerfProperties(std::vector& mob_perfcells, + extractWellPerfProperties(const SolutionState& state, + std::vector& mob_perfcells, std::vector& b_perfcells) const; bool diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index ed89c4481..c9703992b 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -795,18 +795,14 @@ namespace detail { } } - - template - void BlackoilModelBase::computeWellConnectionPressures(const SolutionState& state, - const WellState& xw) + void BlackoilModelBase::computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) { - if( ! localWellsActive() ) return ; - - using namespace Opm::AutoDiffGrid; - // 1. Compute properties required by computeConnectionPressureDelta(). - // Note that some of the complexity of this part is due to the function - // taking std::vector arguments, and not Eigen objects. const int nperf = wells().well_connpos[wells().number_of_wells]; const int nw = wells().number_of_wells; @@ -837,8 +833,6 @@ namespace detail { } const PhaseUsage& pu = fluid_.phaseUsage(); DataBlock b(nperf, pu.num_phases); - std::vector rsmax_perf(nperf, 0.0); - std::vector rvmax_perf(nperf, 0.0); if (pu.phase_used[BlackoilPhases::Aqua]) { const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; @@ -851,6 +845,8 @@ namespace detail { b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; const V rssat = fluidRsSat(avg_press, perf_so, well_cells); rsmax_perf.assign(rssat.data(), rssat.data() + nperf); + } else { + rsmax_perf.assign(nperf, 0.0); } if (pu.phase_used[BlackoilPhases::Vapour]) { const ADB perf_rv = subset(state.rv, well_cells); @@ -858,13 +854,12 @@ namespace detail { b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; const V rvsat = fluidRvSat(avg_press, perf_so, well_cells); rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); + } else { + rvmax_perf.assign(nperf, 0.0); } // b is row major, so can just copy data. - std::vector b_perf(b.data(), b.data() + nperf * pu.num_phases); - // Extract well connection depths. - const V depth = cellCentroidsZToEigen(grid_); - const V pdepth = subset(depth, well_cells); - std::vector perf_depth(pdepth.data(), pdepth.data() + nperf); + b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); + // Surface density. // The compute density segment wants the surface densities as @@ -873,17 +868,43 @@ namespace detail { for (int phase = 1; phase < pu.num_phases; ++phase) { rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); } - std::vector surf_dens_perf(rho.data(), rho.data() + nperf * pu.num_phases); + surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); + } - // Gravity - double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); + template + void BlackoilModelBase::computeWellConnectionPressures(const SolutionState& state, + const WellState& xw) + { + if( ! localWellsActive() ) return ; + + using namespace Opm::AutoDiffGrid; + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + std::vector b_perf; + std::vector rsmax_perf; + std::vector rvmax_perf; + std::vector surf_dens_perf; + asImpl().computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + // 2. Compute densities std::vector cd = WellDensitySegmented::computeConnectionDensities( wells(), xw, fluid_.phaseUsage(), b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector& well_cells = wops_.well_cells; + + // Extract well connection depths. + const V depth = cellCentroidsZToEigen(grid_); + const V pdepth = subset(depth, well_cells); + std::vector perf_depth(pdepth.data(), pdepth.data() + nperf); + + // Gravity + double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); + // 3. Compute pressure deltas std::vector cdp = WellDensitySegmented::computeConnectionPressureDelta( @@ -954,7 +975,7 @@ namespace detail { std::vector mob_perfcells; std::vector b_perfcells; - asImpl().extractWellPerfProperties(mob_perfcells, b_perfcells); + asImpl().extractWellPerfProperties(state, mob_perfcells, b_perfcells); if (param_.solve_welleq_initially_ && initial_assembly) { // solve the well equations as a pre-processing step asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state); @@ -1101,7 +1122,8 @@ namespace detail { template void - BlackoilModelBase::extractWellPerfProperties(std::vector& mob_perfcells, + BlackoilModelBase::extractWellPerfProperties(const SolutionState&, + std::vector& mob_perfcells, std::vector& b_perfcells) const { // If we have wells, extract the mobilities and b-factors for diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp index f58ac64c2..5a2c5fcb2 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -630,7 +630,7 @@ namespace Opm { std::vector mob_perfcells; std::vector b_perfcells; - asImpl().extractWellPerfProperties(mob_perfcells, b_perfcells); + asImpl().extractWellPerfProperties(state, mob_perfcells, b_perfcells); if (param_.solve_welleq_initially_ && initial_assembly) { // solve the well equations as a pre-processing step asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state); diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index eaaa5bf04..176e8a758 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -87,14 +87,6 @@ namespace Opm { ReservoirState& reservoir_state, WellState& well_state); - /// Assemble the residual and Jacobian of the nonlinear system. - /// \param[in] reservoir_state reservoir state variables - /// \param[in, out] well_state well state variables - /// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep - void assemble(const ReservoirState& reservoir_state, - WellState& well_state, - const bool initial_assembly); - protected: @@ -159,6 +151,7 @@ namespace Opm { using Base::updateWellControls; using Base::computeWellConnectionPressures; using Base::addWellControlEq; + using Base::computePropertiesForWellConnectionPressures; std::vector computeRelPerm(const SolutionState& state) const; @@ -212,8 +205,12 @@ namespace Opm { const SolutionState& state, WellState& xw); - void computeWellConnectionPressures(const SolutionState& state, - const WellState& xw); + void computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf); void updateEquationsScaling(); @@ -229,6 +226,11 @@ namespace Opm { const std::vector phaseCondition() const {return this->phaseCondition_;} + void extractWellPerfProperties(const SolutionState& state, + std::vector& mob_perfcells, + std::vector& b_perfcells); + + // compute effective viscosities (mu_eff_) and effective b factors (b_eff_) using the ToddLongstaff model void computeEffectiveProperties(const SolutionState& state); diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index 03dcf8746..65fb18ee6 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -135,7 +135,7 @@ namespace Opm { std::vector vars0 = Base::variableStateInitials(x, xw); assert(int(vars0.size()) == fluid_.numPhases() + 2); - // Initial polymer concentration. + // Initial solvent concentration. if (has_solvent_) { const auto& solvent_saturation = x.getCellData( BlackoilSolventState::SSOL ); const int nc = solvent_saturation.size(); @@ -257,6 +257,10 @@ namespace Opm { BlackoilSolventModel::computeAccum(const SolutionState& state, const int aix ) { + + if (is_miscible_) { + computeEffectiveProperties(state); + } Base::computeAccum(state, aix); // Compute accumulation of the solvent @@ -378,10 +382,13 @@ namespace Opm { } template - void BlackoilSolventModel::computeWellConnectionPressures(const SolutionState& state, - const WellState& xw) + void BlackoilSolventModel::computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) { - if( ! Base::localWellsActive() ) return ; using namespace Opm::AutoDiffGrid; // 1. Compute properties required by computeConnectionPressureDelta(). @@ -417,8 +424,7 @@ namespace Opm { const PhaseUsage& pu = fluid_.phaseUsage(); DataBlock b(nperf, pu.num_phases); - std::vector rsmax_perf(nperf, 0.0); - std::vector rvmax_perf(nperf, 0.0); + const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); if (pu.phase_used[BlackoilPhases::Aqua]) { b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; @@ -435,6 +441,8 @@ namespace Opm { b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; const V rssat = fluidRsSat(avg_press, perf_so, well_cells); rsmax_perf.assign(rssat.data(), rssat.data() + nperf); + } else { + rsmax_perf.assign(0.0, nperf); } V surf_dens_copy = superset(fluid_.surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); for (int phase = 1; phase < pu.num_phases; ++phase) { @@ -492,41 +500,18 @@ namespace Opm { const V rvsat = fluidRvSat(avg_press, perf_so, well_cells); rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); + } else { + rvmax_perf.assign(0.0, nperf); } // b and surf_dens_perf is row major, so can just copy data. - std::vector b_perf(b.data(), b.data() + nperf * pu.num_phases); - std::vector surf_dens_perf(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases); - - // Extract well connection depths. - const V depth = cellCentroidsZToEigen(grid_); - const V pdepth = subset(depth, well_cells); - std::vector perf_depth(pdepth.data(), pdepth.data() + nperf); - - // Gravity - double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); - - // 2. Compute densities - std::vector cd = - WellDensitySegmented::computeConnectionDensities( - wells(), xw, fluid_.phaseUsage(), - b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - // 3. Compute pressure deltas - std::vector cdp = - WellDensitySegmented::computeConnectionPressureDelta( - wells(), perf_depth, cd, grav); - - // 4. Store the results - Base::well_perforation_densities_ = Eigen::Map(cd.data(), nperf); - Base::well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); + b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); + surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases); } - - template void BlackoilSolventModel::updateState(const V& dx, ReservoirState& reservoir_state, @@ -771,6 +756,44 @@ namespace Opm { } + template + void + BlackoilSolventModel::extractWellPerfProperties(const SolutionState& state, + std::vector& mob_perfcells, + std::vector& b_perfcells) + { + Base::extractWellPerfProperties(state, mob_perfcells, b_perfcells); + if (has_solvent_) { + int gas_pos = fluid_.phaseUsage().phase_pos[Gas]; + const std::vector& well_cells = wops_.well_cells; + const int nperf = well_cells.size(); + // Gas and solvent is combinded and solved together + // The input in the well equation is then the + // total gas phase = hydro carbon gas + solvent gas + + // The total mobility is the sum of the solvent and gas mobiliy + mob_perfcells[gas_pos] += subset(rq_[solvent_pos_].mob, well_cells); + + // A weighted sum of the b-factors of gas and solvent are used. + const int nc = Opm::AutoDiffGrid::numCells(grid_); + + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + const ADB zero = ADB::constant(V::Zero(nc)); + const ADB& ss = state.solvent_saturation; + const ADB& sg = (active_[ Gas ] + ? state.saturation[ pu.phase_pos[ Gas ] ] + : zero); + + Selector zero_selector(ss.value() + sg.value(), Selector::Zero); + ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells); + V ones = V::Constant(nperf,1.0); + + b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos]; + b_perfcells[gas_pos] += (F_solvent * subset(rq_[solvent_pos_].b, well_cells)); + + } + } + template void BlackoilSolventModel::computeEffectiveProperties(const SolutionState& state) @@ -965,114 +988,20 @@ namespace Opm { const ADB& ss) const { std::vector pressures = Base::computePressures(po, sw, so, sg); - // The imiscible capillary pressure is evaluated using the total gas saturation (sg + ss) - std::vector pressures_imisc = Base::computePressures(po, sw, so, sg + ss); - - // Pressure effects on capillary pressure miscibility - const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_); - // Only the pcog is effected by the miscibility. Since pg = po + pcog, changing pg is eqvivalent - // to changing the gas pressure directly. - const int nc = cells_.size(); - const V ones = V::Constant(nc, 1.0); - pressures[Gas] = ( pmisc * pressures[Gas] + ((ones - pmisc) * pressures_imisc[Gas])); - return pressures; - } - - - template - void - BlackoilSolventModel::assemble(const ReservoirState& reservoir_state, - WellState& well_state, - const bool initial_assembly) - { - - using namespace Opm::AutoDiffGrid; - - // Possibly switch well controls and updating well state to - // get reasonable initial conditions for the wells - updateWellControls(well_state); - - // Create the primary variables. - SolutionState state = variableState(reservoir_state, well_state); - - if (initial_assembly) { - // Create the (constant, derivativeless) initial state. - SolutionState state0 = state; - makeConstantState(state0); - // Compute initial accumulation contributions - // and well connection pressures. - if (is_miscible_) { - computeEffectiveProperties(state0); - } - - computeAccum(state0, 0); - computeWellConnectionPressures(state0, well_state); - } - if (is_miscible_) { - computeEffectiveProperties(state); - } - - // -------- Mass balance equations -------- - assembleMassBalanceEq(state); - - // -------- Well equations ---------- - if ( ! wellsActive() ) { - return; - } - - V aliveWells; - - const int np = wells().number_of_phases; - std::vector cq_s(np, ADB::null()); - - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); - - std::vector mob_perfcells(np, ADB::null()); - std::vector b_perfcells(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - mob_perfcells[phase] = subset(rq_[phase].mob, well_cells); - b_perfcells[phase] = subset(rq_[phase].b, well_cells); - } if (has_solvent_) { - int gas_pos = fluid_.phaseUsage().phase_pos[Gas]; - // Gas and solvent is combinded and solved together - // The input in the well equation is then the - // total gas phase = hydro carbon gas + solvent gas - - // The total mobility is the sum of the solvent and gas mobiliy - mob_perfcells[gas_pos] += subset(rq_[solvent_pos_].mob, well_cells); - - // A weighted sum of the b-factors of gas and solvent are used. - const int nc = Opm::AutoDiffGrid::numCells(grid_); - - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - const ADB zero = ADB::constant(V::Zero(nc)); - const ADB& ss = state.solvent_saturation; - const ADB& sg = (active_[ Gas ] - ? state.saturation[ pu.phase_pos[ Gas ] ] - : zero); - - Selector zero_selector(ss.value() + sg.value(), Selector::Zero); - ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells); - V ones = V::Constant(nperf,1.0); - - b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos]; - b_perfcells[gas_pos] += (F_solvent * subset(rq_[solvent_pos_].b, well_cells)); + // The imiscible capillary pressure is evaluated using the total gas saturation (sg + ss) + std::vector pressures_imisc = Base::computePressures(po, sw, so, sg + ss); + // Pressure effects on capillary pressure miscibility + const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_); + // Only the pcog is effected by the miscibility. Since pg = po + pcog, changing pg is eqvivalent + // to changing the gas pressure directly. + const int nc = cells_.size(); + const V ones = V::Constant(nc, 1.0); + pressures[Gas] = ( pmisc * pressures[Gas] + ((ones - pmisc) * pressures_imisc[Gas])); } - if (param_.solve_welleq_initially_ && initial_assembly) { - // solve the well equations as a pre-processing step - Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); - } - Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); - Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state); - Base::addWellFluxEq(cq_s, state); - addWellContributionToMassBalanceEq(cq_s, state, well_state); - Base::addWellControlEq(state, well_state, aliveWells); - + return pressures; }