OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Calculate parallel averages in RateConverter.

Browse Source 
Previously, local averages were calculated and used in the
well equations. With this commit we add versions of defineState and
calcAverages that take into account the parallel domain decomposition
and calculate correct averages.

Function  calcAverages has a boolean template parameter
indicating whether this is a parallel run. Additionally we introduce
AverageIncrementCalculator with the same boolean template parameter.
In a parallel run we check whether the cell is owned by the process and
only in this case return an increment bigger than zero. In a sequential run
(no MPI or just one process -> empty boost::any parameter) no overhead is
introduced.
This commit is contained in:
Markus Blatt 2015-11-02 14:42:05 +01:00
parent 2a0051655c
commit f8715e31e7
2 changed files with 124 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

114
opm/autodiff/RateConverter.hpp
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<true>
{
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 othercase.
* If it is no such run then it is not used.
* \tparam is_parallel True if the runi 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.

26
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -365,12 +365,34 @@ 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) )
{
global_number_resv_wells = boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation()).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);