diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 920aa6ea0..865a36641 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -108,6 +108,7 @@ namespace Opm { static const int solventSaturationIdx = Indices::solventSaturationIdx; static constexpr bool has_solvent_ = getPropValue(); static constexpr bool has_polymer_ = getPropValue(); + static constexpr bool has_energy_ = getPropValue(); // TODO: where we should put these types, WellInterface or Well Model? // or there is some other strategy, like TypeTag @@ -509,6 +510,8 @@ namespace Opm { void assignGroupGuideRates(const Group& group, const std::unordered_map& groupGuideRates, data::GroupData& gdata) const; + + void computeWellTemperature(); }; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 7d718daa3..1c6e2d124 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -540,6 +540,9 @@ namespace Opm { if (terminal_output_) { global_deferredLogger.logMessages(); } + + //reporting output temperatures + this->computeWellTemperature(); } @@ -703,12 +706,6 @@ namespace Opm { const SummaryState& summaryState) { std::vector cellPressures(this->local_num_cells_, 0.0); - std::vector cellTemperatures(this->local_num_cells_, 0.0); - const int np = numPhases(); - std::vector> cellInternalEnergy(this->local_num_cells_, std::vector(np)); - std::vector> cellBinv(this->local_num_cells_, std::vector(np)); - std::vector> cellDensity(this->local_num_cells_, std::vector(np)); - ElementContext elemCtx(ebosSimulator_); const auto& gridView = ebosSimulator_.vanguard().gridView(); @@ -725,16 +722,6 @@ namespace Opm { elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); const auto& fs = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0).fluidState(); - - cellTemperatures[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)] = fs.temperature(/*phaseIdx*/0).value(); - - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) - { - cellInternalEnergy[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)][phaseIdx] = fs.enthalpy(phaseIdx).value() - fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value(); - cellBinv[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)][phaseIdx] = fs.invB(phaseIdx).value(); - cellDensity[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)][phaseIdx] = fs.density(phaseIdx).value(); - } - // copy of get perfpressure in Standard well except for value double& perf_pressure = cellPressures[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)]; if (Indices::oilEnabled) { @@ -746,7 +733,7 @@ namespace Opm { } } - well_state_.init(cellPressures, cellTemperatures, cellInternalEnergy, cellBinv, cellDensity, schedule(), wells_ecl_, local_parallel_well_info_, timeStepIdx, + well_state_.init(cellPressures, schedule(), wells_ecl_, local_parallel_well_info_, timeStepIdx, &previous_well_state_, phase_usage_, well_perf_data_, summaryState, globalNumWells); } @@ -3081,4 +3068,57 @@ namespace Opm { prod = xgr.production; inj = xgr.injection; } + + template + void + BlackoilWellModel:: + computeWellTemperature() + { + if (!has_energy_) + return; + + int np = numPhases(); + double cellInternalEnergy; + double cellBinv; + double cellDensity; + double perfPhaseRate; + const int nw = numLocalWells(); + for (auto wellID = 0*nw; wellID < nw; ++wellID) { + const Well& well = wells_ecl_[wellID]; + if (well.isInjector()) + continue; + + int connpos = 0; + for (int i = 0; i < wellID; ++i) { + connpos += well_perf_data_[i].size(); + } + connpos *= np; + double weighted_temperature = 0.0; + double total_weight = 0.0; + + auto& well_info = *local_parallel_well_info_[wellID]; + const int num_perf_this_well = well_info.communication().sum(well_perf_data_[wellID].size()); + + for (int perf = 0; perf < num_perf_this_well; ++perf) { + const int cell_idx = well_perf_data_[wellID][perf].cell_index; + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + double cellTemperatures = fs.temperature(/*phaseIdx*/0).value(); + double weight_factor = 0.0; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + cellInternalEnergy = fs.enthalpy(phaseIdx).value() - fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value(); + cellBinv = fs.invB(phaseIdx).value(); + cellDensity = fs.density(phaseIdx).value(); + perfPhaseRate = well_state_.perfPhaseRates()[connpos + perf*np + phaseIdx ]; + weight_factor += cellDensity * perfPhaseRate/cellBinv * cellInternalEnergy/cellTemperatures; + } + total_weight += weight_factor; + weighted_temperature += weight_factor * cellTemperatures; + } + weighted_temperature = well_info.communication().sum(weighted_temperature); + total_weight = well_info.communication().sum(total_weight); + well_state_.temperature()[wellID] = weighted_temperature/total_weight; + } + } } // namespace Opm diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index e2230f526..0e354ea63 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -39,7 +39,7 @@ namespace Opm { /// The state of a set of wells. - class WellState + class WellState { public: typedef std::array< int, 3 > mapentry_t; @@ -51,11 +51,6 @@ namespace Opm /// perfRates() field is filled with zero, and perfPress() /// with -1e100. void init(const std::vector& cellPressures, - const std::vector& cellTemperatures, - const std::vector>& cellInternalEnergy, - const std::vector>& cellBinv, - const std::vector>& cellDensity, - const std::vector& perforationRates, const std::vector& wells_ecl, const std::vector& parallel_well_info, const PhaseUsage& pu, @@ -78,13 +73,13 @@ namespace Opm bhp_.resize(nw, 0.0); thp_.resize(nw, 0.0); temperature_.resize(nw, 273.15 + 15.56); // standard condition temperature - wellrates_.resize(nw * np, 0.0); + wellrates_.resize(nw * np, 0.0); int connpos = 0; for (int w = 0; w < nw; ++w) { const Well& well = wells_ecl[w]; // Initialize bhp(), thp(), wellRates(), temperature(). - initSingleWell(cellPressures, cellTemperatures, cellInternalEnergy, cellBinv, cellDensity, perforationRates, w, well, *parallel_well_info[w], pu, summary_state); + initSingleWell(cellPressures, w, well, *parallel_well_info[w], pu, summary_state); // Setup wellname -> well index mapping. const int num_perf_this_well = well_perf_data[w].size(); @@ -358,11 +353,6 @@ namespace Opm sizes.data(), displ.data(), 0); } void initSingleWell(const std::vector& cellPressures, - const std::vector& cellTemperatures, - const std::vector>& cellInternalEnergy, - const std::vector>& cellBinv, - const std::vector>& cellDensity, - const std::vector& perforationRates, const int w, const Well& well, const ParallelWellInfo& well_info, @@ -383,32 +373,6 @@ namespace Opm } const int num_perf_this_well = well_info.communication().sum(well_perf_data_[w].size()); - if ( well.isProducer() ) { - int sz = perforationRates.size(); - int connpos = 0; - for (int i = 0; i < w; ++i){ - connpos += well_perf_data_[i].size(); - } - connpos *= np; - double weighted_temperature = 0.0; - double total_weight = 0.0; - if (sz > 0){ - for (int i = 0; i < num_perf_this_well; ++i){ - int perf_cell_index = well_perf_data_[w][i].cell_index; - double weight_factor = 0.0; - for (int p = 0; p < np; ++p) { - weight_factor += cellDensity[perf_cell_index][p] * perforationRates[connpos + i*np + p]/cellBinv[perf_cell_index][p] * cellInternalEnergy[perf_cell_index][p]/cellTemperatures[perf_cell_index]; - } - total_weight += weight_factor; - weighted_temperature += weight_factor * cellTemperatures[perf_cell_index]; - } - temperature_[w] = weighted_temperature/total_weight; - } - else { - temperature_[w]=0; - } - } - if ( num_perf_this_well == 0 ) { // No perforations of the well. Initialize to zero. bhp_[w] = 0.; @@ -525,7 +489,7 @@ namespace Opm protected: std::vector> well_perf_data_; - std::vector parallel_well_info_; + std::vector parallel_well_info_; }; } // namespace Opm diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp index 5b74f6e77..cb5590cae 100644 --- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp +++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp @@ -70,10 +70,6 @@ namespace Opm /// to give useful initial values to the bhp(), wellRates() /// and perfPhaseRates() fields, depending on controls void init(const std::vector& cellPressures, - const std::vector& cellTemperatures, - const std::vector>& cellInternalEnergy, - const std::vector>& cellBinv, - const std::vector>& cellDensity, const Schedule& schedule, const std::vector& wells_ecl, const std::vector& parallel_well_info, @@ -84,10 +80,8 @@ namespace Opm const SummaryState& summary_state, const int globalNumberOfWells) { - std::vector perforationRates; - perforationRates = prevState->perfphaserates_; // call init on base class - BaseType :: init(cellPressures, cellTemperatures, cellInternalEnergy, cellBinv, cellDensity, perforationRates, wells_ecl, parallel_well_info, pu, well_perf_data, summary_state); + BaseType :: init(cellPressures, wells_ecl, parallel_well_info, pu, well_perf_data, summary_state); for (const auto& winfo: parallel_well_info) { @@ -185,7 +179,6 @@ namespace Opm perfRateBrine_.clear(); perfRateBrine_.resize(nperf, 0.0); - // intialize wells that have been there before // order may change so the mapping is based on the well name if (prevState && !prevState->wellMap().empty()) { @@ -355,8 +348,7 @@ namespace Opm const int globalNumWells) { const std::vector tmp(numCells, 0.0); // <- UGLY HACK to pass the size - std::vector> tmp1(numCells, std::vector( numPhases()));// <- UGLY HACK to pass the size - init(tmp, tmp, tmp1, tmp1, tmp1, schedule, wells_ecl, parallel_well_info, 0, nullptr, pu, well_perf_data, summary_state, globalNumWells); + init(tmp, schedule, wells_ecl, parallel_well_info, 0, nullptr, pu, well_perf_data, summary_state, globalNumWells); if (handle_ms_well) { initWellStateMSWell(wells_ecl, pu, nullptr);