From 137ff53ae79b2d5380f6c64ba7c17f98b31aee82 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 12 Jan 2018 15:32:04 +0100 Subject: [PATCH 1/8] Remove output of FIP --- .../SimulatorFullyImplicitBlackoilEbos.hpp | 22 ++++++++++++------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index ad4094ac2..5dac74df8 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -209,7 +209,7 @@ public: perfTimer.start(); if (terminal_output_) { - outputFluidInPlace(timer, originalFluidInPlace_); + //outputFluidInPlace(timer, originalFluidInPlace_); } // No per cell data is written for initial step, but will be @@ -293,21 +293,27 @@ public: if (terminal_output_ ) { - outputFluidInPlace(timer, currentFluidInPlace); - - std::string msg = - "Time step took " + std::to_string(solver_timer.secsSinceStart()) + " seconds; " - "total solver time " + std::to_string(report.solver_time) + " seconds."; - OpmLog::debug(msg); + if (!timer.initialStep()) { + const std::string version = moduleVersionName(); + outputTimestampFIP(timer, version); + } } // write simulation state at the report stage Dune::Timer perfTimer; perfTimer.start(); - const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0; + const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0; output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model(), false, nextstep, report); report.output_write_time += perfTimer.stop(); + if (terminal_output_ ) + { + std::string msg = + "Time step took " + std::to_string(solver_timer.secsSinceStart()) + " seconds; " + "total solver time " + std::to_string(report.solver_time) + " seconds."; + OpmLog::debug(msg); + } + } // Stop timer and create timing report From 386ade39f447ec89fa000a7c56018d1f0d8d0d1b Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Mon, 15 Jan 2018 10:14:13 +0100 Subject: [PATCH 2/8] flow: let core ebos handle the output directory for the result files this adds a new parameter --ecl-output-dir=$FOO to the "virtual" command line arguments of the ebos simulator. --- opm/autodiff/FlowMainEbos.hpp | 30 ++++++++++-------------------- 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index 1f3f4f3bd..0cbf92cda 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -26,9 +26,7 @@ #include - #include -#include #include #include @@ -257,23 +255,7 @@ namespace Opm << std::endl; } - output_to_files_ = output_cout_ && output_ > OUTPUT_NONE; - - // Setup output directory. - auto& ioConfig = eclState().getIOConfig(); - // Default output directory is the directory where the deck is found. - const std::string default_output_dir = ioConfig.getOutputDir(); - output_dir_ = param_.getDefault("output_dir", default_output_dir); - // Override output directory if user specified. - ioConfig.setOutputDir(output_dir_); - - // Write parameters used for later reference. (only if rank is zero) - if (output_to_files_) { - // Create output directory if needed. - ensureDirectoryExists(output_dir_); - // Write simulation parameters. - param_.writeParam(output_dir_ + "/simulation.param"); - } + output_to_files_ = output_cout_ && (output_ != OUTPUT_NONE); } // Setup OpmLog backend with output_dir. @@ -413,9 +395,17 @@ namespace Opm argv.push_back("flow_ebos"); std::string deckFileParam("--ecl-deck-file-name="); - deckFileParam += param_.get("deck_filename"); + const std::string& deckFileName = param_.get("deck_filename"); + deckFileParam += deckFileName; argv.push_back(deckFileParam.c_str()); + const std::string default_output_dir = boost::filesystem::basename(deckFileName); + output_dir_ = param_.getDefault("output_dir", default_output_dir); + + std::string outputDirParam("--ecl-output-dir="); + outputDirParam += output_dir_; + argv.push_back(outputDirParam.c_str()); + #if defined(_OPENMP) std::string numThreadsParam("--threads-per-process="); int numThreads = omp_get_max_threads(); From 36f6b7ad00f8533cc46059dfc3128afa3b442e52 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 15 Jan 2018 13:06:42 +0100 Subject: [PATCH 3/8] Remove writInit() The EGRID and INIT files are written using ebos --- opm/autodiff/FlowMainEbos.hpp | 259 ---------------------------------- 1 file changed, 259 deletions(-) diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index 0cbf92cda..14f418157 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -108,7 +108,6 @@ namespace Opm printPRTHeader(); extractMessages(); runDiagnostics(); - writeInit(); setupOutputWriter(); setupLinearSolver(); createSimulator(); @@ -420,10 +419,6 @@ namespace Opm ebosSimulator_.reset(new EbosSimulator(/*verbose=*/false)); ebosSimulator_->model().applyInitialSolution(); - // Create a grid with a global view. - globalGrid_.reset(new Grid(grid())); - globalGrid_->switchToGlobalView(); - try { if (output_cout_) { MissingFeatures::checkKeywords(deck()); @@ -464,9 +459,6 @@ namespace Opm const Schedule& schedule() const { return ebosSimulator_->gridManager().schedule(); } - const SummaryConfig& summaryConfig() const - { return ebosSimulator_->gridManager().summaryConfig(); } - // Extract messages from parser. // Writes to: // OpmLog singleton. @@ -513,27 +505,6 @@ namespace Opm diagnostic.diagnosis(eclState(), deck(), this->grid()); } - void writeInit() - { - bool output = ( output_ > OUTPUT_LOG_ONLY ); - bool output_ecl = param_.getDefault("output_ecl", true); - auto int_vectors = computeCellRanks(output, output_ecl); - - if( output && output_ecl && grid().comm().rank() == 0 ) - { - exportNncStructure_(); - - const EclipseGrid& inputGrid = eclState().getInputGrid(); - eclIO_.reset(new EclipseIO(eclState(), - UgGridHelpers::createEclipseGrid( this->globalGrid() , inputGrid ), - schedule(), - summaryConfig())); - eclIO_->writeInitial(computeLegacySimProps_(), int_vectors, nnc_); - Problem& problem = ebosProblem(); - problem.setEclIO(std::move(eclIO_)); - } - } - // Setup output writer. // Writes to: // output_writer_ @@ -679,232 +650,6 @@ namespace Opm Grid& grid() { return ebosSimulator_->gridManager().grid(); } - const Grid& globalGrid() - { return *globalGrid_; } - - Problem& ebosProblem() - { return ebosSimulator_->problem(); } - - const Problem& ebosProblem() const - { return ebosSimulator_->problem(); } - - std::shared_ptr materialLawManager() - { return ebosProblem().materialLawManager(); } - - Scalar gravity() const - { return ebosProblem().gravity()[2]; } - - std::map > computeCellRanks(bool output, bool output_ecl) - { - std::map > integerVectors; - - if( output && output_ecl && grid().comm().size() > 1 ) - { - typedef typename Grid::LeafGridView GridView; -#if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6) - using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; -#else - // Get the owner rank number for each cell - using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; -#endif - using Handle = CellOwnerDataHandle; - const Grid& globalGrid = this->globalGrid(); - const auto& globalGridView = globalGrid.leafGridView(); -#if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6) - ElementMapper globalMapper(globalGridView, Dune::mcmgElementLayout()); -#else - ElementMapper globalMapper(globalGridView); -#endif - const auto globalSize = globalGrid.size(0); - std::vector ranks(globalSize, -1); - Handle handle(globalMapper, ranks); - this->grid().gatherData(handle); - integerVectors.emplace("MPI_RANK", ranks); - } - - return integerVectors; - } - - data::Solution computeLegacySimProps_() - { - const int* dims = UgGridHelpers::cartDims(grid()); - const int globalSize = dims[0]*dims[1]*dims[2]; - - data::CellData tranx = {UnitSystem::measure::transmissibility, std::vector( globalSize ), data::TargetType::INIT}; - data::CellData trany = {UnitSystem::measure::transmissibility, std::vector( globalSize ), data::TargetType::INIT}; - data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector( globalSize ), data::TargetType::INIT}; - - for (size_t i = 0; i < tranx.data.size(); ++i) { - tranx.data[0] = 0.0; - trany.data[0] = 0.0; - tranz.data[0] = 0.0; - } - - const Grid& globalGrid = this->globalGrid(); - const auto& globalGridView = globalGrid.leafGridView(); - typedef typename Grid::LeafGridView GridView; -#if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6) - typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; - ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); -#else - typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; - ElementMapper globalElemMapper(globalGridView); -#endif - const auto& cartesianCellIdx = globalGrid.globalCell(); - - const auto* globalTrans = &(ebosSimulator_->gridManager().globalTransmissibility()); - if (grid().comm().size() < 2) { - // in the sequential case we must use the transmissibilites defined by - // the problem. (because in the sequential case, the grid manager does - // not compute "global" transmissibilities for performance reasons. in - // the parallel case, the problem's transmissibilities can't be used - // because this object refers to the distributed grid and we need the - // sequential version here.) - globalTrans = &ebosSimulator_->problem().eclTransmissibilities(); - } - - auto elemIt = globalGridView.template begin(); - const auto& elemEndIt = globalGridView.template end(); - for (; elemIt != elemEndIt; ++ elemIt) { - const auto& elem = *elemIt; - - auto isIt = globalGridView.ibegin(elem); - const auto& isEndIt = globalGridView.iend(elem); - for (; isIt != isEndIt; ++ isIt) { - const auto& is = *isIt; - - if (!is.neighbor()) - { - continue; // intersection is on the domain boundary - } - - unsigned c1 = globalElemMapper.index(is.inside()); - unsigned c2 = globalElemMapper.index(is.outside()); - - if (c1 > c2) - { - continue; // we only need to handle each connection once, thank you. - } - - - int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]); - int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]); - if (gc2 - gc1 == 1) { - tranx.data[gc1] = globalTrans->transmissibility(c1, c2); - } - - if (gc2 - gc1 == dims[0]) { - trany.data[gc1] = globalTrans->transmissibility(c1, c2); - } - - if (gc2 - gc1 == dims[0]*dims[1]) { - tranz.data[gc1] = globalTrans->transmissibility(c1, c2); - } - } - } - - return {{"TRANX" , tranx}, - {"TRANY" , trany} , - {"TRANZ" , tranz}}; - } - - void exportNncStructure_() - { - nnc_ = eclState().getInputNNC(); - int nx = eclState().getInputGrid().getNX(); - int ny = eclState().getInputGrid().getNY(); - //int nz = eclState().getInputGrid().getNZ() - - const Grid& globalGrid = this->globalGrid(); - const auto& globalGridView = globalGrid.leafGridView(); - typedef typename Grid::LeafGridView GridView; -#if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6) - typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; - ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); -#else - typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; - ElementMapper globalElemMapper(globalGridView); -#endif - - const auto* globalTrans = &(ebosSimulator_->gridManager().globalTransmissibility()); - if (grid().comm().size() < 2) { - // in the sequential case we must use the transmissibilites defined by - // the problem. (because in the sequential case, the grid manager does - // not compute "global" transmissibilities for performance reasons. in - // the parallel case, the problem's transmissibilities can't be used - // because this object refers to the distributed grid and we need the - // sequential version here.) - globalTrans = &ebosSimulator_->problem().eclTransmissibilities(); - } - - auto elemIt = globalGridView.template begin(); - const auto& elemEndIt = globalGridView.template end(); - for (; elemIt != elemEndIt; ++ elemIt) { - const auto& elem = *elemIt; - - auto isIt = globalGridView.ibegin(elem); - const auto& isEndIt = globalGridView.iend(elem); - for (; isIt != isEndIt; ++ isIt) { - const auto& is = *isIt; - - if (!is.neighbor()) - { - continue; // intersection is on the domain boundary - } - - unsigned c1 = globalElemMapper.index(is.inside()); - unsigned c2 = globalElemMapper.index(is.outside()); - - if (c1 > c2) - { - continue; // we only need to handle each connection once, thank you. - } - - // TODO (?): use the cartesian index mapper to make this code work - // with grids other than Dune::CpGrid. The problem is that we need - // the a mapper for the sequential grid, not for the distributed one. - int cc1 = globalGrid.globalCell()[c1]; - int cc2 = globalGrid.globalCell()[c2]; - - if (std::abs(cc1 - cc2) != 1 && - std::abs(cc1 - cc2) != nx && - std::abs(cc1 - cc2) != nx*ny) - { - nnc_.addNNC(cc1, cc2, globalTrans->transmissibility(c1, c2)); - } - } - } - } - - /// Convert saturations from a vector of individual phase saturation vectors - /// to an interleaved format where all values for a given cell come before all - /// values for the next cell, all in a single vector. - template - void convertSats(std::vector& sat_interleaved, const std::vector< std::vector >& sat, const PhaseUsage& pu) - { - assert(sat.size() == 3); - const auto nc = sat[0].size(); - const auto np = sat_interleaved.size() / nc; - for (size_t c = 0; c < nc; ++c) { - if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const int opos = pu.phase_pos[BlackoilPhases::Liquid]; - const std::vector& sat_p = sat[ FluidSystem::oilPhaseIdx]; - sat_interleaved[np*c + opos] = sat_p[c]; - } - if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; - const std::vector& sat_p = sat[ FluidSystem::waterPhaseIdx]; - sat_interleaved[np*c + wpos] = sat_p[c]; - } - if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; - const std::vector& sat_p = sat[ FluidSystem::gasPhaseIdx]; - sat_interleaved[np*c + gpos] = sat_p[c]; - } - } - } - - std::unique_ptr ebosSimulator_; int mpi_rank_ = 0; bool output_cout_ = false; @@ -913,15 +658,11 @@ namespace Opm ParameterGroup param_; bool output_to_files_ = false; std::string output_dir_ = std::string("."); - NNC nnc_; - std::unique_ptr eclIO_; std::unique_ptr output_writer_; boost::any parallel_information_; std::unique_ptr fis_solver_; std::unique_ptr simulator_; std::string logFile_; - // Needs to be shared pointer because it gets initialzed before MPI_Init. - std::shared_ptr globalGrid_; }; } // namespace Opm From 6d0c716d76010d1caca489f3295c7e73ed748661 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 15 Jan 2018 13:51:40 +0100 Subject: [PATCH 4/8] Clean-up FIP --- opm/autodiff/BlackoilModelEbos.hpp | 157 +------------ opm/autodiff/BlackoilOutputEbos.hpp | 147 +------------ .../SimulatorFullyImplicitBlackoilEbos.hpp | 207 +----------------- 3 files changed, 9 insertions(+), 502 deletions(-) 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_; From f980c39e6c6b59fd6da9d262b3d024a879b3d837 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 16 Jan 2018 15:50:21 +0100 Subject: [PATCH 5/8] Change directory of parallel restart test --- compareECLFiles.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake index 8405ac70c..5ddd6b44f 100644 --- a/compareECLFiles.cmake +++ b/compareECLFiles.cmake @@ -119,7 +119,7 @@ function(add_test_compare_parallel_restarted_simulation) set(multiValueArgs TEST_ARGS) cmake_parse_arguments(PARAM "$" "${oneValueArgs}" "${multiValueArgs}" ${ARGN} ) - set(RESULT_PATH ${BASE_RESULT_PATH}/restart/${PARAM_SIMULATOR}+${PARAM_CASENAME}) + set(RESULT_PATH ${BASE_RESULT_PATH}/parallelRestart/${PARAM_SIMULATOR}+${PARAM_CASENAME}) set(TEST_ARGS ${OPM_DATA_ROOT}/${PARAM_CASENAME}/${PARAM_FILENAME} ${PARAM_TEST_ARGS}) opm_add_test(compareParallelRestartedSim_${PARAM_SIMULATOR}+${PARAM_FILENAME} NO_COMPILE From b38430ea2ed6fe317d218676ce4dd649cd07da62 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 16 Jan 2018 15:50:43 +0100 Subject: [PATCH 6/8] Output TCPU for every substep. --- opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp index e5947f310..29f8618c7 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp @@ -333,7 +333,7 @@ namespace Opm { perfTimer.start(); bool substep = true; const auto& physicalModel = solver.model(); - outputWriter->writeTimeStep( substepTimer, state, well_state, physicalModel, substep); + outputWriter->writeTimeStep( substepTimer, state, well_state, physicalModel, substep, -1.0, substepReport); report.output_write_time += perfTimer.secsSinceStart(); } From 84a8b8eca6decf6b444fabfe550582171ca02386 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 17 Jan 2018 15:26:06 +0100 Subject: [PATCH 7/8] Pass empty regionSummaryData and blockSummaryData to output --- opm/autodiff/FlowMainEbos.hpp | 5 +++++ opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp | 7 ++----- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index 14f418157..c59e9d632 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -405,6 +405,11 @@ namespace Opm outputDirParam += output_dir_; argv.push_back(outputDirParam.c_str()); + const bool restart_double_si = param_.getDefault("restart_double_si", false); + std::string outputDoublePrecisionParam("--ecl-output-double-precision="); + outputDoublePrecisionParam += restart_double_si; + argv.push_back(outputDoublePrecisionParam.c_str()); + #if defined(_OPENMP) std::string numThreadsParam("--threads-per-process="); int numThreads = omp_get_max_threads(); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp index 3dec4af6a..8d17b3bcd 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp @@ -333,11 +333,6 @@ namespace Opm if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) { std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl; } else { - if ( timer.initialStep() ) - { - // Set the initial OIP - eclIO_->overwriteInitialOIP(simProps); - } // ... insert "extra" data (KR, VISC, ...) const int reportStepForOutput = substep ? timer.reportStepNum() + 1 : timer.reportStepNum(); eclIO_->writeTimeStep(reportStepForOutput, @@ -346,6 +341,8 @@ namespace Opm simProps, wellState.report(phaseUsage_), miscSummaryData, + {}, //regionData + {}, //blockData extraRestartData, restart_double_si_); } From cdeefc3a344d62e7333dac0801f28e9509157c35 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 22 Jan 2018 14:03:15 +0100 Subject: [PATCH 8/8] Fix command line argument for double precision restart output --- opm/autodiff/FlowMainEbos.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index c59e9d632..3410f4673 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -407,7 +407,7 @@ namespace Opm const bool restart_double_si = param_.getDefault("restart_double_si", false); std::string outputDoublePrecisionParam("--ecl-output-double-precision="); - outputDoublePrecisionParam += restart_double_si; + outputDoublePrecisionParam += restart_double_si ? "true" : "false"; argv.push_back(outputDoublePrecisionParam.c_str()); #if defined(_OPENMP)