diff --git a/opm/autodiff/RateConverter.hpp b/opm/autodiff/RateConverter.hpp index 1ec218109..1265c8c2b 100644 --- a/opm/autodiff/RateConverter.hpp +++ b/opm/autodiff/RateConverter.hpp @@ -26,6 +26,7 @@ #include #include #include +#include #include @@ -37,7 +38,6 @@ #include #include #include - /** * \file * Facility for converting component rates at surface conditions to @@ -70,6 +70,52 @@ namespace Opm { }; } // Select + /** + * \brief Computes the temperature, pressure, and counter increment. + * + * In a parallel run only cells owned contribute to the cell average. + * \tparam is_parallel Whether this is a parallel run. + */ + template + struct AverageIncrementCalculator + { + /** + * \brief Computes the temperature, pressure, and counter increment. + * \param pressure The pressure. + * \param temperature The temperature. + * \param cell The current cell index. + * \param ownership A vector indicating whether a cell is owned + * by this process (value 1), or not (value 0). + * \param cell The cell index. + */ + std::tuple + operator()(const std::vector& pressure, + const std::vector& temperature, + const std::vector& ownership, + std::size_t cell){ + if ( ownership[cell] ) + { + return std::make_tuple(pressure[cell], + temperature[cell], 1); + } + else + { + return std::make_tuple(0, 0, 0); + } + } + }; + template<> + struct AverageIncrementCalculator + { + std::tuple + operator()(const std::vector& pressure, + const std::vector& temperature, + const std::vector&, + std::size_t cell){ + return std::make_tuple(pressure[cell], + temperature[cell], 1); + } + }; /** * Provide mapping from Region IDs to user-specified collection * of per-region attributes. @@ -378,14 +424,31 @@ namespace Opm { * reservoir voidage rate. * * \param[in] state Dynamic reservoir state. + * \param[in] any The information and communication utilities + * about/of the parallelization. in any parallel + * it wraps a ParallelISTLInformation. Parameter + * is optional. */ void - defineState(const BlackoilState& state) + defineState(const BlackoilState& state, + const boost::any& info = boost::any()) { - calcAverages(state); +#if HAVE_MPI + if( info.type() == typeid(ParallelISTLInformation) ) + { + const auto& ownership = + boost::any_cast(info) + .updateOwnerMask(state.pressure()); + calcAverages(state, info, ownership); + } + else +#endif + { + std::vector dummyOwnership; // not actually used + calcAverages(state, info, dummyOwnership); + } calcRmax(); } - /** * Region identifier. * @@ -551,10 +614,20 @@ namespace Opm { * Compute average hydrocarbon pressure and temperatures in all * regions. * - * \param[in] state Dynamic reservoir state. + * \param[in] state Dynamic reservoir state. + * \param[in] info The information and communication utilities + * about/of the parallelization. + * \param[in] ownership In a parallel run this is vector containing + * 1 for every owned unknown, zero otherwise. + * Not used in a sequential run. + * \tparam is_parallel True if the run is parallel. In this case + * info has to contain a ParallelISTLInformation + * object. */ + template void - calcAverages(const BlackoilState& state) + calcAverages(const BlackoilState& state, const boost::any& info, + const std::vector& ownerShip) { const auto& press = state.pressure(); const auto& temp = state.temperature(); @@ -567,17 +640,30 @@ namespace Opm { std::size_t n = 0; p = T = 0.0; for (const auto& cell : rmap_.cells(reg)) { - p += press[ cell ]; - T += temp [ cell ]; - - n += 1; + auto increment = Details:: + AverageIncrementCalculator()(press, temp, + ownerShip, + cell); + p += std::get<0>(increment); + T += std::get<1>(increment); + n += std::get<2>(increment); } - - p /= n; - T /= n; + std::size_t global_n = n; + double global_p = p; + double global_T = T; +#if HAVE_MPI + if ( is_parallel ) + { + const auto& real_info = boost::any_cast(info); + global_n = real_info.communicator().sum(n); + global_p = real_info.communicator().sum(p); + global_T = real_info.communicator().sum(T); + } +#endif + p = global_p / global_n; + T = global_T / global_n; } } - /** * Compute maximum dissolution and evaporation ratios at * average hydrocarbon pressure. diff --git a/opm/autodiff/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index 3154f73ed..f43e73607 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -365,12 +365,36 @@ namespace Opm const std::vector& resv_wells = SimFIBODetails::resvWells(wells, step, wmap); + const std::size_t number_resv_wells = resv_wells.size(); + std::size_t global_number_resv_wells = number_resv_wells; +#if HAVE_MPI + if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) + { + const auto& info = + boost::any_cast(solver_.parallelInformation()); + global_number_resv_wells = info.communicator().sum(global_number_resv_wells); + if ( global_number_resv_wells ) + { + // At least one process has resv wells. Therefore rate converter needs + // to calculate averages over regions that might cross process + // borders. This needs to be done by all processes and therefore + // outside of the next if statement. + rateConverter_.defineState(x, boost::any_cast(solver_.parallelInformation())); + } + } + else +#endif + { + if ( global_number_resv_wells ) + { + rateConverter_.defineState(x); + } + } + if (! resv_wells.empty()) { const PhaseUsage& pu = props_.phaseUsage(); const std::vector::size_type np = props_.numPhases(); - rateConverter_.defineState(x); - std::vector distr (np); std::vector hrates(np); std::vector prates(np);