23#ifndef OPM_ECL_GENERIC_WRITER_IMPL_HPP
24#define OPM_ECL_GENERIC_WRITER_IMPL_HPP
26#include <dune/grid/common/mcmgmapper.hh>
28#include <opm/grid/GridHelpers.hpp>
29#include <opm/grid/utility/cartesianToCompressed.hpp>
31#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
32#include <opm/input/eclipse/EclipseState/Grid/RegionSetMatcher.hpp>
33#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
35#include <opm/input/eclipse/Schedule/Action/State.hpp>
36#include <opm/input/eclipse/Schedule/Schedule.hpp>
37#include <opm/input/eclipse/Schedule/SummaryState.hpp>
38#include <opm/input/eclipse/Schedule/UDQ/UDQConfig.hpp>
39#include <opm/input/eclipse/Schedule/UDQ/UDQState.hpp>
40#include <opm/input/eclipse/Schedule/Well/WellMatcher.hpp>
42#include <opm/input/eclipse/Units/UnitSystem.hpp>
44#include <opm/output/eclipse/EclipseIO.hpp>
45#include <opm/output/eclipse/RestartValue.hpp>
46#include <opm/output/eclipse/Summary.hpp>
66#include <unordered_map>
86bool directVerticalNeighbors(
const std::array<int, 3>& cartDims,
87 const std::unordered_map<int,int>& cartesianToActive,
88 int smallGlobalIndex,
int largeGlobalIndex)
90 assert(smallGlobalIndex <= largeGlobalIndex);
91 std::array<int, 3> ijk1, ijk2;
92 auto globalToIjk = [cartDims](
int gc) {
93 std::array<int, 3> ijk;
94 ijk[0] = gc % cartDims[0];
96 ijk[1] = gc % cartDims[1];
97 ijk[2] = gc / cartDims[1];
100 ijk1 = globalToIjk(smallGlobalIndex);
101 ijk2 = globalToIjk(largeGlobalIndex);
102 assert(ijk2[2]>=ijk1[2]);
104 if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1)
106 assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0);
107 for (
int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex;
108 gi += cartDims[0] * cartDims[1] )
110 if ( cartesianToActive.find( gi ) != cartesianToActive.end() )
120std::unordered_map<std::string, Opm::data::InterRegFlowMap>
123 auto maps = std::unordered_map<std::string, Opm::data::InterRegFlowMap>{};
125 const auto& regionNames = map.
names();
127 const auto nmap = regionNames.size();
130 for (
auto mapID = 0*nmap; mapID < nmap; ++mapID) {
131 maps.emplace(regionNames[mapID], std::move(flows[mapID]));
137struct EclWriteTasklet :
public Opm::TaskletInterface
139 Opm::Action::State actionState_;
140 Opm::WellTestState wtestState_;
141 Opm::SummaryState summaryState_;
142 Opm::UDQState udqState_;
143 Opm::EclipseIO& eclIO_;
146 double secondsElapsed_;
147 Opm::RestartValue restartValue_;
148 bool writeDoublePrecision_;
150 explicit EclWriteTasklet(
const Opm::Action::State& actionState,
151 const Opm::WellTestState& wtestState,
152 const Opm::SummaryState& summaryState,
153 const Opm::UDQState& udqState,
154 Opm::EclipseIO& eclIO,
157 double secondsElapsed,
158 Opm::RestartValue restartValue,
159 bool writeDoublePrecision)
160 : actionState_(actionState)
161 , wtestState_(wtestState)
162 , summaryState_(summaryState)
163 , udqState_(udqState)
165 , reportStepNum_(reportStepNum)
166 , isSubStep_(isSubStep)
167 , secondsElapsed_(secondsElapsed)
168 , restartValue_(std::move(restartValue))
169 , writeDoublePrecision_(writeDoublePrecision)
175 this->eclIO_.writeTimeStep(this->actionState_,
179 this->reportStepNum_,
181 this->secondsElapsed_,
182 std::move(this->restartValue_),
183 this->writeDoublePrecision_);
191template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
194 const EclipseState& eclState,
195 const SummaryConfig& summaryConfig,
197 const EquilGrid* equilGrid,
198 const GridView& gridView,
201 bool enableAsyncOutput,
203 : collectOnIORank_(grid,
208 summaryConfig.fip_regions_interreg_flow())
210 , gridView_ (gridView)
211 , schedule_ (schedule)
212 , eclState_ (eclState)
213 , cartMapper_ (cartMapper)
214 , equilCartMapper_(equilCartMapper)
215 , equilGrid_ (equilGrid)
218 this->
eclIO_ = std::make_unique<EclipseIO>
220 UgGridHelpers::createEclipseGrid(*equilGrid,
eclState_.getInputGrid()),
221 this->schedule_, summaryConfig,
"", enableEsmry);
226 int numWorkerThreads = 0;
228 numWorkerThreads = 1;
234template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
242template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
244writeInit(
const std::function<
unsigned int(
unsigned int)>& map)
246 if (collectOnIORank_.isIORank()) {
247 std::map<std::string, std::vector<int>> integerVectors;
248 if (collectOnIORank_.isParallel()) {
249 integerVectors.emplace(
"MPI_RANK", collectOnIORank_.globalRanks());
252 auto cartMap = cartesianToCompressed(equilGrid_->size(0), UgGridHelpers::globalCell(*equilGrid_));
254 eclIO_->writeInitial(computeTrans_(cartMap, map),
256 exportNncStructure_(cartMap, map));
260 if (collectOnIORank_.isParallel()) {
261 const auto& comm = grid_.comm();
268template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
271computeTrans_(
const std::unordered_map<int,int>& cartesianToActive,
272 const std::function<
unsigned int(
unsigned int)>& map)
const
274 const auto& cartMapper = *equilCartMapper_;
275 const auto& cartDims = cartMapper.cartesianDimensions();
277 auto tranx = data::CellData {
278 UnitSystem::measure::transmissibility,
279 std::vector<double>(cartDims[0] * cartDims[1] * cartDims[2], 0.0),
280 data::TargetType::INIT
286 using GlobalGridView =
typename EquilGrid::LeafGridView;
287 using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
288 const GlobalGridView& globalGridView = this->equilGrid_->leafGridView();
289 const GlobElementMapper globalElemMapper { globalGridView, Dune::mcmgElementLayout() };
291 auto isNumAquCell = [numAquCell = this->eclState_.aquifer().hasNumericalAquifer()
292 ? this->eclState_.aquifer().numericalAquifers().allAquiferCellIds()
293 : std::vector<std::size_t>{}]
294 (
const std::size_t cellIdx)
296 return std::binary_search(numAquCell.begin(), numAquCell.end(), cellIdx);
299 for (
const auto& elem : elements(globalGridView)) {
300 for (
const auto& is : intersections(globalGridView, elem)) {
305 unsigned c1 = globalElemMapper.index(is.inside());
306 unsigned c2 = globalElemMapper.index(is.outside());
312 const int cartIdx1 = cartMapper.cartesianIndex( c1 );
313 const int cartIdx2 = cartMapper.cartesianIndex( c2 );
315 if (isNumAquCell(cartIdx1) || isNumAquCell(cartIdx2)) {
324 assert(cartIdx1 <= cartIdx2);
325 int gc1 = std::min(cartIdx1, cartIdx2);
326 int gc2 = std::max(cartIdx1, cartIdx2);
334 if (gc2 - gc1 == 1 && cartDims[0] > 1 ) {
335 tranx.data<
double>()[gc1] = globalTrans().transmissibility(c1, c2);
339 if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) {
340 trany.data<
double>()[gc1] = globalTrans().transmissibility(c1, c2);
344 if ( gc2 - gc1 == cartDims[0]*cartDims[1] ||
345 directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2))
346 tranz.data<
double>()[gc1] = globalTrans().transmissibility(c1, c2);
357template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
359EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
360exportNncStructure_(
const std::unordered_map<int,int>& cartesianToActive,
361 const std::function<
unsigned int(
unsigned int)>& map)
const
363 auto isNumAquCell = [numAquCell = this->eclState_.aquifer().hasNumericalAquifer()
364 ? this->eclState_.aquifer().numericalAquifers().allAquiferCellIds()
365 : std::vector<std::size_t>{}]
366 (
const std::size_t cellIdx)
368 return std::binary_search(numAquCell.begin(), numAquCell.end(), cellIdx);
371 auto isNumAquConn = [&isNumAquCell](
const std::size_t cellIdx1,
372 const std::size_t cellIdx2)
374 return isNumAquCell(cellIdx1) || isNumAquCell(cellIdx2);
377 auto isCartesianNeighbour = [nx = this->eclState_.getInputGrid().getNX(),
378 ny = this->eclState_.getInputGrid().getNY()]
379 (
const std::size_t cellIdx1,
const std::size_t cellIdx2)
381 const auto cellDiff = cellIdx2 - cellIdx1;
383 return (cellDiff == 1)
385 || (cellDiff == nx * ny);
388 auto activeCell = [&cartesianToActive](
const std::size_t cellIdx)
390 auto pos = cartesianToActive.find(cellIdx);
391 return (pos == cartesianToActive.end()) ? -1 : pos->second;
394 const auto& nncData = this->eclState_.getInputNNC().input();
395 const auto& unitSystem = this->eclState_.getDeckUnitSystem();
397 for (
const auto& entry : nncData) {
406 assert (entry.cell2 >= entry.cell1);
408 if (! isCartesianNeighbour(entry.cell1, entry.cell2) ||
409 isNumAquConn(entry.cell1, entry.cell2))
414 const auto c1 = activeCell(entry.cell1);
415 const auto c2 = activeCell(entry.cell2);
417 if ((c1 < 0) || (c2 < 0)) {
423 const auto trans = this->globalTrans().transmissibility(c1, c2);
424 const auto tt = unitSystem
425 .from_si(UnitSystem::measure::transmissibility, trans);
430 if (std::isnormal(tt) && ! (tt < 1.0e-6)) {
431 this->outputNnc_.emplace_back(entry.cell1, entry.cell2, trans);
436 auto isDirectNeighbours = [&isCartesianNeighbour, &cartesianToActive,
437 cartDims = &this->cartMapper_.cartesianDimensions()]
438 (
const std::size_t cellIdx1,
const std::size_t cellIdx2)
440 return isCartesianNeighbour(cellIdx1, cellIdx2)
441 || directVerticalNeighbors(*cartDims, cartesianToActive, cellIdx1, cellIdx2);
444 using GlobalGridView =
typename EquilGrid::LeafGridView;
445 using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
446 const GlobalGridView& globalGridView = this->equilGrid_->leafGridView();
447 const GlobElementMapper globalElemMapper { globalGridView, Dune::mcmgElementLayout() };
450 const auto& equilCartMapper = *equilCartMapper_;
451 for (
const auto& elem : elements(globalGridView)) {
452 for (
const auto& is : intersections(globalGridView, elem)) {
457 unsigned c1 = globalElemMapper.index(is.inside());
458 unsigned c2 = globalElemMapper.index(is.outside());
463 std::size_t cc1 = equilCartMapper.cartesianIndex( c1 );
464 std::size_t cc2 = equilCartMapper.cartesianIndex( c2 );
475 if (isNumAquConn(cc1, cc2) || ! isDirectNeighbours(cc1, cc2)) {
478 auto t = this->globalTrans().transmissibility(c1, c2);
479 auto candidate = std::lower_bound(nncData.begin(), nncData.end(),
480 NNCdata { cc1, cc2, 0.0 });
482 while ((candidate != nncData.end()) &&
483 (candidate->cell1 == cc1) &&
484 (candidate->cell2 == cc2))
486 t -= candidate->trans;
495 const auto tt = unitSystem
496 .from_si(UnitSystem::measure::transmissibility, t);
498 if (std::isnormal(tt) && (tt > 1.0e-12)) {
499 this->outputNnc_.emplace_back(cc1, cc2, t);
505 return this->outputNnc_;
508template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
511 const bool isSubStep,
512 data::Solution&& localCellData,
513 data::Wells&& localWellData,
514 data::GroupAndNetworkValues&& localGroupAndNetworkData,
515 data::Aquifers&& localAquiferData,
516 WellTestState&& localWTestState,
517 const Action::State& actionState,
518 const UDQState& udqState,
519 const SummaryState& summaryState,
520 const std::vector<Scalar>& thresholdPressure,
523 bool doublePrecision,
529 const auto isParallel = this->collectOnIORank_.isParallel();
530 const bool needsReordering = this->collectOnIORank_.doesNeedReordering();
532 RestartValue restartValue {
533 (isParallel || needsReordering)
534 ? this->collectOnIORank_.globalCellData()
535 : std::move(localCellData),
537 isParallel ? this->collectOnIORank_.globalWellData()
538 : std::move(localWellData),
540 isParallel ? this->collectOnIORank_.globalGroupAndNetworkData()
541 : std::move(localGroupAndNetworkData),
543 isParallel ? this->collectOnIORank_.globalAquiferData()
544 : std::move(localAquiferData)
547 if (eclState_.getSimulationConfig().useThresholdPressure()) {
548 restartValue.addExtra(
"THRESHPR", UnitSystem::measure::pressure,
554 restartValue.addExtra(
"OPMEXTRA", std::vector<double>(1, nextStepSize));
559 const auto flowsn_global = isParallel ? this->collectOnIORank_.globalFlowsn() : std::move(flowsn);
560 for (
const auto& flows : flowsn_global) {
561 if (flows.name.empty())
563 if (flows.name ==
"FLOGASN+") {
564 restartValue.addExtra(flows.name, UnitSystem::measure::gas_surface_rate, flows.values);
566 restartValue.addExtra(flows.name, UnitSystem::measure::liquid_surface_rate, flows.values);
571 const auto floresn_global = isParallel ? this->collectOnIORank_.globalFloresn() : std::move(floresn);
572 for (
const auto& flores : floresn_global) {
573 if (flores.name.empty()) {
576 restartValue.addExtra(flores.name, UnitSystem::measure::rate, flores.values);
582 auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(
584 isParallel ? this->collectOnIORank_.globalWellTestState() : std::move(localWTestState),
585 summaryState, udqState, *this->eclIO_,
586 reportStepNum, isSubStep, curTime, std::move(restartValue), doublePrecision);
591 this->taskletRunner_->barrier();
594 this->taskletRunner_->dispatch(std::move(eclWriteTasklet));
597template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
600 const Scalar curTime,
601 const data::Wells& localWellData,
602 const data::WellBlockAveragePressures& localWBPData,
603 const data::GroupAndNetworkValues& localGroupAndNetworkData,
604 const std::map<int,data::AquiferData>& localAquiferData,
605 const std::map<std::pair<std::string, int>,
double>& blockData,
606 const std::map<std::string, double>& miscSummaryData,
607 const std::map<std::string, std::vector<double>>& regionData,
608 const Inplace& inplace,
609 const Inplace& initialInPlace,
611 SummaryState& summaryState,
614 if (collectOnIORank_.isIORank()) {
615 const auto& summary = eclIO_->summary();
617 const auto& wellData = this->collectOnIORank_.isParallel()
618 ? this->collectOnIORank_.globalWellData()
621 const auto& wbpData = this->collectOnIORank_.isParallel()
622 ? this->collectOnIORank_.globalWBPData()
625 const auto& groupAndNetworkData = this->collectOnIORank_.isParallel()
626 ? this->collectOnIORank_.globalGroupAndNetworkData()
627 : localGroupAndNetworkData;
629 const auto& aquiferData = this->collectOnIORank_.isParallel()
630 ? this->collectOnIORank_.globalAquiferData()
633 summary.eval(summaryState,
645 getInterRegFlowsAsMap(interRegFlows));
651 const auto udq_step = reportStepNum - 1;
653 this->schedule_.getUDQConfig(udq_step)
656 this->schedule_.wellMatcher(udq_step),
657 this->schedule_.segmentMatcherFactory(udq_step),
658 [es = std::cref(this->eclState_)]() {
659 return std::make_unique<RegionSetMatcher>
660 (es.get().fipRegionStatistics());
662 summaryState, udqState);
666 if (collectOnIORank_.isParallel()) {
673template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
678 assert (globalTrans_);
679 return *globalTrans_;
Definition: CollectDataOnIORank.hpp:49
bool isIORank() const
Definition: CollectDataOnIORank.hpp:125
Definition: EclGenericWriter.hpp:65
void writeInit(const std::function< unsigned int(unsigned int)> &map)
Definition: EclGenericWriter_impl.hpp:244
const EclipseState & eclState_
Definition: EclGenericWriter.hpp:156
CollectDataOnIORankType collectOnIORank_
Definition: EclGenericWriter.hpp:152
void doWriteOutput(const int reportStepNum, const bool isSubStep, data::Solution &&localCellData, data::Wells &&localWellData, data::GroupAndNetworkValues &&localGroupAndNetworkData, data::Aquifers &&localAquiferData, WellTestState &&localWTestState, const Action::State &actionState, const UDQState &udqState, const SummaryState &summaryState, const std::vector< Scalar > &thresholdPressure, Scalar curTime, Scalar nextStepSize, bool doublePrecision, bool isFlowsn, std::array< FlowsData< double >, 3 > &&flowsn, bool isFloresn, std::array< FlowsData< double >, 3 > &&floresn)
Definition: EclGenericWriter_impl.hpp:510
std::unique_ptr< TaskletRunner > taskletRunner_
Definition: EclGenericWriter.hpp:158
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 > > ®ionData, const Inplace &inplace, const Inplace &initialInPlace, const InterRegFlowMap &interRegFlows, SummaryState &summaryState, UDQState &udqState)
Definition: EclGenericWriter_impl.hpp:599
EclGenericWriter(const Schedule &schedule, const EclipseState &eclState, const SummaryConfig &summaryConfig, const Grid &grid, const EquilGrid *equilGrid, const GridView &gridView, const Dune::CartesianIndexMapper< Grid > &cartMapper, const Dune::CartesianIndexMapper< EquilGrid > *equilCartMapper, bool enableAsyncOutput, bool enableEsmry)
Definition: EclGenericWriter_impl.hpp:193
const EclipseIO & eclIO() const
Definition: EclGenericWriter_impl.hpp:236
const TransmissibilityType & globalTrans() const
Definition: EclGenericWriter_impl.hpp:676
std::unique_ptr< EclipseIO > eclIO_
Definition: EclGenericWriter.hpp:157
Inter-region flow accumulation maps for all region definition arrays.
Definition: InterRegFlows.hpp:179
std::vector< data::InterRegFlowMap > getInterRegFlows() const
const std::vector< std::string > & names() const
Class for serializing and broadcasting data using MPI.
Definition: MPISerializer.hpp:31
void append(T &data, int root=0)
Serialize and broadcast on root process, de-serialize and append on others.
Definition: MPISerializer.hpp:107
void broadcast(T &data, int root=0)
Serialize and broadcast on root process, de-serialize on others.
Definition: MPISerializer.hpp:46
Definition: Transmissibility.hpp:54
Definition: BlackoilPhases.hpp:27