diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index 3ce6dfaa9..111375766 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -625,6 +625,21 @@ namespace Opm { wellModel().updateWellState(dwells, dpMaxRel(), well_state); + for( auto w = 0; w < wells().number_of_wells; ++w ) { + if (wells().type[w] == INJECTOR) { + continue; + } + + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf ) { + int wc = wells().well_cells[perf]; + if ( (ss[wc] + sg[wc]) > 0) { + well_state.solventFraction()[perf] = ss[wc] / (ss[wc] + sg[wc]); + } + } + } + + + // Update phase conditions used for property calculations. updatePhaseCondFromPrimalVariable(reservoir_state); } diff --git a/opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp b/opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp index eff83810b..c2061e910 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp @@ -33,6 +33,37 @@ namespace Opm /// One solvent fraction per well connection std::vector& solventFraction() { return solvent_fraction_; } const std::vector& solventFraction() const { return solvent_fraction_; } + + data::Wells report(const PhaseUsage &pu) const override { + data::Wells res = WellStateFullyImplicitBlackoil::report(pu); + + const int nw = WellState::numWells(); + + // If there are now wells numPhases throws a floating point + // exception. + if (nw == 0) { + return res; + } + + const int np = BaseType::numPhases(); + + assert( np == 3 ); // the solvent model assumes 3 phases in the base model + + // completions aren't supported yet + for( auto w = 0; w < nw; ++w ) { + using rt = data::Rates::opt; + double solvent_well_rate = 0.0; + for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) { + auto solvent_rate_this = BaseType::perfPhaseRates()[np*perf + pu.phase_pos[BlackoilPhases::Vapour]] * solventFraction()[perf]; + solvent_well_rate += solvent_rate_this; + } + + res.at( wells_->name[ w ]).rates.set( rt::solvent, solvent_well_rate ); + + } + + return res; + } private: std::vector solvent_fraction_; };