Merge pull request #531 from blattms/parallel-RateConverter

Calculate parallel averages in RateConverter.

Closes issue #530
This commit is contained in:
Bård Skaflestad 2015-11-13 18:26:36 +01:00
commit 65da8b7e11
2 changed files with 126 additions and 16 deletions

View File

@ -26,6 +26,7 @@
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <Eigen/Core>
@ -37,7 +38,6 @@
#include <unordered_map>
#include <utility>
#include <vector>
/**
* \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<bool is_parallel>
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<double, double, int>
operator()(const std::vector<double>& pressure,
const std::vector<double>& temperature,
const std::vector<double>& 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<false>
{
std::tuple<double, double, int>
operator()(const std::vector<double>& pressure,
const std::vector<double>& temperature,
const std::vector<double>&,
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<const ParallelISTLInformation&>(info)
.updateOwnerMask(state.pressure());
calcAverages<true>(state, info, ownership);
}
else
#endif
{
std::vector<double> dummyOwnership; // not actually used
calcAverages<false>(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<bool is_parallel>
void
calcAverages(const BlackoilState& state)
calcAverages(const BlackoilState& state, const boost::any& info,
const std::vector<double>& 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<is_parallel>()(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<const ParallelISTLInformation&>(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.

View File

@ -365,12 +365,36 @@ namespace Opm
const std::vector<int>& 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<const ParallelISTLInformation&>(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<const ParallelISTLInformation&>(solver_.parallelInformation()));
}
}
else
#endif
{
if ( global_number_resv_wells )
{
rateConverter_.defineState(x);
}
}
if (! resv_wells.empty()) {
const PhaseUsage& pu = props_.phaseUsage();
const std::vector<double>::size_type np = props_.numPhases();
rateConverter_.defineState(x);
std::vector<double> distr (np);
std::vector<double> hrates(np);
std::vector<double> prates(np);