diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 255fae7f0..fe95bcbc1 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -131,7 +131,7 @@ list (APPEND PROGRAM_SOURCE_FILES # originally generated with the command: # find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort list (APPEND PUBLIC_HEADER_FILES -<<<<<<< HEAD + opm/autodiff/AdditionalObjectDeleter.hpp opm/autodiff/AutoDiffBlock.hpp opm/autodiff/AutoDiffHelpers.hpp @@ -197,7 +197,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/WellStateMultiSegment.hpp opm/autodiff/WellMultiSegment.hpp opm/autodiff/StandardWells.hpp - opm/autodiff/StandardWells_impl.hpp + opm/autodiff/StandardWellsSolvent.hpp opm/polymer/CompressibleTpfaPolymer.hpp opm/polymer/GravityColumnSolverPolymer.hpp opm/polymer/GravityColumnSolverPolymer_impl.hpp @@ -227,3 +227,4 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/SimulatorIncompTwophase.hpp ) + diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index b8ceae236..6e27f4c20 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -383,13 +383,6 @@ 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); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 363b5a3d7..323a346bd 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -758,81 +758,8 @@ namespace detail { } } - template - 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) - { - const int nperf = wells().well_connpos[wells().number_of_wells]; - const int nw = wells().number_of_wells; - - // Compute the average pressure in each well block - const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); - V avg_press = perf_press*0; - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; - const double p_avg = (perf_press[perf] + p_above)/2; - avg_press[perf] = p_avg; - } - } - - const std::vector& well_cells = stdWells().wellOps().well_cells; - - // Use cell values for the temperature as the wells don't knows its temperature yet. - const ADB perf_temp = subset(state.temperature, well_cells); - - // Compute b, rsmax, rvmax values for perforations. - // Evaluate the properties using average well block pressures - // and cell values for rs, rv, phase condition and temperature. - const ADB avg_press_ad = ADB::constant(avg_press); - std::vector perf_cond(nperf); - const std::vector& pc = phaseCondition(); - for (int perf = 0; perf < nperf; ++perf) { - perf_cond[perf] = pc[well_cells[perf]]; - } - const PhaseUsage& pu = fluid_.phaseUsage(); - DataBlock b(nperf, pu.num_phases); - 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; - } - assert(active_[Oil]); - const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); - if (pu.phase_used[BlackoilPhases::Liquid]) { - const ADB perf_rs = subset(state.rs, well_cells); - const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); - 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); - const V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); - 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. - b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); - // Surface density. - // The compute density segment wants the surface densities as - // an np * number of wells cells array - V rho = 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) { - rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); - } - surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); - } template diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index 505f7a3a6..02b7dd3ac 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -25,6 +25,7 @@ #include #include #include +#include namespace Opm { @@ -103,6 +104,9 @@ namespace Opm { const bool is_miscible_; std::vector mu_eff_; std::vector b_eff_; + StandardWellsSolvent std_wells_; + const StandardWellsSolvent& stdWells() const { return std_wells_; } + StandardWellsSolvent& stdWells() { return std_wells_; } // Need to declare Base members we want to use here. @@ -130,7 +134,7 @@ namespace Opm { // --------- Protected methods --------- // Need to declare Base members we want to use here. - using Base::stdWells; + // using Base::stdWells; using Base::wells; using Base::variableState; using Base::computeGasPressure; @@ -148,7 +152,7 @@ namespace Opm { using Base::updateWellControls; using Base::computeWellConnectionPressures; using Base::addWellControlEq; - using Base::computePropertiesForWellConnectionPressures; + // using Base::computePropertiesForWellConnectionPressures; std::vector computeRelPerm(const SolutionState& state) const; @@ -202,13 +206,6 @@ namespace Opm { const SolutionState& state, 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(); void diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index c292c0626..6145267df 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -89,7 +89,8 @@ namespace Opm { has_solvent_(has_solvent), solvent_pos_(detail::solventPos(fluid.phaseUsage())), solvent_props_(solvent_props), - is_miscible_(is_miscible) + is_miscible_(is_miscible), + std_wells_(wells_arg, solvent_props) { if (has_solvent_) { @@ -381,132 +382,6 @@ namespace Opm { } } - template - 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) - { - 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; - const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); - - // Compute the average pressure in each well block - const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); - V avg_press = perf_press*0; - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; - const double p_avg = (perf_press[perf] + p_above)/2; - avg_press[perf] = p_avg; - } - } - - // Use cell values for the temperature as the wells don't knows its temperature yet. - const ADB perf_temp = subset(state.temperature, well_cells); - - // Compute b, rsmax, rvmax values for perforations. - // Evaluate the properties using average well block pressures - // and cell values for rs, rv, phase condition and temperature. - const ADB avg_press_ad = ADB::constant(avg_press); - std::vector perf_cond(nperf); - const std::vector& pc = phaseCondition(); - for (int perf = 0; perf < nperf; ++perf) { - perf_cond[perf] = pc[well_cells[perf]]; - } - - const PhaseUsage& pu = fluid_.phaseUsage(); - DataBlock b(nperf, pu.num_phases); - - 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; - } - - assert(active_[Oil]); - assert(active_[Gas]); - const ADB perf_rv = subset(state.rv, well_cells); - const ADB perf_rs = subset(state.rs, well_cells); - const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); - if (pu.phase_used[BlackoilPhases::Liquid]) { - const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); - //const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value(); - 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) { - if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) { - continue; // the gas surface density is added after the solvent is accounted for. - } - surf_dens_copy += superset(fluid_.surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); - } - - if (pu.phase_used[BlackoilPhases::Vapour]) { - // Unclear wether the effective or the pure values should be used for the wells - // the current usage of unmodified properties values gives best match. - //V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value(); - V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); - V rhog = fluid_.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells); - if (has_solvent_) { - - const V bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value(); - //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value(); - - // A weighted sum of the b-factors of gas and solvent are used. - const int nc = Opm::AutoDiffGrid::numCells(grid_); - - 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); - V F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value(); - - V injectedSolventFraction = Eigen::Map(&xw.solventFraction()[0], nperf); - - V isProducer = V::Zero(nperf); - V ones = V::Constant(nperf,1.0); - for (int w = 0; w < nw; ++w) { - if(wells().type[w] == PRODUCER) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - isProducer[perf] = 1; - } - } - } - - F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction; - - bg = bg * (ones - F_solvent); - bg = bg + F_solvent * bs; - - const V& rhos = solvent_props_.solventSurfaceDensity(well_cells); - rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos); - } - b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; - surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases); - - 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. - 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); - } diff --git a/opm/autodiff/StandardWellsSolvent.hpp b/opm/autodiff/StandardWellsSolvent.hpp new file mode 100644 index 000000000..9906eae3b --- /dev/null +++ b/opm/autodiff/StandardWellsSolvent.hpp @@ -0,0 +1,61 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + 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_STANDARDWELLSSOLVENT_HEADER_INCLUDED +#define OPM_STANDARDWELLSSOLVENT_HEADER_INCLUDED + +#include +#include + +namespace Opm { + + + /// Class for handling the standard well model for solvent model + class StandardWellsSolvent : public StandardWells + { + public: + + using Base = StandardWells; + + // --------- Public methods --------- + explicit StandardWellsSolvent(const Wells* wells, const SolventPropsAdFromDeck& solvent_props); + + template + void computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& pc, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf); + protected: + const SolventPropsAdFromDeck& solvent_props_; + + }; + + +} // namespace Opm + +#include "StandardWellsSolvent_impl.hpp" + +#endif diff --git a/opm/autodiff/StandardWellsSolvent_impl.hpp b/opm/autodiff/StandardWellsSolvent_impl.hpp new file mode 100644 index 000000000..38b538675 --- /dev/null +++ b/opm/autodiff/StandardWellsSolvent_impl.hpp @@ -0,0 +1,181 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + 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 +{ + + + + + + + StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg, + const SolventPropsAdFromDeck& solvent_props) + : Base(wells_arg) + , solvent_props_(solvent_props) + { + } + + + + + + template + void + StandardWellsSolvent:: + computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& pc, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) + { + // 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; + + // Compute the average pressure in each well block + const Vector perf_press = Eigen::Map(xw.perfPress().data(), nperf); + Vector avg_press = perf_press*0; + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + avg_press[perf] = p_avg; + } + } + + const std::vector& well_cells = wellOps().well_cells; + + // Use cell values for the temperature as the wells don't knows its temperature yet. + const ADB perf_temp = subset(state.temperature, well_cells); + + // Compute b, rsmax, rvmax values for perforations. + // Evaluate the properties using average well block pressures + // and cell values for rs, rv, phase condition and temperature. + const ADB avg_press_ad = ADB::constant(avg_press); + std::vector perf_cond(nperf); + for (int perf = 0; perf < nperf; ++perf) { + perf_cond[perf] = pc[well_cells[perf]]; + } + + const PhaseUsage& pu = fluid.phaseUsage(); + DataBlock b(nperf, pu.num_phases); + + const Vector 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; + } + + assert(active[Oil]); + assert(active[Gas]); + const ADB perf_rv = subset(state.rv, well_cells); + const ADB perf_rs = subset(state.rs, well_cells); + const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); + if (pu.phase_used[BlackoilPhases::Liquid]) { + const Vector bo = fluid.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); + //const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; + // const Vector rssat = fluidRsSat(avg_press, perf_so, well_cells); + const Vector rssat = fluid.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + 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) { + if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) { + continue; // the gas surface density is added after the solvent is accounted for. + } + surf_dens_copy += superset(fluid.surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + // Unclear wether the effective or the pure values should be used for the wells + // the current usage of unmodified properties values gives best match. + //V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value(); + Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); + Vector rhog = fluid.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells); + // to handle solvent related + { + + const Vector bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value(); + //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value(); + + // number of cells + const int nc = state.pressure.size(); + + const ADB zero = ADB::constant(Vector::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); + Vector F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value(); + + Vector injectedSolventFraction = Eigen::Map(&xw.solventFraction()[0], nperf); + + Vector isProducer = Vector::Zero(nperf); + Vector ones = Vector::Constant(nperf,1.0); + for (int w = 0; w < nw; ++w) { + if(wells().type[w] == PRODUCER) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + isProducer[perf] = 1; + } + } + } + + F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction; + + bg = bg * (ones - F_solvent); + bg = bg + F_solvent * bs; + + const Vector& rhos = solvent_props_.solventSurfaceDensity(well_cells); + rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos); + } + b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; + surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases); + + // const Vector rvsat = fluidRvSat(avg_press, perf_so, well_cells); + const Vector rvsat = fluid.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + 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. + 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); + } + + +}