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/RPTConfig.hpp>
37#include <opm/input/eclipse/Schedule/Schedule.hpp>
38#include <opm/input/eclipse/Schedule/SummaryState.hpp>
39#include <opm/input/eclipse/Schedule/UDQ/UDQConfig.hpp>
40#include <opm/input/eclipse/Schedule/UDQ/UDQState.hpp>
41#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
42#include <opm/input/eclipse/Schedule/Well/WellMatcher.hpp>
44#include <opm/input/eclipse/Units/UnitSystem.hpp>
46#include <opm/output/eclipse/EclipseIO.hpp>
47#include <opm/output/eclipse/RestartValue.hpp>
48#include <opm/output/eclipse/Summary.hpp>
53#include <opm/simulators/utils/MPISerializer.hpp>
68#include <unordered_map>
88bool directVerticalNeighbors(
const std::array<int, 3>& cartDims,
89 const std::unordered_map<int,int>& cartesianToActive,
90 int smallGlobalIndex,
int largeGlobalIndex)
92 assert(smallGlobalIndex <= largeGlobalIndex);
93 std::array<int, 3> ijk1, ijk2;
94 auto globalToIjk = [cartDims](
int gc) {
95 std::array<int, 3> ijk;
96 ijk[0] = gc % cartDims[0];
98 ijk[1] = gc % cartDims[1];
99 ijk[2] = gc / cartDims[1];
102 ijk1 = globalToIjk(smallGlobalIndex);
103 ijk2 = globalToIjk(largeGlobalIndex);
104 assert(ijk2[2]>=ijk1[2]);
106 if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1)
108 assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0);
109 for (
int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex;
110 gi += cartDims[0] * cartDims[1] )
112 if ( cartesianToActive.find( gi ) != cartesianToActive.end() )
122std::unordered_map<std::string, Opm::data::InterRegFlowMap>
125 auto maps = std::unordered_map<std::string, Opm::data::InterRegFlowMap>{};
127 const auto& regionNames = map.
names();
129 const auto nmap = regionNames.size();
132 for (
auto mapID = 0*nmap; mapID < nmap; ++mapID) {
133 maps.emplace(regionNames[mapID], std::move(flows[mapID]));
141 Opm::Action::State actionState_;
147 std::optional<int> timeStepNum_;
149 double secondsElapsed_;
151 bool writeDoublePrecision_;
153 explicit EclWriteTasklet(
const Opm::Action::State& actionState,
159 std::optional<int> timeStepNum,
161 double secondsElapsed,
163 bool writeDoublePrecision)
164 : actionState_(actionState)
165 , wtestState_(wtestState)
166 , summaryState_(summaryState)
167 , udqState_(udqState)
169 , reportStepNum_(reportStepNum)
170 , timeStepNum_(timeStepNum)
171 , isSubStep_(isSubStep)
172 , secondsElapsed_(secondsElapsed)
173 , restartValue_(std::move(restartValue))
174 , writeDoublePrecision_(writeDoublePrecision)
180 this->eclIO_.writeTimeStep(this->actionState_,
184 this->reportStepNum_,
186 this->secondsElapsed_,
187 std::move(this->restartValue_),
188 this->writeDoublePrecision_,
198template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
199EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
200EclGenericWriter(
const Schedule& schedule,
201 const EclipseState& eclState,
204 const EquilGrid* equilGrid,
205 const GridView& gridView,
210 : collectOnIORank_(grid,
217 , gridView_ (gridView)
218 , schedule_ (schedule)
219 , eclState_ (eclState)
222 , equilGrid_ (equilGrid)
224 if (this->collectOnIORank_.isIORank()) {
225 this->eclIO_ = std::make_unique<EclipseIO>
227 UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()),
228 this->schedule_, summaryConfig,
"", enableEsmry);
233 int numWorkerThreads = 0;
234 if (enableAsyncOutput && collectOnIORank_.isIORank()) {
235 numWorkerThreads = 1;
238 this->taskletRunner_.reset(
new TaskletRunner(numWorkerThreads));
241template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
242const EclipseIO& EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
249template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
250void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
253 if (collectOnIORank_.isIORank()) {
254 std::map<std::string, std::vector<int>> integerVectors;
255 if (collectOnIORank_.isParallel()) {
256 integerVectors.emplace(
"MPI_RANK", collectOnIORank_.globalRanks());
259 eclIO_->writeInitial(*this->outputTrans_,
262 this->outputTrans_.reset();
265template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
267EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
268extractOutputTransAndNNC(
const std::function<
unsigned int(
unsigned int)>& map)
270 if (collectOnIORank_.isIORank()) {
271 auto cartMap = cartesianToCompressed(equilGrid_->size(0), UgGridHelpers::globalCell(*equilGrid_));
272 computeTrans_(cartMap, map);
273 exportNncStructure_(cartMap, map);
277 if (collectOnIORank_.isParallel()) {
278 const auto& comm = grid_.comm();
280 ser.broadcast(outputNnc_);
285template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
287EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
288computeTrans_(
const std::unordered_map<int,int>& cartesianToActive,
289 const std::function<
unsigned int(
unsigned int)>& map)
const
292 outputTrans_ = std::make_unique<data::Solution>();
295 const auto& cartMapper = *equilCartMapper_;
296 const auto& cartDims = cartMapper.cartesianDimensions();
298 auto createCellData = [&cartDims]() {
299 return data::CellData{
300 UnitSystem::measure::transmissibility,
301 std::vector<double>(cartDims[0] * cartDims[1] * cartDims[2], 0.0),
302 data::TargetType::INIT
306 outputTrans_->clear();
307 outputTrans_->emplace(
"TRANX", createCellData());
308 outputTrans_->emplace(
"TRANY", createCellData());
309 outputTrans_->emplace(
"TRANZ", createCellData());
311 auto& tranx = this->outputTrans_->at(
"TRANX");
312 auto& trany = this->outputTrans_->at(
"TRANY");
313 auto& tranz = this->outputTrans_->at(
"TRANZ");
315 using GlobalGridView =
typename EquilGrid::LeafGridView;
316 using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
317 const GlobalGridView& globalGridView = this->equilGrid_->leafGridView();
318 const GlobElementMapper globalElemMapper { globalGridView, Dune::mcmgElementLayout() };
320 auto isNumAquCell = [numAquCell = this->eclState_.aquifer().hasNumericalAquifer()
321 ? this->eclState_.aquifer().numericalAquifers().allAquiferCellIds()
322 : std::vector<std::size_t>{}]
323 (
const std::size_t cellIdx)
325 return std::binary_search(numAquCell.begin(), numAquCell.end(), cellIdx);
328 for (
const auto& elem : elements(globalGridView)) {
329 for (
const auto& is : intersections(globalGridView, elem)) {
334 unsigned c1 = globalElemMapper.index(is.inside());
335 unsigned c2 = globalElemMapper.index(is.outside());
341 const int cartIdx1 = cartMapper.cartesianIndex( c1 );
342 const int cartIdx2 = cartMapper.cartesianIndex( c2 );
344 if (isNumAquCell(cartIdx1) || isNumAquCell(cartIdx2)) {
353 assert(cartIdx1 <= cartIdx2);
354 int gc1 = std::min(cartIdx1, cartIdx2);
355 int gc2 = std::max(cartIdx1, cartIdx2);
363 if (gc2 - gc1 == 1 && cartDims[0] > 1 ) {
364 tranx.template data<double>()[gc1] = globalTrans().transmissibility(c1, c2);
368 if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) {
369 trany.template data<double>()[gc1] = globalTrans().transmissibility(c1, c2);
373 if ( gc2 - gc1 == cartDims[0]*cartDims[1] ||
374 directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2))
375 tranz.template data<double>()[gc1] = globalTrans().transmissibility(c1, c2);
380template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
382EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
383exportNncStructure_(
const std::unordered_map<int,int>& cartesianToActive,
384 const std::function<
unsigned int(
unsigned int)>& map)
const
386 auto isNumAquCell = [numAquCell = this->eclState_.aquifer().hasNumericalAquifer()
387 ? this->eclState_.aquifer().numericalAquifers().allAquiferCellIds()
388 : std::vector<std::size_t>{}]
389 (
const std::size_t cellIdx)
391 return std::binary_search(numAquCell.begin(), numAquCell.end(), cellIdx);
394 auto isNumAquConn = [&isNumAquCell](
const std::size_t cellIdx1,
395 const std::size_t cellIdx2)
397 return isNumAquCell(cellIdx1) || isNumAquCell(cellIdx2);
400 auto isCartesianNeighbour = [nx = this->eclState_.getInputGrid().getNX(),
401 ny = this->eclState_.getInputGrid().getNY()]
402 (
const std::size_t cellIdx1,
const std::size_t cellIdx2)
404 const auto cellDiff = cellIdx2 - cellIdx1;
406 return (cellDiff == 1)
408 || (cellDiff == nx * ny);
411 auto activeCell = [&cartesianToActive](
const std::size_t cellIdx)
413 auto pos = cartesianToActive.find(cellIdx);
414 return (pos == cartesianToActive.end()) ? -1 : pos->second;
417 const auto& nncData = this->eclState_.getInputNNC().input();
418 const auto& unitSystem = this->eclState_.getDeckUnitSystem();
420 for (
const auto& entry : nncData) {
429 assert (entry.cell2 >= entry.cell1);
431 if (! isCartesianNeighbour(entry.cell1, entry.cell2) ||
432 isNumAquConn(entry.cell1, entry.cell2))
437 const auto c1 = activeCell(entry.cell1);
438 const auto c2 = activeCell(entry.cell2);
440 if ((c1 < 0) || (c2 < 0)) {
446 const auto trans = this->globalTrans().transmissibility(c1, c2);
447 const auto tt = unitSystem
448 .from_si(UnitSystem::measure::transmissibility, trans);
453 if (std::isnormal(tt) && ! (tt < 1.0e-6)) {
454 this->outputNnc_.emplace_back(entry.cell1, entry.cell2, trans);
459 auto isDirectNeighbours = [&isCartesianNeighbour, &cartesianToActive,
460 cartDims = &this->cartMapper_.cartesianDimensions()]
461 (
const std::size_t cellIdx1,
const std::size_t cellIdx2)
463 return isCartesianNeighbour(cellIdx1, cellIdx2)
464 || directVerticalNeighbors(*cartDims, cartesianToActive, cellIdx1, cellIdx2);
467 using GlobalGridView =
typename EquilGrid::LeafGridView;
468 using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
469 const GlobalGridView& globalGridView = this->equilGrid_->leafGridView();
470 const GlobElementMapper globalElemMapper { globalGridView, Dune::mcmgElementLayout() };
473 const auto& equilCartMapper = *equilCartMapper_;
474 for (
const auto& elem : elements(globalGridView)) {
475 for (
const auto& is : intersections(globalGridView, elem)) {
480 unsigned c1 = globalElemMapper.index(is.inside());
481 unsigned c2 = globalElemMapper.index(is.outside());
486 std::size_t cc1 = equilCartMapper.cartesianIndex( c1 );
487 std::size_t cc2 = equilCartMapper.cartesianIndex( c2 );
498 if (isNumAquConn(cc1, cc2) || ! isDirectNeighbours(cc1, cc2)) {
501 auto t = this->globalTrans().transmissibility(c1, c2);
502 auto candidate = std::lower_bound(nncData.begin(), nncData.end(),
503 NNCdata { cc1, cc2, 0.0 });
505 while ((candidate != nncData.end()) &&
506 (candidate->cell1 == cc1) &&
507 (candidate->cell2 == cc2))
509 t -= candidate->trans;
518 const auto tt = unitSystem
519 .from_si(UnitSystem::measure::transmissibility, t);
521 if (std::isnormal(tt) && (tt > 1.0e-12)) {
522 this->outputNnc_.emplace_back(cc1, cc2, t);
528 return this->outputNnc_;
531template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
532void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
533doWriteOutput(
const int reportStepNum,
534 const std::optional<int> timeStepNum,
535 const bool isSubStep,
536 data::Solution&& localCellData,
537 data::Wells&& localWellData,
538 data::GroupAndNetworkValues&& localGroupAndNetworkData,
539 data::Aquifers&& localAquiferData,
540 WellTestState&& localWTestState,
541 const Action::State& actionState,
542 const UDQState& udqState,
543 const SummaryState& summaryState,
544 const std::vector<Scalar>& thresholdPressure,
547 bool doublePrecision,
549 std::array<FlowsData<double>, 3>&& flowsn,
551 std::array<FlowsData<double>, 3>&& floresn)
553 const auto isParallel = this->collectOnIORank_.isParallel();
554 const bool needsReordering = this->collectOnIORank_.doesNeedReordering();
556 RestartValue restartValue {
557 (isParallel || needsReordering)
558 ? this->collectOnIORank_.globalCellData()
559 : std::move(localCellData),
561 isParallel ? this->collectOnIORank_.globalWellData()
562 : std::move(localWellData),
564 isParallel ? this->collectOnIORank_.globalGroupAndNetworkData()
565 : std::move(localGroupAndNetworkData),
567 isParallel ? this->collectOnIORank_.globalAquiferData()
568 : std::move(localAquiferData)
571 if (eclState_.getSimulationConfig().useThresholdPressure()) {
572 restartValue.addExtra(
"THRESHPR", UnitSystem::measure::pressure,
578 restartValue.addExtra(
"OPMEXTRA", std::vector<double>(1, nextStepSize));
583 const auto flowsn_global = isParallel ? this->collectOnIORank_.globalFlowsn() : std::move(flowsn);
584 for (
const auto& flows : flowsn_global) {
585 if (flows.name.empty())
587 if (flows.name ==
"FLOGASN+") {
588 restartValue.addExtra(flows.name, UnitSystem::measure::gas_surface_rate, flows.values);
590 restartValue.addExtra(flows.name, UnitSystem::measure::liquid_surface_rate, flows.values);
595 const auto floresn_global = isParallel ? this->collectOnIORank_.globalFloresn() : std::move(floresn);
596 for (
const auto& flores : floresn_global) {
597 if (flores.name.empty()) {
600 restartValue.addExtra(flores.name, UnitSystem::measure::rate, flores.values);
606 this->taskletRunner_->barrier();
609 if (this->taskletRunner_->failure()) {
610 throw std::runtime_error(
"Failure in the TaskletRunner while writing output.");
614 auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(
616 isParallel ? this->collectOnIORank_.globalWellTestState() : std::move(localWTestState),
617 summaryState, udqState, *this->eclIO_,
618 reportStepNum, timeStepNum, isSubStep, curTime, std::move(restartValue), doublePrecision);
621 this->taskletRunner_->dispatch(std::move(eclWriteTasklet));
624template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
625void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
626evalSummary(
const int reportStepNum,
627 const Scalar curTime,
628 const data::Wells& localWellData,
629 const data::WellBlockAveragePressures& localWBPData,
630 const data::GroupAndNetworkValues& localGroupAndNetworkData,
631 const std::map<int,data::AquiferData>& localAquiferData,
632 const std::map<std::pair<std::string, int>,
double>& blockData,
633 const std::map<std::string, double>& miscSummaryData,
634 const std::map<std::string, std::vector<double>>& regionData,
635 const Inplace& inplace,
636 const Inplace& initialInPlace,
637 const InterRegFlowMap& interRegFlows,
638 SummaryState& summaryState,
641 if (collectOnIORank_.isIORank()) {
642 const auto& summary = eclIO_->summary();
644 const auto& wellData = this->collectOnIORank_.isParallel()
645 ? this->collectOnIORank_.globalWellData()
648 const auto& wbpData = this->collectOnIORank_.isParallel()
649 ? this->collectOnIORank_.globalWBPData()
652 const auto& groupAndNetworkData = this->collectOnIORank_.isParallel()
653 ? this->collectOnIORank_.globalGroupAndNetworkData()
654 : localGroupAndNetworkData;
656 const auto& aquiferData = this->collectOnIORank_.isParallel()
657 ? this->collectOnIORank_.globalAquiferData()
660 summary.eval(summaryState,
672 getInterRegFlowsAsMap(interRegFlows));
678 const auto udq_step = reportStepNum - 1;
680 this->schedule_[udq_step].udq()
682 this->schedule_.wellMatcher(udq_step),
683 this->schedule_[udq_step].group_order(),
684 this->schedule_.segmentMatcherFactory(udq_step),
685 [es = std::cref(this->eclState_)]() {
686 return std::make_unique<RegionSetMatcher>
687 (es.get().fipRegionStatistics());
689 summaryState, udqState);
693 if (collectOnIORank_.isParallel()) {
694 Parallel::MpiSerializer ser(grid_.comm());
695 ser.append(summaryState);
700template<
class Gr
id,
class EquilGr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
701const typename EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::TransmissibilityType&
702EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
705 assert (globalTrans_);
706 return *globalTrans_;
Collects necessary output values and pass it to opm-common's ECL output.
Definition CollectDataOnIORank.hpp:49
Inter-region flow accumulation maps for all region definition arrays.
Definition InterRegFlows.hpp:179
const std::vector< std::string > & names() const
Names of all applicable region definition arrays.
Definition InterRegFlows.cpp:182
std::vector< data::InterRegFlowMap > getInterRegFlows() const
Get read-only access to the underlying CSR representation.
Definition InterRegFlows.cpp:188
Class for serializing and broadcasting data using MPI.
Definition MPISerializer.hpp:31
The base class for tasklets.
Definition tasklets.hpp:44
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242