Prepare for Revised Implementation of WBPn

This initial commit changes the API of the CollectDataToIORank
class' handling of WBPn values from collecting a set cell pressures
into communicating pre-computed WBPn values through the new
WellBlockAveragePressures container class.  This is in preparation
of moving the WBPn calculation to the simulator side for greater
parallelism.  For now we do not compute any of the actual WBPn
values.  That will be the subject of follow-up commits.

While here, also split a number of very long lines for readability.
This commit is contained in:
Bård Skaflestad 2023-06-06 16:51:47 +02:00
parent b1ffc68853
commit 9d75915e4b
8 changed files with 307 additions and 257 deletions

View File

@ -531,55 +531,48 @@ public:
class PackUnPackWBPData : public P2PCommunicatorType::DataHandleInterface
{
const std::map<std::size_t, double>& localWBPData_;
std::map<std::size_t, double>& globalWBPValues_;
const data::WellBlockAveragePressures& localWBPData_;
data::WellBlockAveragePressures& globalWBPValues_;
public:
PackUnPackWBPData(const std::map<std::size_t, double>& localWBPData,
std::map<std::size_t, double>& globalWBPValues,
bool isIORank)
: localWBPData_(localWBPData)
PackUnPackWBPData(const data::WellBlockAveragePressures& localWBPData,
data::WellBlockAveragePressures& globalWBPValues,
const bool isIORank)
: localWBPData_ (localWBPData)
, globalWBPValues_(globalWBPValues)
{
if (isIORank) {
MessageBufferType buffer;
pack(0, buffer);
// pass a dummyLink to satisfy virtual class
int dummyLink = -1;
unpack(dummyLink, buffer);
if (! isIORank) {
return;
}
MessageBufferType buffer;
pack(0, buffer);
// Pass a dummy link to satisfy base class API requirement
const int dummyLink = -1;
unpack(dummyLink, buffer);
}
// pack all data associated with link
// Pack all data associated with link
void pack(int link, MessageBufferType& buffer)
{
// we should only get one link
if (link != 0)
throw std::logic_error("link in method pack is not 0 as expected");
// write all block data
unsigned int size = localWBPData_.size();
buffer.write(size);
for (const auto& [global_index, wbp_value] : localWBPData_) {
buffer.write(global_index);
buffer.write(wbp_value);
// We should only get one link
if (link != 0) {
throw std::logic_error {
"link (" + std::to_string(link) +
") in method pack() is not 0 as expected"
};
}
// Write all local, per-well, WBP data
this->localWBPData_.write(buffer);
}
// unpack all data associated with link
void unpack(int /*link*/, MessageBufferType& buffer)
// Unpack all data associated with link
void unpack([[maybe_unused]] const int link,
MessageBufferType& buffer)
{
// read all block data
unsigned int size = 0;
buffer.read(size);
for (size_t i = 0; i < size; ++i) {
std::size_t idx;
double data;
buffer.read(idx);
buffer.read(data);
globalWBPValues_[idx] = data;
}
this->globalWBPValues_.read(buffer);
}
};
@ -984,21 +977,21 @@ CollectDataToIORank(const Grid& grid, const EquilGrid* equilGrid,
template <class Grid, class EquilGrid, class GridView>
void CollectDataToIORank<Grid,EquilGrid,GridView>::
collect(const data::Solution& localCellData,
collect(const data::Solution& localCellData,
const std::map<std::pair<std::string, int>, double>& localBlockData,
const std::map<std::size_t, double>& localWBPData,
const data::Wells& localWellData,
const data::GroupAndNetworkValues& localGroupAndNetworkData,
const data::Aquifers& localAquiferData,
const WellTestState& localWellTestState,
const EclInterRegFlowMap& localInterRegFlows,
const data::Wells& localWellData,
const data::WellBlockAveragePressures& localWBPData,
const data::GroupAndNetworkValues& localGroupAndNetworkData,
const data::Aquifers& localAquiferData,
const WellTestState& localWellTestState,
const EclInterRegFlowMap& localInterRegFlows,
const std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>& localFlowsn,
const std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>& localFloresn)
{
globalCellData_ = {};
globalBlockData_.clear();
globalWBPData_.clear();
globalWellData_.clear();
globalWBPData_.values.clear();
globalGroupAndNetworkData_.clear();
globalAquiferData_.clear();
globalWellTestState_.clear();
@ -1055,8 +1048,8 @@ collect(const data::Solution& localCellData,
PackUnPackWBPData packUnpackWBPData {
localWBPData,
this->globalWBPData_,
this->isIORank()
this->globalWBPData_,
this->isIORank()
};
PackUnPackAquiferData packUnpackAquiferData {
@ -1110,28 +1103,37 @@ template <class Grid, class EquilGrid, class GridView>
int CollectDataToIORank<Grid,EquilGrid,GridView>::
localIdxToGlobalIdx(unsigned localIdx) const
{
if (!isParallel())
if (!isParallel()) {
return localIdx;
}
if (localIdxToGlobalIdx_.empty())
if (this->localIdxToGlobalIdx_.empty()) {
throw std::logic_error("index map is not created on this rank");
}
if (localIdx > localIdxToGlobalIdx_.size())
if (localIdx > this->localIdxToGlobalIdx_.size()) {
throw std::logic_error("local index is outside map range");
}
return localIdxToGlobalIdx_[localIdx];
return this->localIdxToGlobalIdx_[localIdx];
}
template <class Grid, class EquilGrid, class GridView>
bool CollectDataToIORank<Grid,EquilGrid,GridView>::
isCartIdxOnThisRank(int cartIdx) const
{
if (!isParallel())
if (! this->isParallel()) {
return true;
}
assert(!needsReordering);
auto candidate = std::lower_bound(sortedCartesianIdx_.begin(), sortedCartesianIdx_.end(), cartIdx);
return (candidate != sortedCartesianIdx_.end() && *candidate == cartIdx);
assert (! needsReordering);
auto candidate = std::lower_bound(this->sortedCartesianIdx_.begin(),
this->sortedCartesianIdx_.end(),
cartIdx);
return (candidate != sortedCartesianIdx_.end())
&& (*candidate == cartIdx);
}

View File

@ -25,15 +25,21 @@
#include <opm/output/data/Aquifer.hpp>
#include <opm/output/data/Cells.hpp>
#include <opm/output/data/Groups.hpp>
#include <opm/output/data/Solution.hpp>
#include <opm/output/data/Wells.hpp>
#include <opm/output/data/Groups.hpp>
#include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
#include <opm/grid/common/p2pcommunicator.hh>
#include <ebos/eclinterregflows.hh>
#include <array>
#include <map>
#include <set>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
@ -67,20 +73,17 @@ public:
const std::set<std::string>& fipRegionsInterregFlow = {});
// gather solution to rank 0 for EclipseWriter
void collect(const data::Solution& localCellData,
void collect(const data::Solution& localCellData,
const std::map<std::pair<std::string, int>, double>& localBlockData,
const std::map<std::size_t, double>& localWBPData,
const data::Wells& localWellData,
const data::GroupAndNetworkValues& localGroupAndNetworkData,
const data::Aquifers& localAquiferData,
const WellTestState& localWellTestState,
const EclInterRegFlowMap& interRegFlows,
const data::Wells& localWellData,
const data::WellBlockAveragePressures& localWBPData,
const data::GroupAndNetworkValues& localGroupAndNetworkData,
const data::Aquifers& localAquiferData,
const WellTestState& localWellTestState,
const EclInterRegFlowMap& interRegFlows,
const std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>& localFlowsn,
const std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>& localFloresn);
const std::map<std::size_t, double>& globalWBPData() const
{ return this->globalWBPData_; }
const std::map<std::pair<std::string, int>, double>& globalBlockData() const
{ return globalBlockData_; }
@ -90,8 +93,8 @@ public:
const data::Wells& globalWellData() const
{ return globalWellData_; }
const WellTestState& globalWellTestState() const
{ return this->globalWellTestState_; }
const data::WellBlockAveragePressures& globalWBPData() const
{ return this->globalWBPData_; }
const data::GroupAndNetworkValues& globalGroupAndNetworkData() const
{ return globalGroupAndNetworkData_; }
@ -99,6 +102,9 @@ public:
const data::Aquifers& globalAquiferData() const
{ return globalAquiferData_; }
const WellTestState& globalWellTestState() const
{ return this->globalWellTestState_; }
EclInterRegFlowMap& globalInterRegFlows()
{ return this->globalInterRegFlows_; }
@ -144,8 +150,8 @@ protected:
std::vector<int> globalRanks_;
data::Solution globalCellData_;
std::map<std::pair<std::string, int>, double> globalBlockData_;
std::map<std::size_t, double> globalWBPData_;
data::Wells globalWellData_;
data::WellBlockAveragePressures globalWBPData_;
data::GroupAndNetworkValues globalGroupAndNetworkData_;
data::Aquifers globalAquiferData_;
WellTestState globalWellTestState_;

View File

@ -234,11 +234,6 @@ public:
return this->interRegionFlows_.wantInterRegflowSummary();
}
const std::map<std::size_t, double>& getWBPData() const
{
return this->wbpData_;
}
const std::map<std::pair<std::string, int>, double>& getBlockData()
{
return blockData_;
@ -552,7 +547,6 @@ protected:
std::map<size_t, Scalar> waterConnectionSaturations_;
std::map<size_t, Scalar> gasConnectionSaturations_;
std::map<std::pair<std::string, int>, double> blockData_;
std::map<std::size_t , double> wbpData_;
std::optional<Inplace> initialInplace_;
bool local_data_valid_;

View File

@ -22,12 +22,14 @@
*/
#include <config.h>
#include <ebos/eclgenericwriter.hh>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/cpgrid/GridHelpers.hpp>
#include <opm/grid/polyhedralgrid.hh>
#include <opm/grid/utility/cartesianToCompressed.hpp>
#if HAVE_DUNE_ALUGRID
#include "eclalugridvanguard.hh"
#include <dune/alugrid/grid.hh>
@ -40,13 +42,14 @@
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#include <opm/input/eclipse/Schedule/Action/State.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Schedule/SummaryState.hpp>
#include <opm/input/eclipse/Schedule/UDQ/UDQConfig.hpp>
#include <opm/input/eclipse/Schedule/UDQ/UDQState.hpp>
#include <opm/input/eclipse/Schedule/Well/PAvgCalculatorCollection.hpp>
#include <opm/input/eclipse/Schedule/Well/WellMatcher.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <dune/grid/common/mcmgmapper.hh>
@ -65,7 +68,12 @@
#include <mpi.h>
#endif
#include <algorithm>
#include <array>
#include <cassert>
#include <functional>
#include <map>
#include <memory>
#include <string>
#include <unordered_map>
#include <utility>
@ -169,22 +177,22 @@ struct EclWriteTasklet : public Opm::TaskletInterface
, reportStepNum_(reportStepNum)
, isSubStep_(isSubStep)
, secondsElapsed_(secondsElapsed)
, restartValue_(restartValue)
, restartValue_(std::move(restartValue))
, writeDoublePrecision_(writeDoublePrecision)
{ }
{}
// callback to eclIO serial writeTimeStep method
void run()
{
eclIO_.writeTimeStep(actionState_,
wtestState_,
summaryState_,
udqState_,
reportStepNum_,
isSubStep_,
secondsElapsed_,
restartValue_,
writeDoublePrecision_);
this->eclIO_.writeTimeStep(this->actionState_,
this->wtestState_,
this->summaryState_,
this->udqState_,
this->reportStepNum_,
this->isSubStep_,
this->secondsElapsed_,
std::move(this->restartValue_),
this->writeDoublePrecision_);
}
};
@ -210,38 +218,30 @@ EclGenericWriter(const Schedule& schedule,
cartMapper,
equilCartMapper,
summaryConfig.fip_regions_interreg_flow())
, grid_(grid)
, gridView_(gridView)
, schedule_(schedule)
, eclState_(eclState)
, summaryConfig_(summaryConfig)
, cartMapper_(cartMapper)
, grid_ (grid)
, gridView_ (gridView)
, schedule_ (schedule)
, eclState_ (eclState)
, summaryConfig_ (summaryConfig)
, cartMapper_ (cartMapper)
, equilCartMapper_(equilCartMapper)
, equilGrid_(equilGrid)
, equilGrid_ (equilGrid)
{
if (collectToIORank_.isIORank()) {
eclIO_.reset(new EclipseIO(eclState_,
UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()),
schedule_,
summaryConfig_, "", enableEsmry));
if (this->collectToIORank_.isIORank()) {
this->eclIO_ = std::make_unique<EclipseIO>
(this->eclState_,
UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()),
this->schedule_, this->summaryConfig_, "", enableEsmry);
}
const auto& wbp_calculators = eclIO_->summary().wbp_calculators( schedule.size() - 1 );
wbp_index_list_ = wbp_calculators.index_list();
}
if (collectToIORank_.isParallel()) {
const auto& comm = grid_.comm();
unsigned long size = wbp_index_list_.size();
comm.broadcast(&size, 1, collectToIORank_.ioRank);
if (!collectToIORank_.isIORank())
wbp_index_list_.resize( size );
comm.broadcast(wbp_index_list_.data(), size, collectToIORank_.ioRank);
}
// create output thread if enabled and rank is I/O rank
// async output is enabled by default if pthread are enabled
int numWorkerThreads = 0;
if (enableAsyncOutput && collectToIORank_.isIORank())
if (enableAsyncOutput && collectToIORank_.isIORank()) {
numWorkerThreads = 1;
taskletRunner_.reset(new TaskletRunner(numWorkerThreads));
}
this->taskletRunner_.reset(new TaskletRunner(numWorkerThreads));
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
@ -257,12 +257,18 @@ void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
writeInit(const std::function<unsigned int(unsigned int)>& map)
{
if (collectToIORank_.isIORank()) {
std::map<std::string, std::vector<int> > integerVectors;
if (collectToIORank_.isParallel())
std::map<std::string, std::vector<int>> integerVectors;
if (collectToIORank_.isParallel()) {
integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks());
}
auto cartMap = cartesianToCompressed(equilGrid_->size(0), UgGridHelpers::globalCell(*equilGrid_));
eclIO_->writeInitial(computeTrans_(cartMap, map), integerVectors, exportNncStructure_(cartMap, map));
eclIO_->writeInitial(computeTrans_(cartMap, map),
integerVectors,
exportNncStructure_(cartMap, map));
}
#if HAVE_MPI
if (collectToIORank_.isParallel()) {
const auto& comm = grid_.comm();
@ -450,8 +456,9 @@ doWriteOutput(const int reportStepNum,
const bool needsReordering = this->collectToIORank_.doesNeedReordering();
RestartValue restartValue {
(isParallel || needsReordering) ? this->collectToIORank_.globalCellData()
: std::move(localCellData),
(isParallel || needsReordering)
? this->collectToIORank_.globalCellData()
: std::move(localCellData),
isParallel ? this->collectToIORank_.globalWellData()
: std::move(localWellData),
@ -515,32 +522,32 @@ doWriteOutput(const int reportStepNum,
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
evalSummary(const int reportStepNum,
const Scalar curTime,
const std::map<std::size_t, double>& wbpData,
const data::Wells& localWellData,
const data::GroupAndNetworkValues& localGroupAndNetworkData,
const std::map<int,data::AquiferData>& localAquiferData,
evalSummary(const int reportStepNum,
const Scalar curTime,
const data::Wells& localWellData,
const data::WellBlockAveragePressures& localWBPData,
const data::GroupAndNetworkValues& localGroupAndNetworkData,
const std::map<int,data::AquiferData>& localAquiferData,
const std::map<std::pair<std::string, int>, double>& blockData,
const std::map<std::string, double>& miscSummaryData,
const std::map<std::string, std::vector<double>>& regionData,
const Inplace& inplace,
const Inplace& initialInPlace,
const EclInterRegFlowMap& interRegionFlowMap,
SummaryState& summaryState,
UDQState& udqState)
const std::map<std::string, double>& miscSummaryData,
const std::map<std::string, std::vector<double>>& regionData,
const Inplace& inplace,
const Inplace& initialInPlace,
const EclInterRegFlowMap& interRegFlows,
SummaryState& summaryState,
UDQState& udqState)
{
if (collectToIORank_.isIORank()) {
const auto& summary = eclIO_->summary();
auto wbp_calculators = summary.wbp_calculators(reportStepNum);
for (const auto& [global_index, pressure] : wbpData)
wbp_calculators.add_pressure( global_index, pressure );
const auto& wellData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalWellData()
: localWellData;
const auto& wbpData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalWBPData()
: localWBPData;
const auto& groupAndNetworkData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalGroupAndNetworkData()
: localGroupAndNetworkData;
@ -553,24 +560,29 @@ evalSummary(const int reportStepNum,
reportStepNum,
curTime,
wellData,
wbpData,
groupAndNetworkData,
miscSummaryData,
initialInPlace,
inplace,
wbp_calculators,
regionData,
blockData,
aquiferData,
getInterRegFlowsAsMap(interRegionFlowMap));
getInterRegFlowsAsMap(interRegFlows));
// Off-by-one-fun: The reportStepNum argument corresponds to the
// report step these results will be written to, whereas the
// argument to UDQ function evaluation corresponds to the report
// step we are currently on.
auto udq_step = reportStepNum - 1;
const auto& udq_config = schedule_.getUDQConfig(udq_step);
udq_config.eval( udq_step, schedule_.wellMatcher(udq_step), summaryState, udqState);
const auto udq_step = reportStepNum - 1;
this->schedule_.getUDQConfig(udq_step)
.eval(udq_step,
this->schedule_.wellMatcher(udq_step),
summaryState,
udqState);
}
#if HAVE_MPI
if (collectToIORank_.isParallel()) {
EclMpiSerializer ser(grid_.comm());

View File

@ -116,10 +116,10 @@ protected:
data::GroupAndNetworkValues&& localGroupAndNetworkData,
data::Aquifers&& localAquiferData,
WellTestState&& localWTestState,
const Action::State& actionState,
const UDQState& udqState,
const SummaryState& summaryState,
const std::vector<Scalar>& thresholdPressure,
const Action::State& actionState,
const UDQState& udqState,
const SummaryState& summaryState,
const std::vector<Scalar>& thresholdPressure,
Scalar curTime,
Scalar nextStepSize,
bool doublePrecision,
@ -128,20 +128,20 @@ protected:
bool isFloresn,
std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>&& floresn);
void evalSummary(int reportStepNum,
Scalar curTime,
const std::map<std::size_t, double>& wbpData,
const data::Wells& localWellData,
const data::GroupAndNetworkValues& localGroupAndNetworkData,
const std::map<int,data::AquiferData>& localAquiferData,
void evalSummary(int reportStepNum,
Scalar curTime,
const data::Wells& localWellData,
const data::WellBlockAveragePressures& localWBPData,
const data::GroupAndNetworkValues& localGroupAndNetworkData,
const std::map<int,data::AquiferData>& localAquiferData,
const std::map<std::pair<std::string, int>, double>& blockData,
const std::map<std::string, double>& miscSummaryData,
const std::map<std::string, std::vector<double>>& regionData,
const Inplace& inplace,
const Inplace& initialInPlace,
const EclInterRegFlowMap& interRegionFlowMap,
SummaryState& summaryState,
UDQState& udqState);
const std::map<std::string, double>& miscSummaryData,
const std::map<std::string, std::vector<double>>& regionData,
const Inplace& inplace,
const Inplace& initialInPlace,
const EclInterRegFlowMap& interRegFlows,
SummaryState& summaryState,
UDQState& udqState);
CollectDataToIORankType collectToIORank_;
const Grid& grid_;
@ -156,7 +156,6 @@ protected:
const Dune::CartesianIndexMapper<Grid>& cartMapper_;
const Dune::CartesianIndexMapper<EquilGrid>* equilCartMapper_;
const EquilGrid* equilGrid_;
std::vector<std::size_t> wbp_index_list_;
SimulatorReportSingle sub_step_report_;
SimulatorReport simulation_report_;
mutable std::vector<NNCdata> outputNnc_;

View File

@ -138,7 +138,6 @@ class EclOutputBlackOilModule : public EclGenericOutputBlackoilModule<GetPropTyp
public:
template <class CollectDataToIORankType>
EclOutputBlackOilModule(const Simulator& simulator,
const std::vector<std::size_t>& wbp_index_list,
const CollectDataToIORankType& collectToIORank)
: BaseType(simulator.vanguard().eclState(),
simulator.vanguard().schedule(),
@ -159,22 +158,17 @@ public:
this->createLocalRegion_(region_pair.second);
}
auto isCartIdxOnThisRank = [&collectToIORank](int idx)
{
return collectToIORank.isCartIdxOnThisRank(idx);
};
auto isCartIdxOnThisRank = [&collectToIORank](const int idx) {
return collectToIORank.isCartIdxOnThisRank(idx);
};
this->setupBlockData(isCartIdxOnThisRank);
for (const auto& global_index : wbp_index_list) {
if (collectToIORank.isCartIdxOnThisRank(global_index - 1)) {
this->wbpData_.emplace(global_index, 0.0);
}
}
this->forceDisableFipOutput_ =
EWOMS_GET_PARAM(TypeTag, bool, ForceDisableFluidInPlaceOutput);
this->forceDisableFipOutput_ = EWOMS_GET_PARAM(TypeTag, bool, ForceDisableFluidInPlaceOutput);
this->forceDisableFipresvOutput_ = EWOMS_GET_PARAM(TypeTag, bool, ForceDisableResvFluidInPlaceOutput);
this->forceDisableFipresvOutput_ =
EWOMS_GET_PARAM(TypeTag, bool, ForceDisableResvFluidInPlaceOutput);
}
/*!
@ -199,25 +193,32 @@ public:
* write to ECL output files
*/
void
allocBuffers(unsigned bufferSize, unsigned reportStepNum, const bool substep, const bool log, const bool isRestart)
allocBuffers(const unsigned bufferSize,
const unsigned reportStepNum,
const bool substep,
const bool log,
const bool isRestart)
{
if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value)
if (! std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
return;
}
const auto& problem = this->simulator_.problem();
this->doAllocBuffers(bufferSize,
reportStepNum,
substep,
log,
isRestart,
simulator_.problem().vapparsActive(std::max(simulator_.episodeIndex(), 0)),
simulator_.problem().materialLawManager()->enableHysteresis(),
simulator_.problem().tracerModel().numTracers(),
simulator_.problem().eclWriter()->getOutputNnc().size());
problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)),
problem.materialLawManager()->enableHysteresis(),
problem.tracerModel().numTracers(),
problem.eclWriter()->getOutputNnc().size());
}
/*!
* \brief Modify the internal buffers according to the intensive quanties relevant
* for an element
* \brief Modify the internal buffers according to the intensive
* quanties relevant for an element
*/
void processElement(const ElementContext& elemCtx)
{
@ -550,7 +551,6 @@ public:
= FluidSystem::viscosity(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
}
// Adding Well RFT data
const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
if (this->oilConnectionPressures_.count(cartesianIdx) > 0) {
@ -562,25 +562,17 @@ public:
if (this->gasConnectionSaturations_.count(cartesianIdx) > 0) {
this->gasConnectionSaturations_[cartesianIdx] = getValue(fs.saturation(gasPhaseIdx));
}
if (this->wbpData_.count(cartesianIdx) > 0) {
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
this->wbpData_[cartesianIdx] = getValue(fs.pressure(oilPhaseIdx));
} else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
this->wbpData_[cartesianIdx] = getValue(fs.pressure(gasPhaseIdx));
} else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
this->wbpData_[cartesianIdx] = getValue(fs.pressure(waterPhaseIdx));
}
}
// tracers
const auto& tracerModel = simulator_.problem().tracerModel();
if (!this->tracerConcentrations_.empty()) {
for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); tracerIdx++) {
if (this->tracerConcentrations_[tracerIdx].empty())
if (! this->tracerConcentrations_.empty()) {
for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
if (this->tracerConcentrations_[tracerIdx].empty()) {
continue;
}
this->tracerConcentrations_[tracerIdx][globalDofIdx]
= tracerModel.tracerConcentration(tracerIdx, globalDofIdx);
this->tracerConcentrations_[tracerIdx][globalDofIdx] =
tracerModel.tracerConcentration(tracerIdx, globalDofIdx);
}
}
}

View File

@ -148,22 +148,28 @@ public:
simulator.vanguard().eclState(),
simulator.vanguard().summaryConfig(),
simulator.vanguard().grid(),
simulator.vanguard().grid().comm().rank() == 0 ? &simulator.vanguard().equilGrid() : nullptr,
((simulator.vanguard().grid().comm().rank() == 0)
? &simulator.vanguard().equilGrid()
: nullptr),
simulator.vanguard().gridView(),
simulator.vanguard().cartesianIndexMapper(),
simulator.vanguard().grid().comm().rank() == 0 ? &simulator.vanguard().equilCartesianIndexMapper() : nullptr,
EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput), EWOMS_GET_PARAM(TypeTag, bool, EnableEsmry))
((simulator.vanguard().grid().comm().rank() == 0)
? &simulator.vanguard().equilCartesianIndexMapper()
: nullptr),
EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput),
EWOMS_GET_PARAM(TypeTag, bool, EnableEsmry))
, simulator_(simulator)
{
#ifdef HAVE_DAMARIS
this->damarisUpdate_ = enableDamarisOutput_();
#endif
this->eclOutputModule_ = std::make_unique<EclOutputBlackOilModule<TypeTag>>(simulator, this->wbp_index_list_, this->collectToIORank_);
this->wbp_index_list_.clear();
this->eclOutputModule_ = std::make_unique<EclOutputBlackOilModule<TypeTag>>
(simulator, this->collectToIORank_);
}
~EclWriter()
{ }
{}
const EquilGrid& globalGrid() const
{
@ -204,6 +210,7 @@ public:
simulator_.vanguard().setupTime();
const auto localWellData = simulator_.problem().wellModel().wellData();
const auto localWBP = data::WellBlockAveragePressures{};
const auto localGroupAndNetworkData = simulator_.problem().wellModel()
.groupAndNetworkData(reportStepNum);
@ -220,8 +227,8 @@ public:
this->collectToIORank_.collect({},
eclOutputModule_->getBlockData(),
eclOutputModule_->getWBPData(),
localWellData,
localWBP,
localGroupAndNetworkData,
localAquiferData,
localWellTestState,
@ -283,26 +290,34 @@ public:
if (this->simulation_report_.success.total_newton_iterations != 0) {
miscSummaryData["MSUMNEWT"] = this->simulation_report_.success.total_newton_iterations;
}
{
OPM_TIMEBLOCK(evalSummary);
this->evalSummary(reportStepNum, curTime,
this->collectToIORank_.isParallel() ?
this->collectToIORank_.globalWBPData() :
this->eclOutputModule_->getWBPData(),
localWellData,
localGroupAndNetworkData,
localAquiferData,
this->collectToIORank_.isParallel() ?
this->collectToIORank_.globalBlockData() :
this->eclOutputModule_->getBlockData(),
miscSummaryData, regionData,
inplace,
eclOutputModule_->initialInplace(),
this->collectToIORank_.isParallel()
? this->collectToIORank_.globalInterRegFlows()
: this->eclOutputModule_->getInterRegFlows(),
summaryState(), udqState());
OPM_TIMEBLOCK(evalSummary);
const auto& blockData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalBlockData()
: this->eclOutputModule_->getBlockData();
const auto& interRegFlows = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalInterRegFlows()
: this->eclOutputModule_->getInterRegFlows();
this->evalSummary(reportStepNum,
curTime,
localWellData,
localWBP,
localGroupAndNetworkData,
localAquiferData,
blockData,
miscSummaryData,
regionData,
inplace,
this->eclOutputModule_->initialInplace(),
interRegFlows,
this->summaryState(),
this->udqState());
}
{
OPM_TIMEBLOCK(outputXXX);
eclOutputModule_->outputProdLog(reportStepNum, isSubStep, forceDisableProdOutput);
@ -314,9 +329,11 @@ public:
void writeOutput(bool isSubStep)
{
OPM_TIMEBLOCK(writeOutput);
const int reportStepNum = simulator_.episodeIndex() + 1;
this->prepareLocalCellData(isSubStep, reportStepNum);
this->eclOutputModule_->outputErrorLog(simulator_.gridView().comm());
#ifdef HAVE_DAMARIS
if (EWOMS_GET_PARAM(TypeTag, bool, EnableDamarisOutput)) {
// N.B. damarisUpdate_ should be set to true if at any time the model geometry changes
@ -340,35 +357,44 @@ public:
}
}
#endif
// output using eclWriter if enabled
auto localWellData = simulator_.problem().wellModel().wellData();
auto localGroupAndNetworkData = simulator_.problem().wellModel()
.groupAndNetworkData(reportStepNum);
auto localAquiferData = simulator_.problem().aquiferModel().aquiferData();
auto localWellTestState = simulator_.problem().wellModel().wellTestState();
auto flowsn = this->eclOutputModule_->getFlowsn();
const bool isFlowsn = this->eclOutputModule_->hasFlowsn();
auto floresn = this->eclOutputModule_->getFloresn();
auto flowsn = this->eclOutputModule_->getFlowsn();
const bool isFloresn = this->eclOutputModule_->hasFloresn();
auto floresn = this->eclOutputModule_->getFloresn();
data::Solution localCellData = {};
if (! isSubStep) {
this->eclOutputModule_->assignToSolution(localCellData);
// add cell data to perforations for Rft output
// Add cell data to perforations for RFT output
this->eclOutputModule_->addRftDataToWells(localWellData, reportStepNum);
}
if (this->collectToIORank_.isParallel()|| this->collectToIORank_.doesNeedReordering()) {
if (this->collectToIORank_.isParallel() ||
this->collectToIORank_.doesNeedReordering())
{
// Note: We don't need WBP (well-block averaged pressures) or
// inter-region flow rate values in order to create restart file
// output. There's consequently no need to collect those
// properties on the I/O rank.
this->collectToIORank_.collect(localCellData,
eclOutputModule_->getBlockData(),
eclOutputModule_->getWBPData(),
this->eclOutputModule_->getBlockData(),
localWellData,
/* wbpData = */ {},
localGroupAndNetworkData,
localAquiferData,
localWellTestState,
{},
/* interRegFlows = */ {},
flowsn,
floresn);
}
@ -376,6 +402,7 @@ public:
if (this->collectToIORank_.isIORank()) {
const Scalar curTime = simulator_.time() + simulator_.timeStepSize();
const Scalar nextStepSize = simulator_.problem().nextTimeStepSize();
this->doWriteOutput(reportStepNum, isSubStep,
std::move(localCellData),
std::move(localWellData),
@ -385,13 +412,11 @@ public:
this->actionState(),
this->udqState(),
this->summaryState(),
simulator_.problem().thresholdPressure().getRestartVector(),
this->simulator_.problem().thresholdPressure().getRestartVector(),
curTime, nextStepSize,
EWOMS_GET_PARAM(TypeTag, bool, EclOutputDoublePrecision),
isFlowsn,
std::move(flowsn),
isFloresn,
std::move(floresn));
isFlowsn, std::move(flowsn),
isFloresn, std::move(floresn));
}
}
@ -487,7 +512,7 @@ public:
Scalar restartTimeStepSize() const
{ return restartTimeStepSize_; }
template<class Serializer>
template <class Serializer>
void serializeOp(Serializer& serializer)
{
serializer(*eclOutputModule_);
@ -496,10 +521,12 @@ public:
private:
static bool enableEclOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
#ifdef HAVE_DAMARIS
static bool enableDamarisOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EnableDamarisOutput); }
#endif
const EclipseState& eclState() const
{ return simulator_.vanguard().eclState(); }
@ -519,7 +546,8 @@ private:
const int reportStepNum)
{
OPM_TIMEBLOCK(prepareLocalCellData);
if (eclOutputModule_->localDataValid()) {
if (this->eclOutputModule_->localDataValid()) {
return;
}
@ -527,54 +555,68 @@ private:
const int numElements = gridView.size(/*codim=*/0);
const bool log = this->collectToIORank_.isIORank();
eclOutputModule_->allocBuffers(numElements, reportStepNum,
isSubStep, log, /*isRestart*/ false);
this->eclOutputModule_->
allocBuffers(numElements, reportStepNum,
isSubStep, log, /*isRestart*/ false);
ElementContext elemCtx(simulator_);
OPM_BEGIN_PARALLEL_TRY_CATCH();
{
OPM_TIMEBLOCK(prepareCellBasedData);
for (const auto& elem : elements(gridView)) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
eclOutputModule_->processElement(elemCtx);
OPM_BEGIN_PARALLEL_TRY_CATCH();
{
OPM_TIMEBLOCK(prepareCellBasedData);
for (const auto& elem : elements(gridView)) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
this->eclOutputModule_->processElement(elemCtx);
}
}
}
if(!simulator_.model().linearizer().getFlowsInfo().empty()){
if (! this->simulator_.model().linearizer().getFlowsInfo().empty()) {
OPM_TIMEBLOCK(prepareFlowsData);
for (const auto& elem : elements(gridView)) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
eclOutputModule_->processElementFlows(elemCtx);
this->eclOutputModule_->processElementFlows(elemCtx);
}
}
{
OPM_TIMEBLOCK(prepareBlockData);
for (const auto& elem : elements(gridView)) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
eclOutputModule_->processElementBlockData(elemCtx);
}
OPM_TIMEBLOCK(prepareBlockData);
for (const auto& elem : elements(gridView)) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
this->eclOutputModule_->processElementBlockData(elemCtx);
}
}
{
OPM_TIMEBLOCK(prepareFluidInPlace);
OPM_TIMEBLOCK(prepareFluidInPlace);
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (int dofIdx=0; dofIdx < numElements; ++dofIdx){
const auto& intQuants = *(simulator_.model().cachedIntensiveQuantities(dofIdx, /*timeIdx=*/0));
for (int dofIdx = 0; dofIdx < numElements; ++dofIdx) {
const auto& intQuants = *simulator_.model().cachedIntensiveQuantities(dofIdx, /*timeIdx=*/0);
const auto totVolume = simulator_.model().dofTotalVolume(dofIdx);
eclOutputModule_->updateFluidInPlace(dofIdx, intQuants, totVolume);
this->eclOutputModule_->updateFluidInPlace(dofIdx, intQuants, totVolume);
}
}
}
eclOutputModule_->validateLocalData();
OPM_END_PARALLEL_TRY_CATCH("EclWriter::prepareLocalCellData() failed: ", simulator_.vanguard().grid().comm());
this->eclOutputModule_->validateLocalData();
OPM_END_PARALLEL_TRY_CATCH("EclWriter::prepareLocalCellData() failed: ",
this->simulator_.vanguard().grid().comm());
}
void captureLocalFluxData()
{
OPM_TIMEBLOCK(captureLocalData);
const auto& gridView = this->simulator_.vanguard().gridView();
const auto timeIdx = 0u;
@ -612,6 +654,7 @@ private:
Simulator& simulator_;
std::unique_ptr<EclOutputBlackOilModule<TypeTag>> eclOutputModule_;
Scalar restartTimeStepSize_;
#ifdef HAVE_DAMARIS
bool damarisUpdate_ = false; ///< Whenever this is true writeOutput() will set up Damaris offsets of model fields
#endif

View File

@ -211,6 +211,8 @@ TEST_FOR_TYPE_NAMED(data::SegmentPressures, SegmentPressures)
TEST_FOR_TYPE_NAMED(data::Solution, Solution)
TEST_FOR_TYPE_NAMED(data::Well, dataWell)
TEST_FOR_TYPE_NAMED(data::Wells, Wells)
TEST_FOR_TYPE_NAMED(data::WellBlockAvgPress, dataWBPObject)
TEST_FOR_TYPE_NAMED(data::WellBlockAveragePressures, dataWBPCollection)
TEST_FOR_TYPE(Deck)
TEST_FOR_TYPE(DeckItem)
TEST_FOR_TYPE(DeckKeyword)