diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 62aa71d38..96facfa87 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -67,6 +67,7 @@ list (APPEND MAIN_SOURCE_FILES opm/polymer/fullyimplicit/PolymerPropsAd.cpp opm/simulators/SimulatorCompressibleTwophase.cpp opm/simulators/SimulatorIncompTwophase.cpp + opm/simulators/WellSwitchingLogger.cpp ) @@ -88,6 +89,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_solventprops_ad.cpp tests/test_multisegmentwells.cpp # tests/test_thresholdpressure.cpp + tests/test_wellswitchlogger.cpp ) list (APPEND TEST_DATA_FILES @@ -241,7 +243,9 @@ list (APPEND PUBLIC_HEADER_FILES opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp + opm/simulators/ParallelFileMerger.hpp opm/simulators/SimulatorCompressibleTwophase.hpp opm/simulators/SimulatorIncompTwophase.hpp + opm/simulators/WellSwitchingLogger.hpp ) diff --git a/debian/control b/debian/control index b2bd821ca..7f7dc4e12 100644 --- a/debian/control +++ b/debian/control @@ -7,7 +7,7 @@ Build-Depends: build-essential, debhelper (>= 9), libboost-filesystem-dev, libdune-common-dev, libdune-istl-dev, cmake, libtinyxml-dev, bc, libert.ecl-dev, git, zlib1g-dev, libtool, doxygen, texlive-latex-extra, texlive-latex-recommended, ghostscript, - libopm-core-dev, libeigen3-dev (>= 3.2.0), libopm-material-dev, + libopm-core-dev, libopm-material-dev, libopm-parser-dev, libboost-iostreams-dev, libopm-common-dev, libopm-grid-dev, libdune-grid-dev, libug-dev, libopm-output-dev, libtrilinos-zoltan-dev, libopenmpi-dev, mpi-default-bin diff --git a/debian/mk_orig_source b/debian/mk_orig_source new file mode 100755 index 000000000..809b8d1e9 --- /dev/null +++ b/debian/mk_orig_source @@ -0,0 +1,17 @@ +#!/bin/bash + +ORIGDIR=`pwd` +TOPDIR=`dirname $0` +tmp=`mktemp -d` +pushd . +cd $TOPDIR/.. +git archive --prefix opm-simulators-$1/ -o $tmp/opm-simulators.tar.gz $2 +cd $tmp +wget https://github.com/OPM/eigen3/archive/master.tar.gz -O eigen3-master.tar.gz +tar zxvf opm-simulators.tar.gz +cd opm-simulators-$1 +tar zxvf ../eigen3-master.tar.gz +cd .. +tar zcvf $ORIGDIR/opm-simulators-$1.tar.gz opm-simulators-$1/ +popd +rm $tmp -Rf diff --git a/debian/rules b/debian/rules index 20fd372c2..b0d296980 100755 --- a/debian/rules +++ b/debian/rules @@ -21,7 +21,10 @@ override_dh_auto_build: # consider using -DUSE_VERSIONED_DIR=ON if backporting override_dh_auto_configure: - dh_auto_configure --buildsystem=cmake -- -DCMAKE_BUILD_TYPE=Release -DUSE_MPI=1 -DBUILD_SHARED_LIBS=1 -DCMAKE_INSTALL_DOCDIR=share/doc/libopm-simulators1 -DWHOLE_PROG_OPTIM=OFF -DUSE_RUNPATH=OFF + mkdir build-eigen + cd build-eigen && cmake $(CURDIR)/eigen3-master -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=/tmp/install-eigen3 + cd build-eigen && make install + dh_auto_configure --buildsystem=cmake -- -DCMAKE_BUILD_TYPE=Release -DUSE_MPI=1 -DBUILD_SHARED_LIBS=1 -DCMAKE_INSTALL_DOCDIR=share/doc/libopm-simulators1 -DWHOLE_PROG_OPTIM=OFF -DUSE_RUNPATH=OFF -DCMAKE_PREFIX_PATH=/tmp/install-eigen3 override_dh_auto_install: dh_auto_install -- install-html diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 171cba63e..95da59289 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -92,9 +92,20 @@ namespace Opm { struct SimulatorData { SimulatorData(int num_phases); + + enum FipId { + FIP_AQUA = Opm::Phases::Water, + FIP_LIQUID = Opm::Phases::Oil, + FIP_VAPOUR = Opm::Phases::Gas, + FIP_DISSOLVED_GAS = 3, + FIP_VAPORIZED_OIL = 4, + FIP_PV = 5, //< Pore volume + FIP_WEIGHTED_PRESSURE = 6 + }; std::vector rq; - ADB rsSat; - ADB rvSat; + ADB rsSat; // Saturated gas-oil ratio + ADB rvSat; // Saturated oil-gas ratio + std::array fip; }; typedef typename ModelTraits::ReservoirState ReservoirState; diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index d5826f86f..5d31ccd21 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -424,6 +424,7 @@ typedef Eigen::Array fip(5, V::Zero(nc)); for (int phase = 0; phase < maxnp; ++phase) { if (active_[ phase ]) { const int pos = pu.phase_pos[ phase ]; const auto& b = asImpl().fluidReciprocFVF(phase, canonical_phase_pressures[phase], temperature, rs, rv, cond); - fip[phase] = ((pv_mult * b * saturation[pos] * pv).value()); + sd_.fip[phase] = ((pv_mult * b * saturation[pos] * pv).value()); } } if (active_[ Oil ] && active_[ Gas ]) { // Account for gas dissolved in oil and vaporized oil - const int po = pu.phase_pos[Oil]; - const int pg = pu.phase_pos[Gas]; - fip[3] = rs.value() * fip[po]; - fip[4] = rv.value() * fip[pg]; + sd_.fip[SimulatorData::FIP_DISSOLVED_GAS] = rs.value() * sd_.fip[SimulatorData::FIP_LIQUID]; + sd_.fip[SimulatorData::FIP_VAPORIZED_OIL] = rv.value() * sd_.fip[SimulatorData::FIP_VAPOUR]; } // For a parallel run this is just a local maximum and needs to be updated later @@ -2191,23 +2196,57 @@ typedef Eigen::Array 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/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 304b24203..c704072a9 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include @@ -85,6 +86,7 @@ #include #include + #ifdef _OPENMP #include #endif @@ -150,12 +152,23 @@ namespace Opm asImpl().setupOutputWriter(); asImpl().setupLinearSolver(); asImpl().createSimulator(); - + // Run. - return asImpl().runSimulator(); + auto ret = asImpl().runSimulator(); + + asImpl().mergeParallelLogFiles(); + + return ret; } catch (const std::exception &e) { - std::cerr << "Program threw an exception: " << e.what() << "\n"; + std::ostringstream message; + message << "Program threw an exception: " << e.what(); + + if( output_cout_ ) + { + OpmLog::error(message.str()); + } + return EXIT_FAILURE; } @@ -183,6 +196,7 @@ namespace Opm // members first occur. // setupParallelism() + int mpi_rank_ = 0; bool output_cout_ = false; bool must_distribute_ = false; // setupParameters() @@ -232,9 +246,9 @@ namespace Opm // For a build without MPI the Dune::FakeMPIHelper is used, so rank will // be 0 and size 1. const Dune::MPIHelper& mpi_helper = Dune::MPIHelper::instance(argc, argv); - const int mpi_rank = mpi_helper.rank(); + mpi_rank_ = mpi_helper.rank(); const int mpi_size = mpi_helper.size(); - output_cout_ = ( mpi_rank == 0 ); + output_cout_ = ( mpi_rank_ == 0 ); must_distribute_ = ( mpi_size > 1 ); #ifdef _OPENMP @@ -253,7 +267,7 @@ namespace Opm if (mpi_size == 1) { std::cout << "OpenMP using " << num_omp_threads << " threads." << std::endl; } else { - std::cout << "OpenMP using " << num_omp_threads << " threads on MPI rank " << mpi_rank << "." << std::endl; + std::cout << "OpenMP using " << num_omp_threads << " threads on MPI rank " << mpi_rank_ << "." << std::endl; } } #endif @@ -340,10 +354,14 @@ namespace Opm { // Write parameters used for later reference. (only if rank is zero) output_to_files_ = output_cout_ && param_.getDefault("output", true); + // Always read output_dir as it will be set unconditionally later. + // Not doing this might cause files to be created in the current + // directory. + output_dir_ = + param_.getDefault("output_dir", std::string(".")); + if (output_to_files_) { // Create output directory if needed. - output_dir_ = - param_.getDefault("output_dir", std::string(".")); boost::filesystem::path fpath(output_dir_); if (!is_directory(fpath)) { try { @@ -370,36 +388,80 @@ namespace Opm using boost::filesystem::path; path fpath(deck_filename); std::string baseName; - std::string debugFile; + std::ostringstream debugFileStream; + std::ostringstream logFileStream; + if (boost::to_upper_copy(path(fpath.extension()).string()) == ".DATA") { baseName = path(fpath.stem()).string(); } else { baseName = path(fpath.filename()).string(); } if (param_.has("output_dir")) { - logFile_ = output_dir_ + "/" + baseName + ".PRT"; - debugFile = output_dir_ + "/." + baseName + ".DEBUG"; - } else { - logFile_ = baseName + ".PRT"; - debugFile = "." + baseName + ".DEBUG"; + logFileStream << output_dir_ << "/"; + debugFileStream << output_dir_ + "/"; } - std::shared_ptr prtLog = std::make_shared(logFile_ , Log::NoDebugMessageTypes); + + logFileStream << baseName; + debugFileStream << "." << baseName; + + if ( must_distribute_ && mpi_rank_ != 0 ) + { + // Added rank to log file for non-zero ranks. + // This prevents message loss. + debugFileStream << "."<< mpi_rank_; + // If the following file appears then there is a bug. + logFileStream << "." << mpi_rank_; + } + logFileStream << ".PRT"; + debugFileStream << ".DEBUG"; + + std::string debugFile = debugFileStream.str(); + logFile_ = logFileStream.str(); + + std::shared_ptr prtLog = std::make_shared(logFile_ , Log::NoDebugMessageTypes, false, output_cout_); std::shared_ptr streamLog = std::make_shared(std::cout, Log::StdoutMessageTypes); OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog ); OpmLog::addBackend( "STREAMLOG", streamLog); - std::shared_ptr debugLog = std::make_shared(debugFile, Log::DefaultMessageTypes); + std::shared_ptr debugLog = std::make_shared(debugFile, Log::DefaultMessageTypes, false, output_cout_); OpmLog::addBackend( "DEBUGLOG" , debugLog); prtLog->setMessageFormatter(std::make_shared(false)); streamLog->setMessageLimiter(std::make_shared(10)); streamLog->setMessageFormatter(std::make_shared(true)); + // Read parameters. - OpmLog::debug("\n--------------- Reading parameters ---------------\n"); + if ( output_cout_ ) + { + OpmLog::debug("\n--------------- Reading parameters ---------------\n"); + } } + void mergeParallelLogFiles() + { + // force closing of all log files. + OpmLog::removeAllBackends(); + + if( mpi_rank_ != 0 || !must_distribute_ ) + { + return; + } + + namespace fs = boost::filesystem; + fs::path output_path("."); + if ( param_.has("output_dir") ) + { + output_path = fs::path(output_dir_); + } + + fs::path deck_filename(param_.get("deck_filename")); + + std::for_each(fs::directory_iterator(output_path), + fs::directory_iterator(), + detail::ParallelFileMerger(output_path, deck_filename.stem().string())); + } // Parser the input and creates the Deck and EclipseState objects. // Writes to: @@ -418,7 +480,12 @@ namespace Opm ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }}); deck_ = parser->parseFile(deck_filename, parseContext); checkDeck(deck_, parser); - MissingFeatures::checkKeywords(*deck_); + + if ( output_cout_) + { + MissingFeatures::checkKeywords(*deck_); + } + eclipse_state_.reset(new EclipseState(*deck_, parseContext)); auto ioConfig = eclipse_state_->getIOConfig(); ioConfig->setOutputDir(output_dir_); @@ -619,6 +686,11 @@ namespace Opm // OpmLog singleton. void extractMessages() { + if ( !output_cout_ ) + { + return; + } + auto extractMessage = [](const Message& msg) { auto log_type = detail::convertMessageType(msg.mtype); const auto& location = msg.location; @@ -649,6 +721,11 @@ namespace Opm // OpmLog singleton. void runDiagnostics() { + if( ! output_cout_ ) + { + return; + } + // Run relperm diagnostics RelpermDiagnostics diagnostic; diagnostic.diagnosis(eclipse_state_, deck_, grid_init_->grid()); diff --git a/opm/autodiff/MultisegmentWells.hpp b/opm/autodiff/MultisegmentWells.hpp index 0d8e5aece..144428962 100644 --- a/opm/autodiff/MultisegmentWells.hpp +++ b/opm/autodiff/MultisegmentWells.hpp @@ -22,6 +22,8 @@ #ifndef OPM_MULTISEGMENTWELLS_HEADER_INCLUDED #define OPM_MULTISEGMENTWELLS_HEADER_INCLUDED +#include + #include #include #include @@ -41,6 +43,7 @@ #include #include +#include @@ -78,6 +81,9 @@ namespace Opm { Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>; + using Communication = + Dune::CollectiveCommunication; // --------- Public methods --------- // TODO: using a vector of WellMultiSegmentConstPtr for now diff --git a/opm/autodiff/MultisegmentWells_impl.hpp b/opm/autodiff/MultisegmentWells_impl.hpp index 65c732beb..51262420e 100644 --- a/opm/autodiff/MultisegmentWells_impl.hpp +++ b/opm/autodiff/MultisegmentWells_impl.hpp @@ -824,9 +824,10 @@ namespace Opm MultisegmentWells:: updateWellControls(WellState& xw) const { + wellhelpers::WellSwitchingLogger logger; + if( msWells().empty() ) return ; - std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; // Find, for each well, if any constraints are broken. If so, // switch control to first broken constraint. const int np = numPhases(); @@ -860,9 +861,9 @@ namespace Opm if (ctrl_index != nwc) { // Constraint number ctrl_index was broken, switch to it. // Each well is only active on one process. Therefore we always print the sitch info. - std::cout << "Switching control mode for well " << msWells()[w]->name() - << " from " << modestring[well_controls_iget_type(wc, current)] - << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; + logger.wellSwitched(msWells()[w]->name(), + well_controls_iget_type(wc, current), + well_controls_iget_type(wc, ctrl_index)); xw.currentControls()[w] = ctrl_index; current = xw.currentControls()[w]; } diff --git a/opm/autodiff/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index cd0186c6d..06588e8fd 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -112,7 +112,12 @@ namespace Opm std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; if( param_.getDefault("timestep.adaptive", true ) ) { - adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) ); + + if (param_.getDefault("use_TUNING", false)) { + adaptiveTimeStepping.reset( new AdaptiveTimeStepping( schedule->getTuning(), timer.currentStepNum(), param_, terminal_output_ ) ); + } else { + adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) ); + } } std::string restorefilename = param_.getDefault("restorefile", std::string("") ); @@ -276,11 +281,14 @@ namespace Opm FIPUnitConvert(eclipse_state_->getUnits(), COIP); V OOIP_totals = FIPTotals(OOIP, state); V COIP_totals = FIPTotals(COIP, state); - outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0); - for (size_t reg = 0; reg < OOIP.size(); ++reg) { - outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1); - } + if ( terminal_output_ ) + { + outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0); + for (size_t reg = 0; reg < OOIP.size(); ++reg) { + outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1); + } + } // accumulate total time stime += st; @@ -299,7 +307,11 @@ namespace Opm step_report.total_newton_iterations = solver->nonlinearIterations(); step_report.total_linear_iterations = solver->linearIterations(); step_report.total_linearizations = solver->linearizations(); - step_report.reportParam(tstep_os); + + if ( output_writer_.isIORank() ) + { + step_report.reportParam(tstep_os); + } } // Increment timer, remember well state. diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index 68a2ddef8..2efb79eaf 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -273,6 +273,12 @@ namespace Opm /** \brief return true if output is enabled */ bool output () const { return output_; } + /** \brief Whether this process does write to disk */ + bool isIORank () const + { + parallelOutput_->isIORank(); + } + void restore(SimulatorTimerInterface& timer, BlackoilState& state, WellStateFullyImplicitBlackoil& wellState, @@ -423,33 +429,44 @@ namespace Opm namespace detail { /** - * Converts an ADB into a standard vector by copy + * Converts an ADB::V into a standard vector by copy */ - inline std::vector adbToDoubleVector(const Opm::AutoDiffBlock& adb) { - const auto& adb_v = adb.value(); + inline std::vector adbVToDoubleVector(const Opm::AutoDiffBlock::V& adb_v) { std::vector vec(adb_v.data(), adb_v.data() + adb_v.size()); return vec; } + /** + * Converts an ADB into a standard vector by copy + */ + inline std::vector adbToDoubleVector(const Opm::AutoDiffBlock& adb) { + return adbVToDoubleVector(adb.value()); + } + + + /** + * Returns the data requested in the restartConfig + */ template - std::vector getCellData( + void getRestartData( + std::vector& output, const Opm::PhaseUsage& phaseUsage, - const Model& model, + const Model& physicalModel, const RestartConfig& restartConfig, - const int reportStepNum) { + const int reportStepNum, + const bool log) { + typedef Opm::AutoDiffBlock ADB; - std::vector simProps; + const typename Model::SimulatorData& sd = physicalModel.getSimulatorData(); - //Get the value of each of the keys - std::map outKeywords = restartConfig.getRestartKeywords(reportStepNum); - for (auto& keyValue : outKeywords) { + //Get the value of each of the keys for the restart keywords + std::map rstKeywords = restartConfig.getRestartKeywords(reportStepNum); + for (auto& keyValue : rstKeywords) { keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum); } - const typename Model::SimulatorData& sd = model.getSimulatorData(); - //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]; @@ -463,161 +480,305 @@ namespace Opm /** * Formation volume factors for water, oil, gas */ - if (aqua_active && outKeywords["BW"] > 0) { - outKeywords["BW"] = 0; - simProps.emplace_back(data::CellData{ + if (aqua_active && rstKeywords["BW"] > 0) { + rstKeywords["BW"] = 0; + output.emplace_back(data::CellData{ "1OVERBW", Opm::UnitSystem::measure::water_inverse_formation_volume_factor, - std::move(adbToDoubleVector(sd.rq[aqua_idx].b))}); + adbToDoubleVector(sd.rq[aqua_idx].b), + true}); } - if (liquid_active && outKeywords["BO"] > 0) { - outKeywords["BO"] = 0; - simProps.emplace_back(data::CellData{ + if (liquid_active && rstKeywords["BO"] > 0) { + rstKeywords["BO"] = 0; + output.emplace_back(data::CellData{ "1OVERBO", Opm::UnitSystem::measure::oil_inverse_formation_volume_factor, - std::move(adbToDoubleVector(sd.rq[liquid_idx].b))}); + adbToDoubleVector(sd.rq[liquid_idx].b), + true}); } - if (vapour_active && outKeywords["BG"] > 0) { - outKeywords["BG"] = 0; - simProps.emplace_back(data::CellData{ + if (vapour_active && rstKeywords["BG"] > 0) { + rstKeywords["BG"] = 0; + output.emplace_back(data::CellData{ "1OVERBG", Opm::UnitSystem::measure::gas_inverse_formation_volume_factor, - std::move(adbToDoubleVector(sd.rq[vapour_idx].b))}); + adbToDoubleVector(sd.rq[vapour_idx].b), + true}); } /** * Densities for water, oil gas */ - if (outKeywords["DEN"] > 0) { - outKeywords["DEN"] = 0; + if (rstKeywords["DEN"] > 0) { + rstKeywords["DEN"] = 0; if (aqua_active) { - simProps.emplace_back(data::CellData{ + output.emplace_back(data::CellData{ "WAT_DEN", Opm::UnitSystem::measure::density, - std::move(adbToDoubleVector(sd.rq[aqua_idx].rho))}); + adbToDoubleVector(sd.rq[aqua_idx].rho), + true}); } if (liquid_active) { - simProps.emplace_back(data::CellData{ + output.emplace_back(data::CellData{ "OIL_DEN", Opm::UnitSystem::measure::density, - std::move(adbToDoubleVector(sd.rq[liquid_idx].rho))}); + adbToDoubleVector(sd.rq[liquid_idx].rho), + true}); } if (vapour_active) { - simProps.emplace_back(data::CellData{ + output.emplace_back(data::CellData{ "GAS_DEN", Opm::UnitSystem::measure::density, - std::move(adbToDoubleVector(sd.rq[vapour_idx].rho))}); + adbToDoubleVector(sd.rq[vapour_idx].rho), + true}); } } /** * Viscosities for water, oil gas */ - if (outKeywords["VISC"] > 0) { - outKeywords["VISC"] = 0; + if (rstKeywords["VISC"] > 0) { + rstKeywords["VISC"] = 0; if (aqua_active) { - simProps.emplace_back(data::CellData{ + output.emplace_back(data::CellData{ "WAT_VISC", Opm::UnitSystem::measure::viscosity, - std::move(adbToDoubleVector(sd.rq[aqua_idx].mu))}); + adbToDoubleVector(sd.rq[aqua_idx].mu), + true}); } if (liquid_active) { - simProps.emplace_back(data::CellData{ + output.emplace_back(data::CellData{ "OIL_VISC", Opm::UnitSystem::measure::viscosity, - std::move(adbToDoubleVector(sd.rq[liquid_idx].mu))}); + adbToDoubleVector(sd.rq[liquid_idx].mu), + true}); } if (vapour_active) { - simProps.emplace_back(data::CellData{ + output.emplace_back(data::CellData{ "GAS_VISC", Opm::UnitSystem::measure::viscosity, - std::move(adbToDoubleVector(sd.rq[vapour_idx].mu))}); + adbToDoubleVector(sd.rq[vapour_idx].mu), + true}); } } /** * Relative permeabilities for water, oil, gas */ - if (aqua_active && outKeywords["KRW"] > 0) { + if (aqua_active && rstKeywords["KRW"] > 0) { if (sd.rq[aqua_idx].kr.size() > 0) { - outKeywords["KRW"] = 0; - simProps.emplace_back(data::CellData{ + rstKeywords["KRW"] = 0; + output.emplace_back(data::CellData{ "WATKR", Opm::UnitSystem::measure::permeability, - std::move(adbToDoubleVector(sd.rq[aqua_idx].kr))}); + adbToDoubleVector(sd.rq[aqua_idx].kr), + true}); } else { - Opm::OpmLog::warning("Empty:WATKR", - "Not emitting empty Water Rel-Perm"); + if ( log ) + { + Opm::OpmLog::warning("Empty:WATKR", + "Not emitting empty Water Rel-Perm"); + } } } - if (liquid_active && outKeywords["KRO"] > 0) { + if (liquid_active && rstKeywords["KRO"] > 0) { if (sd.rq[liquid_idx].kr.size() > 0) { - outKeywords["KRO"] = 0; - simProps.emplace_back(data::CellData{ + rstKeywords["KRO"] = 0; + output.emplace_back(data::CellData{ "OILKR", Opm::UnitSystem::measure::permeability, - std::move(adbToDoubleVector(sd.rq[liquid_idx].kr))}); + adbToDoubleVector(sd.rq[liquid_idx].kr), + true}); } else { - Opm::OpmLog::warning("Empty:OILKR", - "Not emitting empty Oil Rel-Perm"); + if ( log ) + { + Opm::OpmLog::warning("Empty:OILKR", + "Not emitting empty Oil Rel-Perm"); + } } } - if (vapour_active && outKeywords["KRG"] > 0) { + if (vapour_active && rstKeywords["KRG"] > 0) { if (sd.rq[vapour_idx].kr.size() > 0) { - outKeywords["KRG"] = 0; - simProps.emplace_back(data::CellData{ + rstKeywords["KRG"] = 0; + output.emplace_back(data::CellData{ "GASKR", Opm::UnitSystem::measure::permeability, - std::move(adbToDoubleVector(sd.rq[vapour_idx].kr))}); + adbToDoubleVector(sd.rq[vapour_idx].kr), + true}); } else { - Opm::OpmLog::warning("Empty:GASKR", - "Not emitting empty Gas Rel-Perm"); + if ( log ) + { + Opm::OpmLog::warning("Empty:GASKR", + "Not emitting empty Gas Rel-Perm"); + } } } /** * Vaporized and dissolved gas/oil ratio */ - if (vapour_active && liquid_active && outKeywords["RSSAT"] > 0) { - outKeywords["RSSAT"] = 0; - simProps.emplace_back(data::CellData{ + if (vapour_active && liquid_active && rstKeywords["RSSAT"] > 0) { + rstKeywords["RSSAT"] = 0; + output.emplace_back(data::CellData{ "RSSAT", Opm::UnitSystem::measure::gas_oil_ratio, - std::move(adbToDoubleVector(sd.rsSat))}); + adbToDoubleVector(sd.rsSat), + true}); } - if (vapour_active && liquid_active && outKeywords["RVSAT"] > 0) { - outKeywords["RVSAT"] = 0; - simProps.emplace_back(data::CellData{ + if (vapour_active && liquid_active && rstKeywords["RVSAT"] > 0) { + rstKeywords["RVSAT"] = 0; + output.emplace_back(data::CellData{ "RVSAT", Opm::UnitSystem::measure::oil_gas_ratio, - std::move(adbToDoubleVector(sd.rvSat))}); + adbToDoubleVector(sd.rvSat), + true}); } /** * Bubble point and dew point pressures */ - if (vapour_active && liquid_active && outKeywords["PBPD"] > 0) { - outKeywords["PBPD"] = 0; + if (log && vapour_active && + liquid_active && rstKeywords["PBPD"] > 0) { + rstKeywords["PBPD"] = 0; Opm::OpmLog::warning("Bubble/dew point pressure output unsupported", "Writing bubble points and dew points (PBPD) to file is unsupported, " "as the simulator does not use these internally."); } //Warn for any unhandled keyword - for (auto& keyValue : outKeywords) { - if (keyValue.second > 0) { - std::string logstring = "Keyword '"; - logstring.append(keyValue.first); - logstring.append("' is unhandled for output to file."); - Opm::OpmLog::warning("Unhandled output keyword", logstring); + if (log) { + for (auto& keyValue : rstKeywords) { + if (keyValue.second > 0) { + std::string logstring = "Keyword '"; + logstring.append(keyValue.first); + logstring.append("' is unhandled for output to file."); + Opm::OpmLog::warning("Unhandled output keyword", logstring); + } } } + } - return simProps; + + + + /** + * 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( + std::vector& output, + const Opm::PhaseUsage& phaseUsage, + const Model& physicalModel, + const SummaryConfig& summaryConfig) { + + typedef Opm::AutoDiffBlock ADB; + + const typename Model::SimulatorData& sd = physicalModel.getSimulatorData(); + + //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.emplace_back(data::CellData{ + "WIP", + Opm::UnitSystem::measure::volume, + adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_AQUA]), + false}); + } + if (liquid_active) { + //Oil in place (liquid phase only) + if (hasFRBKeyword(summaryConfig, "OIPL")) { + output.emplace_back(data::CellData{ + "OIPL", + Opm::UnitSystem::measure::volume, + adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_LIQUID]), + false}); + } + //Oil in place (gas phase only) + if (hasFRBKeyword(summaryConfig, "OIPG")) { + output.emplace_back(data::CellData{ + "OIPG", + Opm::UnitSystem::measure::volume, + adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_VAPORIZED_OIL]), + false}); + } + // Oil in place (in liquid and gas phases) + if (hasFRBKeyword(summaryConfig, "OIP")) { + ADB::V oip = sd.fip[Model::SimulatorData::FIP_LIQUID] + + sd.fip[Model::SimulatorData::FIP_VAPORIZED_OIL]; + output.emplace_back(data::CellData{ + "OIP", + Opm::UnitSystem::measure::volume, + adbVToDoubleVector(oip), + false}); + } + } + if (vapour_active) { + // Gas in place (gas phase only) + if (hasFRBKeyword(summaryConfig, "GIPG")) { + output.emplace_back(data::CellData{ + "GIPG", + Opm::UnitSystem::measure::volume, + adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_VAPOUR]), + false}); + } + // Gas in place (liquid phase only) + if (hasFRBKeyword(summaryConfig, "GIPL")) { + output.emplace_back(data::CellData{ + "GIPL", + Opm::UnitSystem::measure::volume, + adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_DISSOLVED_GAS]), + false}); + } + // Gas in place (in both liquid and gas phases) + if (hasFRBKeyword(summaryConfig, "GIP")) { + ADB::V gip = sd.fip[Model::SimulatorData::FIP_VAPOUR] + + sd.fip[Model::SimulatorData::FIP_DISSOLVED_GAS]; + output.emplace_back(data::CellData{ + "GIP", + Opm::UnitSystem::measure::volume, + adbVToDoubleVector(gip), + false}); + } + } + // Cell pore volume in reservoir conditions + if (hasFRBKeyword(summaryConfig, "RPV")) { + output.emplace_back(data::CellData{ + "RPV", + Opm::UnitSystem::measure::volume, + adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_PV]), + false}); + } + // Pressure averaged value (hydrocarbon pore volume weighted) + if (summaryConfig.hasKeyword("FPRH") || summaryConfig.hasKeyword("RPRH")) { + output.emplace_back(data::CellData{ + "PRH", + Opm::UnitSystem::measure::pressure, + adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_WEIGHTED_PRESSURE]), + false}); + } } } @@ -635,8 +796,15 @@ namespace Opm bool substep) { const RestartConfig& restartConfig = eclipseState_->getRestartConfig(); + const SummaryConfig& summaryConfig = eclipseState_->getSummaryConfig(); const int reportStepNum = timer.reportStepNum(); - std::vector cellData = detail::getCellData( phaseUsage_, physicalModel, restartConfig, reportStepNum ); + bool logMessages = output_ && parallelOutput_->isIORank(); + + std::vector cellData; + detail::getRestartData( cellData, phaseUsage_, physicalModel, + restartConfig, reportStepNum, logMessages ); + detail::getSummaryData( cellData, phaseUsage_, physicalModel, summaryConfig ); + writeTimeStepWithCellProperties(timer, localState, localWellState, cellData, substep); } } diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 7268453ae..52030526f 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -22,6 +22,8 @@ #ifndef OPM_STANDARDWELLS_HEADER_INCLUDED #define OPM_STANDARDWELLS_HEADER_INCLUDED +#include + #include #include @@ -40,7 +42,7 @@ #include #include #include - +#include namespace Opm { @@ -58,7 +60,9 @@ namespace Opm { // --------- Types --------- using ADB = AutoDiffBlock; using Vector = ADB::V; - using V = ADB::V; + using Communication = + Dune::CollectiveCommunication; // copied from BlackoilModelBase // should put to somewhere better diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 579655b28..719aa97fa 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -72,7 +72,7 @@ namespace Opm - StandardWells::StandardWells(const Wells* wells_arg) +StandardWells::StandardWells(const Wells* wells_arg) : wells_active_(wells_arg!=nullptr) , wells_(wells_arg) , wops_(wells_arg) @@ -708,9 +708,10 @@ namespace Opm StandardWells:: updateWellControls(WellState& xw) const { + wellhelpers::WellSwitchingLogger logger; + if( !localWellsActive() ) return ; - std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; // Find, for each well, if any constraints are broken. If so, // switch control to first broken constraint. const int np = wells().number_of_phases; @@ -745,11 +746,9 @@ namespace Opm // Constraint number ctrl_index was broken, switch to it. // We disregard terminal_ouput here as with it only messages // for wells on one process will be printed. - std::ostringstream ss; - ss << "Switching control mode for well " << wells().name[w] - << " from " << modestring[well_controls_iget_type(wc, current)] - << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; - OpmLog::info(ss.str()); + logger.wellSwitched(wells().name[w], + well_controls_iget_type(wc, current), + well_controls_iget_type(wc, ctrl_index)); xw.currentControls()[w] = ctrl_index; current = xw.currentControls()[w]; } diff --git a/opm/autodiff/WellHelpers.hpp b/opm/autodiff/WellHelpers.hpp index 808df2ab6..6e0728c88 100644 --- a/opm/autodiff/WellHelpers.hpp +++ b/opm/autodiff/WellHelpers.hpp @@ -35,6 +35,8 @@ namespace Opm { namespace wellhelpers { + + inline double rateToCompare(const std::vector& well_phase_flow_rate, const int well, 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_; }; diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 68147b5b0..72a1e401d 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -33,6 +33,7 @@ #include #include #include +#include #include #include #include @@ -555,7 +556,7 @@ namespace { const std::vector& sat = state.saturation; const std::vector cond = phaseCondition(); - std::vector pressure = computePressures(state); + std::vector pressure = computePressures(state); const ADB pv_mult = poroMult(press); const V& pv = geo_.poreVolume(); @@ -596,6 +597,18 @@ namespace { } } + // Fluid in place is not implemented in this class. + // See BlackoilModelBase::computeFluidInPlace(...) for how it's implemented there + // FIXME: This following code has not been tested. + if (sd_.fip[0].size() == 0) { + OpmLog::warning("NOT_COMPUTING_FIP", + "Computing fluid in place is not implemented for summary files."); + for (int i = 0; i < 7; ++i) { + sd_.fip[i] = V::Zero(nc); + } + } + + return values; } diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp index 7a5694530..24aad8c1b 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp @@ -80,11 +80,23 @@ namespace Opm { : rq(num_phases) , rsSat(ADB::null()) , rvSat(ADB::null()) + , fip() { } + + enum FipId { + FIP_AQUA = BlackoilPropsAdInterface::Water, + FIP_LIQUID = BlackoilPropsAdInterface::Oil, + FIP_VAPOUR = BlackoilPropsAdInterface::Gas, + FIP_DISSOLVED_GAS = 3, + FIP_VAPORIZED_OIL = 4, + FIP_PV = 5, //< Pore volume + FIP_WEIGHTED_PRESSURE = 6 + }; std::vector rq; ADB rsSat; ADB rvSat; + std::array fip; }; /// Construct a solver. It will retain references to the diff --git a/opm/simulators/ParallelFileMerger.hpp b/opm/simulators/ParallelFileMerger.hpp new file mode 100644 index 000000000..58ddd0d6e --- /dev/null +++ b/opm/simulators/ParallelFileMerger.hpp @@ -0,0 +1,132 @@ +/* + Copyright 2016 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 STATOIL AS. + + 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_PARALLELFILEMERGER_HEADER_INCLUDED +#define OPM_PARALLELFILEMERGER_HEADER_INCLUDED + +#include + +#include +#include +#include + +namespace Opm +{ +namespace detail +{ + +namespace fs = boost::filesystem; + +/// \brief A functor that merges multiple files of a parallel run to one file. +/// +/// Without care multiple processes might log messages in a parallel run. +/// Non-root processes will do that to seperate files +/// .. debugStream_; + /// \brief Stream to *.PRT file + std::unique_ptr logStream_; +}; +} // end namespace detail +} // end namespace OPM +#endif // end header guard diff --git a/opm/simulators/WellSwitchingLogger.cpp b/opm/simulators/WellSwitchingLogger.cpp new file mode 100644 index 000000000..d2e98ab5b --- /dev/null +++ b/opm/simulators/WellSwitchingLogger.cpp @@ -0,0 +1,205 @@ +/* + Copyright 2016 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 Statoil AS + + 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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include + +namespace Opm +{ +namespace wellhelpers +{ + +#if HAVE_MPI +int WellSwitchingLogger::calculateMessageSize(std::vector& well_name_lengths) +{ + + // Each process will send a message to the root process with + // the following data: + // total number of switches, for each switch the length of the + // well name, for each switch the well name and the two controls. + well_name_lengths.reserve(switchMap_.size()); + + for(const auto& switchEntry : switchMap_) + { + int length = switchEntry.first.size() +1; //we write an additional \0 + well_name_lengths.push_back(length); + } + + // compute the message size + int message_size = 0; + int increment = 0; + // number of switches + MPI_Pack_size(1, MPI_INT, MPI_COMM_WORLD, &message_size); + // const char* length include delimiter for each switch + MPI_Pack_size(switchMap_.size(), MPI_INT, MPI_COMM_WORLD, &increment); + message_size += increment; + + // for each well the name + two controls in one write + for(const auto& length : well_name_lengths) + { + // well name + MPI_Pack_size(length, MPI_CHAR, MPI_COMM_WORLD, &increment); + message_size += increment; + // controls + MPI_Pack_size(2, MPI_CHAR, MPI_COMM_WORLD, &increment); + message_size += increment; + } + return message_size; +} + +void WellSwitchingLogger::packData(std::vector& well_name_lengths, + std::vector& buffer) +{ + // Pack the data + // number of switches + int offset = 0; + int no_switches = switchMap_.size(); + MPI_Pack(&no_switches, 1, MPI_INT, buffer.data(), buffer.size(), + &offset, MPI_COMM_WORLD); + MPI_Pack(well_name_lengths.data(), well_name_lengths.size(), + MPI_INT, buffer.data(), buffer.size(), + &offset, MPI_COMM_WORLD); + + for(const auto& switchEntry : switchMap_) + { + // well name + auto& well_name = switchEntry.first; + MPI_Pack(const_cast(well_name.c_str()), well_name.size()+1, + MPI_CHAR, buffer.data(), buffer.size(), + &offset, MPI_COMM_WORLD); + + // controls + MPI_Pack(const_cast(switchEntry.second.data()), 2 , MPI_CHAR, + buffer.data(), buffer.size(), &offset, MPI_COMM_WORLD); + } +} + +void WellSwitchingLogger::unpackDataAndLog(std::vector& recv_buffer, + const std::vector& displ) +{ + for(int p=1; p < cc_.size(); ++p) + { + int offset = displ[p]; + int no_switches = 0; + MPI_Unpack(recv_buffer.data(), recv_buffer.size(), &offset, + &no_switches, 1, MPI_INT, MPI_COMM_WORLD); + + if ( no_switches == 0 ) + { + continue; + } + + std::vector well_name_lengths(no_switches); + + MPI_Unpack(recv_buffer.data(), recv_buffer.size(), &offset, + well_name_lengths.data(), well_name_lengths.size(), + MPI_INT, MPI_COMM_WORLD); + + std::vector well_name; + for ( int i = 0; i < no_switches; ++i ) + { + well_name.resize(well_name_lengths[i]); + MPI_Unpack(recv_buffer.data(), recv_buffer.size(), &offset, + well_name.data(), well_name_lengths[i], MPI_CHAR, + MPI_COMM_WORLD); + + std::array fromto{{}}; + MPI_Unpack(recv_buffer.data(), recv_buffer.size(), &offset, + fromto.data(), 2, MPI_CHAR, MPI_COMM_WORLD); + + logSwitch(well_name.data(), fromto, p); + } + } +} + +void WellSwitchingLogger::logSwitch(const char* name, std::array fromto, + int rank) +{ + std::ostringstream ss; + ss << "Switching control mode for well " << name + << " from " << modestring[WellControlType(fromto[0])] + << " to " << modestring[WellControlType(fromto[1])] + << " on rank " << rank << std::endl; + OpmLog::info(ss.str()); +} +#endif + +void WellSwitchingLogger::gatherDataAndLog() +{ + +#if HAVE_MPI + if(cc_.size() == 1) + { + return; + } + + std::vector message_sizes; + std::vector well_name_lengths; + int message_size = calculateMessageSize(well_name_lengths); + + if ( cc_.rank() == 0 ){ + for(const auto& entry : switchMap_) + { + logSwitch(entry.first.c_str(), entry.second,0); + } + + message_sizes.resize(cc_.size()); + } + + MPI_Gather(&message_size, 1, MPI_INT, message_sizes.data(), + 1, MPI_INT, 0, MPI_COMM_WORLD); + + std::vector buffer(message_size); + packData(well_name_lengths, buffer); + + std::vector displ; + + if ( cc_.rank() == 0){ + // last entry will be total size of + displ.resize(cc_.size() + 1, 0); + std::partial_sum(message_sizes.begin(), message_sizes.end(), + displ.begin()+1); + } + std::vector recv_buffer; + if ( cc_.rank() == 0 ){ + recv_buffer.resize(displ[cc_.size()]); + } + + MPI_Gatherv(buffer.data(), buffer.size(), MPI_PACKED, + recv_buffer.data(), message_sizes.data(), + displ.data(), MPI_PACKED, 0, MPI_COMM_WORLD); + if ( cc_.rank() == 0 ) + { + unpackDataAndLog(recv_buffer, displ); + } +#endif +} + + +WellSwitchingLogger::~WellSwitchingLogger() +{ + gatherDataAndLog(); +} +} // end namespace wellhelpers +} // end namespace Opm diff --git a/opm/simulators/WellSwitchingLogger.hpp b/opm/simulators/WellSwitchingLogger.hpp new file mode 100644 index 000000000..46fb40974 --- /dev/null +++ b/opm/simulators/WellSwitchingLogger.hpp @@ -0,0 +1,114 @@ +/* + Copyright 2016 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 Statoil AS + + 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_WELLSWITCHINGLOGGER_HEADER_INCLUDED +#define OPM_WELLSWITCHINGLOGGER_HEADER_INCLUDED + +#include +#include +#include +#include + +#include + +#include + +#include + +namespace Opm +{ +namespace wellhelpers +{ + +/// \brief Utility class to handle the log messages about well switching. +/// +/// In parallel all the messages will be send to a root processor +/// and logged there. +class WellSwitchingLogger +{ + typedef std::map > SwitchMap; + +public: + /// \brief The type of the collective communication used. + typedef Dune::CollectiveCommunication + Communication; + + /// \brief Constructor. + /// + /// \param cc The collective communication to use. + explicit WellSwitchingLogger(const Communication& cc = + Dune::MPIHelper::getCollectiveCommunication()) + : cc_(cc) + {} + + /// \brief Log that a well switched. + /// \param name The name of the well. + /// \param from The control of the well before the switch. + /// \param to The control of the well after the switch. + void wellSwitched(std::string name, + WellControlType from, + WellControlType to) + { + if( cc_.size() > 1 ) + { + using Pair = typename SwitchMap::value_type; + switchMap_.insert(Pair(name, {char(from), char(to)})); + } + else + { + std::ostringstream ss; + ss << "Switching control mode for well " << name + << " from " << modestring[from] + << " to " << modestring[to] << std::endl; + OpmLog::info(ss.str()); + } + } + + /// \brief Destructor send does the actual logging. + ~WellSwitchingLogger(); + +private: + +#if HAVE_MPI + void unpackDataAndLog(std::vector& recv_buffer, + const std::vector& displ); + + void packData(std::vector& well_name_length, + std::vector& buffer); + + int calculateMessageSize(std::vector& well_name_length); + + void logSwitch(const char* name, std::array fromto, + int rank); + +#endif // HAVE_MPI + + void gatherDataAndLog(); + + /// \brief A map containing the local switches + SwitchMap switchMap_; + /// \brief Collective communication object. + Communication cc_; + /// \brief The strings for printing. + const std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; +}; +} // end namespace wellhelpers +} // end namespace Opm +#endif diff --git a/opm/simulators/thresholdPressures.hpp b/opm/simulators/thresholdPressures.hpp index b71f88ae4..6c7b5ff09 100644 --- a/opm/simulators/thresholdPressures.hpp +++ b/opm/simulators/thresholdPressures.hpp @@ -187,7 +187,7 @@ void computeMaxDp(std::map, double>& maxDp, int pvtRegionIdx = pvtRegion[cellIdx]; double T = initialState.temperature()[cellIdx]; - double p = initialState.pressure()[cellIdx]; + double p = phasePressure[wpos][cellIdx]; double b = pvtw.inverseFormationVolumeFactor(pvtRegionIdx, T, p); rho[wpos][cellIdx] = surfaceDensity[pvtRegionIdx][wpos]*b; @@ -201,21 +201,22 @@ void computeMaxDp(std::map, double>& maxDp, int pvtRegionIdx = pvtRegion[cellIdx]; double T = initialState.temperature()[cellIdx]; - double p = initialState.pressure()[cellIdx]; + double p = phasePressure[opos][cellIdx]; double Rs = initialState.gasoilratio()[cellIdx]; double RsSat = pvto.saturatedGasDissolutionFactor(pvtRegionIdx, T, p); + double b; if (Rs >= RsSat) { - double b = pvto.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p); - rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b; + b = pvto.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p); } else { - double b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); - rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b; - if (pu.phase_used[BlackoilPhases::Vapour]) { - int gpos = pu.phase_pos[BlackoilPhases::Vapour]; - rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b; - } + b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); + } + + rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b; + if (pu.phase_used[BlackoilPhases::Vapour]) { + int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b; } } } @@ -227,22 +228,21 @@ void computeMaxDp(std::map, double>& maxDp, int pvtRegionIdx = pvtRegion[cellIdx]; double T = initialState.temperature()[cellIdx]; - double p = initialState.pressure()[cellIdx]; + double p = phasePressure[gpos][cellIdx]; double Rv = initialState.rv()[cellIdx]; double RvSat = pvtg.saturatedOilVaporizationFactor(pvtRegionIdx, T, p); + double b; if (Rv >= RvSat) { - double b = pvtg.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p); - rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b; + b = pvtg.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p); } else { - double b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv); - rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b; - - if (pu.phase_used[BlackoilPhases::Liquid]) { - int opos = pu.phase_pos[BlackoilPhases::Liquid]; - rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b; - } + b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv); + } + rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b; + if (pu.phase_used[BlackoilPhases::Liquid]) { + int opos = pu.phase_pos[BlackoilPhases::Liquid]; + rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b; } } } @@ -270,7 +270,7 @@ void computeMaxDp(std::map, double>& maxDp, // update the maximum pressure potential difference between the two // regions - const auto barrierId = std::make_pair(eq1, eq2); + const auto barrierId = std::make_pair(std::min(eq1, eq2), std::max(eq1, eq2)); if (maxDp.count(barrierId) == 0) { maxDp[barrierId] = 0.0; } @@ -278,7 +278,6 @@ void computeMaxDp(std::map, double>& maxDp, for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { const double z1 = UgGridHelpers::cellCenterDepth(grid, c1); const double z2 = UgGridHelpers::cellCenterDepth(grid, c2); - const double zAvg = (z1 + z2)/2; // average depth const double rhoAvg = (rho[phaseIdx][c1] + rho[phaseIdx][c2])/2; @@ -289,8 +288,8 @@ void computeMaxDp(std::map, double>& maxDp, const double sResid2 = minSat[numPhases*c2 + phaseIdx]; // compute gravity corrected pressure potentials at the average depth - const double p1 = phasePressure[phaseIdx][c1] + rhoAvg*gravity*(zAvg - z1); - const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(zAvg - z2); + const double p1 = phasePressure[phaseIdx][c1]; + const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(z1 - z2); if ((p1 > p2 && s1 > sResid1) || (p2 > p1 && s2 > sResid2)) maxDp[barrierId] = std::max(maxDp[barrierId], std::abs(p1 - p2)); @@ -356,8 +355,11 @@ void computeMaxDp(std::map, double>& maxDp, // set the threshold pressure for faces of PVT regions where the third item // has been defaulted to the maximum pressure potential difference between // these regions - const auto barrierId = std::make_pair(eq1, eq2); - thpres_vals[face] = maxDp.at(barrierId); + const auto barrierId = std::make_pair(std::min(eq1, eq2), std::max(eq1, eq2)); + if (maxDp.count(barrierId) > 0) + thpres_vals[face] = maxDp.at(barrierId); + else + thpres_vals[face] = 0.0; } } diff --git a/redhat/opm-simulators.spec b/redhat/opm-simulators.spec index 80f527ab1..dd6bc1200 100644 --- a/redhat/opm-simulators.spec +++ b/redhat/opm-simulators.spec @@ -16,7 +16,7 @@ BuildRequires: blas-devel lapack-devel dune-common-devel opm-output-devel BuildRequires: git suitesparse-devel doxygen bc BuildRequires: opm-parser-devel opm-core-devel opm-grid-devel BuildRequires: tinyxml-devel dune-istl-devel eigen3-devel ert.ecl-devel -%{?el6:BuildRequires: cmake28 devtoolset-2 boost148-devel} +%{?el6:BuildRequires: cmake28 devtoolset-3-toolchain boost148-devel} %{!?el6:BuildRequires: cmake gcc gcc-c++ boost-devel} BuildRoot: %{_tmppath}/%{name}-%{version}-build Requires: libopm-simulators1 = %{version} @@ -59,8 +59,8 @@ This package contains the applications for opm-simulators %setup -q -n %{name}-release-%{version}-%{tag} %build -%{?el6:scl enable devtoolset-2 bash} -%{?el6:cmake28} %{?!el6:cmake} -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix} -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF %{?el6:-DCMAKE_CXX_COMPILER=/opt/rh/devtoolset-2/root/usr/bin/g++ -DCMAKE_C_COMPILER=/opt/rh/devtoolset-2/root/usr/bin/gcc -DBOOST_LIBRARYDIR=%{_libdir}/boost148 -DBOOST_INCLUDEDIR=%{_includedir}/boost148} +%{?el6:scl enable devtoolset-3 bash} +%{?el6:cmake28} %{?!el6:cmake} -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix} -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF %{?el6:-DCMAKE_CXX_COMPILER=/opt/rh/devtoolset-3/root/usr/bin/g++ -DCMAKE_C_COMPILER=/opt/rh/devtoolset-3/root/usr/bin/gcc -DBOOST_LIBRARYDIR=%{_libdir}/boost148 -DBOOST_INCLUDEDIR=%{_includedir}/boost148} make %install diff --git a/tests/test_multisegmentwells.cpp b/tests/test_multisegmentwells.cpp index 3f7d74dd3..e353c55a6 100644 --- a/tests/test_multisegmentwells.cpp +++ b/tests/test_multisegmentwells.cpp @@ -26,6 +26,7 @@ #endif #define BOOST_TEST_MODULE MultisegmentWellsTest +#define BOOST_TEST_NO_MAIN #include #include @@ -167,3 +168,16 @@ BOOST_AUTO_TEST_CASE(testStructure) BOOST_CHECK_EQUAL(0, ms_wells->topWellSegments()[0]); BOOST_CHECK_EQUAL(1, ms_wells->topWellSegments()[1]); } + +bool +init_unit_test_func() +{ + return true; +} + +int main(int argc, char** argv) +{ + Dune::MPIHelper::instance(argc, argv); + boost::unit_test::unit_test_main(&init_unit_test_func, + argc, argv); +} diff --git a/tests/test_wellswitchlogger.cpp b/tests/test_wellswitchlogger.cpp new file mode 100644 index 000000000..6e1252111 --- /dev/null +++ b/tests/test_wellswitchlogger.cpp @@ -0,0 +1,65 @@ +#include +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif +#define BOOST_TEST_MODULE DistributedCpGridTests +#define BOOST_TEST_NO_MAIN + +#include + +#include + +#if HAVE_MPI +class MPIError { +public: + /** @brief Constructor. */ + MPIError(std::string s, int e) : errorstring(s), errorcode(e){} + /** @brief The error string. */ + std::string errorstring; + /** @brief The mpi error code. */ + int errorcode; +}; + +void MPI_err_handler(MPI_Comm *, int *err_code, ...){ + char *err_string=new char[MPI_MAX_ERROR_STRING]; + int err_length; + MPI_Error_string(*err_code, err_string, &err_length); + std::string s(err_string, err_length); + std::cerr << "An MPI Error ocurred:"<