From 7b81facfb053217e4848b025817ca6adab142417 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 16 Mar 2016 15:35:26 +0100 Subject: [PATCH] Make use of extractWellPerfProperties to avoid code duplication The following changes are done in order to remove the duplicated code in assemble(). - extractWellPerfProperties takes SolutionState as input (only used in the solvent model) - the computation of effective parameters is moved to computeAccum() With these changes the solvent model can use assemble() from the base model. --- opm/autodiff/BlackoilModelBase.hpp | 3 +- opm/autodiff/BlackoilModelBase_impl.hpp | 5 +- opm/autodiff/BlackoilSolventModel.hpp | 13 +- opm/autodiff/BlackoilSolventModel_impl.hpp | 141 +++++++-------------- 4 files changed, 53 insertions(+), 109 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 2551bd7fb..6a6a5a84f 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -392,7 +392,8 @@ namespace Opm { void extractWellPerfProperties(std::vector& mob_perfcells, - std::vector& b_perfcells) const; + std::vector& b_perfcells, + const SolutionState&) const; bool solveWellEq(const std::vector& mob_perfcells, diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index ed89c4481..fa6d2f38d 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -954,7 +954,7 @@ namespace detail { std::vector mob_perfcells; std::vector b_perfcells; - asImpl().extractWellPerfProperties(mob_perfcells, b_perfcells); + asImpl().extractWellPerfProperties(mob_perfcells, b_perfcells, state); 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); @@ -1102,7 +1102,8 @@ namespace detail { template void BlackoilModelBase::extractWellPerfProperties(std::vector& mob_perfcells, - std::vector& b_perfcells) const + std::vector& b_perfcells, + const SolutionState&) const { // If we have wells, extract the mobilities and b-factors for // the well-perforated cells. diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index eaaa5bf04..987a93db1 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: @@ -229,6 +221,11 @@ namespace Opm { const std::vector phaseCondition() const {return this->phaseCondition_;} + void extractWellPerfProperties(std::vector& mob_perfcells, + std::vector& b_perfcells, + const SolutionState& state); + + // 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..34e9f1b54 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 @@ -771,6 +775,44 @@ namespace Opm { } + template + void + BlackoilSolventModel::extractWellPerfProperties(std::vector& mob_perfcells, + std::vector& b_perfcells, + const SolutionState& state) + { + Base::extractWellPerfProperties(mob_perfcells, b_perfcells, state); + 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) @@ -979,103 +1021,6 @@ namespace Opm { } - 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)); - - } - 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); - - } - - }