diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 60dcb2627..7a655fdc3 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -117,8 +117,6 @@ namespace Opm { typedef ISTLSolver< MatrixBlockType, VectorBlockType, Indices::pressureSwitchIdx > ISTLSolverType; //typedef typename SolutionVector :: value_type PrimaryVariables ; - typedef Opm::FIPData FIPDataType; - // --------- Public methods --------- /// Construct the model. It will retain references to the @@ -999,144 +997,16 @@ namespace Opm { return computeFluidInPlace(fipnum); } + /// Should not be called std::vector > - computeFluidInPlace(const std::vector& fipnum) const - { - const auto& comm = grid_.comm(); - const auto& gridView = ebosSimulator().gridView(); - const int nc = gridView.size(/*codim=*/0); - int ntFip = *std::max_element(fipnum.begin(), fipnum.end()); - ntFip = comm.max(ntFip); - - std::vector tpv(ntFip, 0.0); - std::vector hcpv(ntFip, 0.0); - - std::vector > regionValues(ntFip, std::vector(FIPDataType::fipValues,0.0)); - - for (int i = 0; i(); - for (auto elemIt = gridView.template begin(); - elemIt != elemEndIt; - ++elemIt) - { - elemCtx.updatePrimaryStencil(*elemIt); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - - const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& fs = intQuants.fluidState(); - - const int regionIdx = fipnum[cellIdx] - 1; - if (regionIdx < 0) { - // the given cell is not attributed to any region - continue; - } - - // calculate the pore volume of the current cell. Note that the porosity - // returned by the intensive quantities is defined as the ratio of pore - // space to total cell volume and includes all pressure dependent (-> - // rock compressibility) and static modifiers (MULTPV, MULTREGP, NTG, - // PORV, MINPV and friends). Also note that because of this, the porosity - // returned by the intensive quantities can be outside of the physical - // range [0, 1] in pathetic cases. - const double pv = - ebosSimulator_.model().dofTotalVolume(cellIdx) - * intQuants.porosity().value(); - - for (unsigned phase = 0; phase < FluidSystem::numPhases; ++phase) { - if (!FluidSystem::phaseIsActive(phase)) { - continue; - } - - const double b = fs.invB(phase).value(); - const double s = fs.saturation(phase).value(); - const unsigned flowCanonicalPhaseIdx = ebosPhaseToFlowCanonicalPhaseIdx(phase); - - fip_.fip[flowCanonicalPhaseIdx][cellIdx] = b * s * pv; - regionValues[regionIdx][flowCanonicalPhaseIdx] += fip_.fip[flowCanonicalPhaseIdx][cellIdx]; - } - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - // Account for gas dissolved in oil and vaporized oil - fip_.fip[FIPDataType::FIP_DISSOLVED_GAS][cellIdx] = fs.Rs().value() * fip_.fip[FIPDataType::FIP_LIQUID][cellIdx]; - fip_.fip[FIPDataType::FIP_VAPORIZED_OIL][cellIdx] = fs.Rv().value() * fip_.fip[FIPDataType::FIP_VAPOUR][cellIdx]; - - regionValues[regionIdx][FIPData::FIP_DISSOLVED_GAS] += fip_.fip[FIPData::FIP_DISSOLVED_GAS][cellIdx]; - regionValues[regionIdx][FIPData::FIP_VAPORIZED_OIL] += fip_.fip[FIPData::FIP_VAPORIZED_OIL][cellIdx]; - } - - const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value(); - - tpv[regionIdx] += pv; - hcpv[regionIdx] += pv * hydrocarbon; - } - - // sum tpv (-> total pore volume of the regions) and hcpv (-> pore volume of the - // the regions that is occupied by hydrocarbons) - comm.sum(tpv.data(), tpv.size()); - comm.sum(hcpv.data(), hcpv.size()); - - for (auto elemIt = gridView.template begin(); - elemIt != elemEndIt; - ++elemIt) - { - const auto& elem = *elemIt; - - elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - - unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - const int regionIdx = fipnum[cellIdx] - 1; - if (regionIdx < 0) { - // the cell is not attributed to any region. ignore it! - continue; - } - - const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& fs = intQuants.fluidState(); - - // calculate the pore volume of the current cell. Note that the - // porosity returned by the intensive quantities is defined as the - // ratio of pore space to total cell volume and includes all pressure - // dependent (-> rock compressibility) and static modifiers (MULTPV, - // MULTREGP, NTG, PORV, MINPV and friends). Also note that because of - // this, the porosity returned by the intensive quantities can be - // outside of the physical range [0, 1] in pathetic cases. - const double pv = - ebosSimulator_.model().dofTotalVolume(cellIdx) - * intQuants.porosity().value(); - - fip_.fip[FIPDataType::FIP_PV][cellIdx] = pv; - const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value(); - - //Compute hydrocarbon pore volume weighted average pressure. - //If we have no hydrocarbon in region, use pore volume weighted average pressure instead - if (hcpv[regionIdx] > 1e-10) { - fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[regionIdx]; - } else { - fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() / tpv[regionIdx]; - } - - regionValues[regionIdx][FIPDataType::FIP_PV] += fip_.fip[FIPDataType::FIP_PV][cellIdx]; - regionValues[regionIdx][FIPDataType::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx]; - } - - // sum the results over all processes - for(int regionIdx=0; regionIdx < ntFip; ++regionIdx) { - comm.sum(regionValues[regionIdx].data(), regionValues[regionIdx].size()); - } - + computeFluidInPlace(const std::vector& /*fipnum*/) const + { + //assert(true) + //return an empty vector + std::vector > regionValues(0, std::vector(0,0.0)); return regionValues; } - const FIPDataType& getFIPData() const { - return fip_; - } - const Simulator& ebosSimulator() const { return ebosSimulator_; } @@ -1176,7 +1046,6 @@ namespace Opm { std::vector> residual_norms_history_; double current_relaxation_; BVector dx_old_; - mutable FIPDataType fip_; public: /// return the StandardWells object @@ -1186,20 +1055,6 @@ namespace Opm { const BlackoilWellModel& wellModel() const { return well_model_; } - int ebosPhaseToFlowCanonicalPhaseIdx( const int phaseIdx ) const - { - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::waterPhaseIdx == phaseIdx) - return Water; - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::oilPhaseIdx == phaseIdx) - return Oil; - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::gasPhaseIdx == phaseIdx) - return Gas; - - assert(phaseIdx < 3); - // for other phases return the index - return phaseIdx; - } - void beginReportStep() { ebosSimulator_.problem().beginEpisode(); diff --git a/opm/autodiff/BlackoilOutputEbos.hpp b/opm/autodiff/BlackoilOutputEbos.hpp index fccf45500..7250bab59 100644 --- a/opm/autodiff/BlackoilOutputEbos.hpp +++ b/opm/autodiff/BlackoilOutputEbos.hpp @@ -136,13 +136,8 @@ namespace Opm const double nextstep = -1.0, const SimulatorReport& simulatorReport = SimulatorReport()) { - data::Solution fip{}; - if( output_ ) { - // Get FIP dat - getSummaryData( fip, phaseUsage_, physicalModel, summaryConfig() ); - // Add TCPU if simulatorReport is not defaulted. const double totalSolverTime = simulatorReport.solver_time; @@ -169,7 +164,7 @@ namespace Opm const WellStateFullyImplicitBlackoil& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState; // The writeOutput expects a local data::solution vector and a global data::well vector. - ebosSimulator_.problem().writeOutput( wellState.report(phaseUsage_), timer.simulationTimeElapsed(), substep, totalSolverTime, nextstep, fip); + ebosSimulator_.problem().writeOutput( wellState.report(phaseUsage_), timer.simulationTimeElapsed(), substep, totalSolverTime, nextstep); } } @@ -222,157 +217,19 @@ namespace Opm } } - bool requireFIPNUM() const - { return summaryConfig().requireFIPNUM(); } - const Grid& grid() { return ebosSimulator_.gridManager().grid(); } const Schedule& schedule() const { return ebosSimulator_.gridManager().schedule(); } - const SummaryConfig& summaryConfig() const - { return ebosSimulator_.gridManager().summaryConfig(); } - const EclipseState& eclState() const { return ebosSimulator_.gridManager().eclState(); } bool isRestart() const { const auto& initconfig = eclState().getInitConfig(); return initconfig.restartRequested(); - } - - private: - /** - * Checks if the summaryConfig has a keyword with the standardized field, region, or block prefixes. - */ - inline bool hasFRBKeyword(const SummaryConfig& summaryConfig, const std::string keyword) { - std::string field_kw = "F" + keyword; - std::string region_kw = "R" + keyword; - std::string block_kw = "B" + keyword; - return summaryConfig.hasKeyword(field_kw) - || summaryConfig.hasKeyword(region_kw) - || summaryConfig.hasKeyword(block_kw); - } - - - /** - * Returns the data as asked for in the summaryConfig - */ - template - void getSummaryData(data::Solution& output, - const Opm::PhaseUsage& phaseUsage, - const Model& physicalModel, - const SummaryConfig& summaryConfig) { - - typedef typename Model::FIPDataType FIPDataType; - typedef typename FIPDataType::VectorType VectorType; - - FIPDataType fd = physicalModel.getFIPData(); - - //Get shorthands for water, oil, gas - const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua]; - const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid]; - const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour]; - - /** - * Now process all of the summary config files - */ - // Water in place - if (aqua_active && hasFRBKeyword(summaryConfig, "WIP")) { - output.insert("WIP", - Opm::UnitSystem::measure::volume, - std::move( fd.fip[ FIPDataType::FIP_AQUA ] ), - data::TargetType::SUMMARY ); - } - if (liquid_active) { - const VectorType& oipl = fd.fip[FIPDataType::FIP_LIQUID]; - VectorType oip ( oipl ); - const size_t size = oip.size(); - - const VectorType& oipg = vapour_active ? fd.fip[FIPDataType::FIP_VAPORIZED_OIL] : VectorType(size, 0.0); - if( vapour_active ) - { - // oip = oipl + oipg - for( size_t i=0; i #include #include -#include #include #include @@ -119,7 +118,6 @@ public: is_parallel_run_ = ( info.communicator().size() > 1 ); } #endif - createLocalFipnum(); } /// Run the simulation. @@ -198,20 +196,11 @@ public: auto solver = createSolver(well_model); - // Compute orignal fluid in place if this has not been done yet - if (originalFluidInPlace_.data.empty()) { - originalFluidInPlace_ = computeFluidInPlace(*solver); - } - // write the inital state at the report stage if (timer.initialStep()) { Dune::Timer perfTimer; perfTimer.start(); - if (terminal_output_) { - //outputFluidInPlace(timer, originalFluidInPlace_); - } - // No per cell data is written for initial step, but will be // for subsequent steps, when we have started simulating output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model() ); @@ -250,8 +239,7 @@ public: events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) || events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum()); - stepReport = adaptiveTimeStepping->step( timer, *solver, dummy_state, wellStateDummy, event, output_writer_, - output_writer_.requireFIPNUM() ? &fipnum_ : nullptr ); + stepReport = adaptiveTimeStepping->step( timer, *solver, dummy_state, wellStateDummy, event, output_writer_, nullptr ); report += stepReport; failureReport_ += adaptiveTimeStepping->failureReport(); } @@ -288,8 +276,6 @@ public: // Increment timer, remember well state. ++timer; - // Compute current fluid in place. - const auto currentFluidInPlace = computeFluidInPlace(*solver); if (terminal_output_ ) { @@ -343,150 +329,6 @@ protected: return std::unique_ptr(new Solver(solver_param_, std::move(model))); } - - void createLocalFipnum() - { - const std::vector& fipnum_global = eclState().get3DProperties().getIntGridProperty("FIPNUM").getData(); - // Get compressed cell fipnum. - fipnum_.resize(Opm::UgGridHelpers::numCells(grid())); - if (fipnum_global.empty()) { - std::fill(fipnum_.begin(), fipnum_.end(), 0); - } else { - for (size_t c = 0; c < fipnum_.size(); ++c) { - fipnum_[c] = fipnum_global[Opm::UgGridHelpers::globalCell(grid())[c]]; - } - } - } - - - void FIPUnitConvert(const UnitSystem& units, - std::vector>& fip) - { - for (size_t i = 0; i < fip.size(); ++i) { - FIPUnitConvert(units, fip[i]); - } - } - - - void FIPUnitConvert(const UnitSystem& units, - std::vector& fip) - { - if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) { - fip[0] = unit::convert::to(fip[0], unit::stb); - fip[1] = unit::convert::to(fip[1], unit::stb); - fip[2] = unit::convert::to(fip[2], 1000*unit::cubic(unit::feet)); - fip[3] = unit::convert::to(fip[3], 1000*unit::cubic(unit::feet)); - fip[4] = unit::convert::to(fip[4], unit::stb); - fip[5] = unit::convert::to(fip[5], unit::stb); - fip[6] = unit::convert::to(fip[6], unit::psia); - } - else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) { - fip[6] = unit::convert::to(fip[6], unit::barsa); - } - else { - OPM_THROW(std::runtime_error, "Unsupported unit type for fluid in place output."); - } - } - - - std::vector FIPTotals(const std::vector>& fip) - { - std::vector totals(7,0.0); - for (int i = 0; i < 5; ++i) { - for (size_t reg = 0; reg < fip.size(); ++reg) { - totals[i] += fip[reg][i]; - } - } - - const auto& gridView = ebosSimulator_.gridManager().gridView(); - const auto& comm = gridView.comm(); - double pv_hydrocarbon_sum = 0.0; - double p_pv_hydrocarbon_sum = 0.0; - - ElementContext elemCtx(ebosSimulator_); - const auto& elemEndIt = gridView.template end(); - for (auto elemIt = gridView.template begin(); - elemIt != elemEndIt; - ++elemIt) - { - const auto& elem = *elemIt; - if (elem.partitionType() != Dune::InteriorEntity) { - continue; - } - - elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - - const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& fs = intQuants.fluidState(); - - const double p = fs.pressure(FluidSystem::oilPhaseIdx).value(); - const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value(); - - // calculate the pore volume of the current cell. Note that the - // porosity returned by the intensive quantities is defined as the - // ratio of pore space to total cell volume and includes all pressure - // dependent (-> rock compressibility) and static modifiers (MULTPV, - // MULTREGP, NTG, PORV, MINPV and friends). Also note that because of - // this, the porosity returned by the intensive quantities can be - // outside of the physical range [0, 1] in pathetic cases. - const double pv = - ebosSimulator_.model().dofTotalVolume(cellIdx) - * intQuants.porosity().value(); - - totals[5] += pv; - pv_hydrocarbon_sum += pv*hydrocarbon; - p_pv_hydrocarbon_sum += p*pv*hydrocarbon; - } - - pv_hydrocarbon_sum = comm.sum(pv_hydrocarbon_sum); - p_pv_hydrocarbon_sum = comm.sum(p_pv_hydrocarbon_sum); - totals[5] = comm.sum(totals[5]); - totals[6] = (p_pv_hydrocarbon_sum / pv_hydrocarbon_sum); - - return totals; - } - - - struct FluidInPlace - { - std::vector> data; - std::vector totals; - }; - - - FluidInPlace computeFluidInPlace(const Solver& solver) - { - FluidInPlace fip; - fip.data = solver.computeFluidInPlace(fipnum_); - fip.totals = FIPTotals(fip.data); - FIPUnitConvert(eclState().getUnits(), fip.data); - FIPUnitConvert(eclState().getUnits(), fip.totals); - return fip; - } - - - void outputFluidInPlace(const SimulatorTimer& timer, - const FluidInPlace& currentFluidInPlace) - { - if (!timer.initialStep()) { - const std::string version = moduleVersionName(); - outputTimestampFIP(timer, version); - } - outputRegionFluidInPlace(originalFluidInPlace_.totals, - currentFluidInPlace.totals, - eclState().getUnits(), - 0); - for (size_t reg = 0; reg < originalFluidInPlace_.data.size(); ++reg) { - outputRegionFluidInPlace(originalFluidInPlace_.data[reg], - currentFluidInPlace.data[reg], - eclState().getUnits(), - reg+1); - } - } - - void outputTimestampFIP(const SimulatorTimer& timer, const std::string version) { std::ostringstream ss; @@ -501,50 +343,6 @@ protected: OpmLog::note(ss.str()); } - - void outputRegionFluidInPlace(const std::vector& oip, const std::vector& cip, const UnitSystem& units, const int reg) - { - std::ostringstream ss; - if (!reg) { - ss << " ===================================================\n" - << " : Field Totals :\n"; - } else { - ss << " ===================================================\n" - << " : FIPNUM report region " - << std::setw(2) << reg << " :\n"; - } - if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) { - ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n" - << std::fixed << std::setprecision(0) - << " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n"; - if (!reg) { - ss << " : Pressure is weighted by hydrocarbon pore volume :\n" - << " : Porv volumes are taken at reference conditions :\n"; - } - ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n"; - } - if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) { - ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :\n" - << std::fixed << std::setprecision(0) - << " : PORV =" << std::setw(14) << cip[5] << " RB :\n"; - if (!reg) { - ss << " : Pressure is weighted by hydrocarbon pore volume :\n" - << " : Pore volumes are taken at reference conditions :\n"; - } - ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n"; - } - ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n" - << ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n" - << ":Currently in place :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) << ":" - << std::setw(13) << cip[0] << " :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) << ":\n" - << ":------------------------:------------------------------------------:----------------:------------------------------------------:\n" - << ":Originally in place :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) << ":" - << std::setw(13) << oip[0] << " :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) << ":\n" - << ":========================:==========================================:================:==========================================:\n"; - OpmLog::note(ss.str()); - } - - const EclipseState& eclState() const { return ebosSimulator_.gridManager().eclState(); } @@ -556,9 +354,6 @@ protected: // Data. Simulator& ebosSimulator_; - std::vector fipnum_; - FluidInPlace originalFluidInPlace_; - typedef typename Solver::SolverParameters SolverParameters; SimulatorReport failureReport_;