EclWriter.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef OPM_ECL_WRITER_HPP
29#define OPM_ECL_WRITER_HPP
30
31#include <dune/grid/common/partitionset.hh>
32
33#include <opm/common/TimingMacros.hpp> // OPM_TIMEBLOCK
34#include <opm/common/OpmLog/OpmLog.hpp>
35#include <opm/input/eclipse/Schedule/RPTConfig.hpp>
36
37#include <opm/input/eclipse/Units/UnitSystem.hpp>
38#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
39
40#include <opm/output/eclipse/Inplace.hpp>
41#include <opm/output/eclipse/RestartValue.hpp>
42
43#include <opm/models/blackoil/blackoilproperties.hh> // Properties::EnableMech, EnableTemperature, EnableSolvent
44#include <opm/models/common/multiphasebaseproperties.hh> // Properties::FluidSystem
45
54
55#include <boost/date_time/posix_time/posix_time.hpp>
56
57#include <limits>
58#include <map>
59#include <memory>
60#include <optional>
61#include <stdexcept>
62#include <string>
63#include <utility>
64#include <vector>
65
66namespace Opm::Parameters {
67
68// If available, write the ECL output in a non-blocking manner
69struct EnableAsyncEclOutput { static constexpr bool value = true; };
70
71// By default, use single precision for the ECL formated results
72struct EclOutputDoublePrecision { static constexpr bool value = false; };
73
74// Write all solutions for visualization, not just the ones for the
75// report steps...
76struct EnableWriteAllSolutions { static constexpr bool value = false; };
77
78// Write ESMRY file for fast loading of summary data
79struct EnableEsmry { static constexpr bool value = false; };
80
81} // namespace Opm::Parameters
82
83namespace Opm::Action {
84 class State;
85} // namespace Opm::Action
86
87namespace Opm {
88 class EclipseIO;
89 class UDQState;
90} // namespace Opm
91
92namespace Opm {
108template <class TypeTag, class OutputModule>
109class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>,
110 GetPropType<TypeTag, Properties::EquilGrid>,
111 GetPropType<TypeTag, Properties::GridView>,
112 GetPropType<TypeTag, Properties::ElementMapper>,
113 GetPropType<TypeTag, Properties::Scalar>>
114{
123 using Element = typename GridView::template Codim<0>::Entity;
125 using ElementIterator = typename GridView::template Codim<0>::Iterator;
127
128 typedef Dune::MultipleCodimMultipleGeomTypeMapper< GridView > VertexMapper;
129
130 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
131 enum { enableMech = getPropValue<TypeTag, Properties::EnableMech>() };
132 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
133 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
134
135public:
136
137 static void registerParameters()
138 {
139 OutputModule::registerParameters();
140
141 Parameters::Register<Parameters::EnableAsyncEclOutput>
142 ("Write the ECL-formated results in a non-blocking way "
143 "(i.e., using a separate thread).");
144 Parameters::Register<Parameters::EnableEsmry>
145 ("Write ESMRY file for fast loading of summary data.");
146 }
147
148 // The Simulator object should preferably have been const - the
149 // only reason that is not the case is due to the SummaryState
150 // object owned deep down by the vanguard.
151 explicit EclWriter(Simulator& simulator)
152 : BaseType(simulator.vanguard().schedule(),
153 simulator.vanguard().eclState(),
154 simulator.vanguard().summaryConfig(),
155 simulator.vanguard().grid(),
156 ((simulator.vanguard().grid().comm().rank() == 0)
157 ? &simulator.vanguard().equilGrid()
158 : nullptr),
159 simulator.vanguard().gridView(),
160 simulator.vanguard().cartesianIndexMapper(),
161 ((simulator.vanguard().grid().comm().rank() == 0)
162 ? &simulator.vanguard().equilCartesianIndexMapper()
163 : nullptr),
164 Parameters::Get<Parameters::EnableAsyncEclOutput>(),
165 Parameters::Get<Parameters::EnableEsmry>())
166 , simulator_(simulator)
167 {
168#if HAVE_MPI
169 if (this->simulator_.vanguard().grid().comm().size() > 1) {
170 auto smryCfg = (this->simulator_.vanguard().grid().comm().rank() == 0)
171 ? this->eclIO_->finalSummaryConfig()
172 : SummaryConfig{};
173
174 eclBroadcast(this->simulator_.vanguard().grid().comm(), smryCfg);
175
176 this->outputModule_ = std::make_unique<OutputModule>
177 (simulator, smryCfg, this->collectOnIORank_);
178 }
179 else
180#endif
181 {
182 this->outputModule_ = std::make_unique<OutputModule>
183 (simulator, this->eclIO_->finalSummaryConfig(), this->collectOnIORank_);
184 }
185
186 this->rank_ = this->simulator_.vanguard().grid().comm().rank();
187
188 this->simulator_.vanguard().eclState().computeFipRegionStatistics();
189 }
190
192 {}
193
194 const EquilGrid& globalGrid() const
195 {
196 return simulator_.vanguard().equilGrid();
197 }
198
202 void evalSummaryState(bool isSubStep)
203 {
204 OPM_TIMEBLOCK(evalSummaryState);
205 const int reportStepNum = simulator_.episodeIndex() + 1;
206
207 /*
208 The summary data is not evaluated for timestep 0, that is
209 implemented with a:
210
211 if (time_step == 0)
212 return;
213
214 check somewhere in the summary code. When the summary code was
215 split in separate methods Summary::eval() and
216 Summary::add_timestep() it was necessary to pull this test out
217 here to ensure that the well and group related keywords in the
218 restart file, like XWEL and XGRP were "correct" also in the
219 initial report step.
220
221 "Correct" in this context means unchanged behavior, might very
222 well be more correct to actually remove this if test.
223 */
224
225 if (reportStepNum == 0)
226 return;
227
228 const Scalar curTime = simulator_.time() + simulator_.timeStepSize();
229 const Scalar totalCpuTime =
230 simulator_.executionTimer().realTimeElapsed() +
231 simulator_.setupTimer().realTimeElapsed() +
232 simulator_.vanguard().setupTime();
233
234 const auto localWellData = simulator_.problem().wellModel().wellData();
235 const auto localWBP = simulator_.problem().wellModel().wellBlockAveragePressures();
236 const auto localGroupAndNetworkData = simulator_.problem().wellModel()
237 .groupAndNetworkData(reportStepNum);
238
239 const auto localAquiferData = simulator_.problem().aquiferModel().aquiferData();
240 const auto localWellTestState = simulator_.problem().wellModel().wellTestState();
241 this->prepareLocalCellData(isSubStep, reportStepNum);
242
243 if (this->outputModule_->needInterfaceFluxes(isSubStep)) {
244 this->captureLocalFluxData();
245 }
246
247 if (this->collectOnIORank_.isParallel()) {
249
250 std::map<std::pair<std::string,int>,double> dummy;
251 this->collectOnIORank_.collect({},
252 outputModule_->getBlockData(),
253 dummy,
254 localWellData,
255 localWBP,
256 localGroupAndNetworkData,
257 localAquiferData,
258 localWellTestState,
259 this->outputModule_->getInterRegFlows(),
260 {},
261 {});
262
263 if (this->collectOnIORank_.isIORank()) {
264 auto& iregFlows = this->collectOnIORank_.globalInterRegFlows();
265
266 if (! iregFlows.readIsConsistent()) {
267 throw std::runtime_error {
268 "Inconsistent inter-region flow "
269 "region set names in parallel"
270 };
271 }
272
273 iregFlows.compress();
274 }
275
276 OPM_END_PARALLEL_TRY_CATCH("Collect to I/O rank: ",
277 this->simulator_.vanguard().grid().comm());
278 }
279
280
281 std::map<std::string, double> miscSummaryData;
282 std::map<std::string, std::vector<double>> regionData;
283 Inplace inplace;
284
285 {
286 OPM_TIMEBLOCK(outputFipLogAndFipresvLog);
287
288 inplace = outputModule_->calc_inplace(miscSummaryData, regionData, simulator_.gridView().comm());
289
290 if (this->collectOnIORank_.isIORank()){
291 inplace_ = inplace;
292 }
293 }
294
295 // Add TCPU
296 if (totalCpuTime != 0.0) {
297 miscSummaryData["TCPU"] = totalCpuTime;
298 }
300 miscSummaryData["NEWTON"] = this->sub_step_report_.total_newton_iterations;
301 }
303 miscSummaryData["MLINEARS"] = this->sub_step_report_.total_linear_iterations;
304 }
306 miscSummaryData["NLINEARS"] = static_cast<float>(this->sub_step_report_.total_linear_iterations) / this->sub_step_report_.total_newton_iterations;
307 }
308 if (this->sub_step_report_.min_linear_iterations != std::numeric_limits<unsigned int>::max()) {
309 miscSummaryData["NLINSMIN"] = this->sub_step_report_.min_linear_iterations;
310 }
312 miscSummaryData["NLINSMAX"] = this->sub_step_report_.max_linear_iterations;
313 }
315 miscSummaryData["MSUMLINS"] = this->simulation_report_.success.total_linear_iterations;
316 }
318 miscSummaryData["MSUMNEWT"] = this->simulation_report_.success.total_newton_iterations;
319 }
320
321 {
322 OPM_TIMEBLOCK(evalSummary);
323
324 const auto& blockData = this->collectOnIORank_.isParallel()
326 : this->outputModule_->getBlockData();
327
328 const auto& interRegFlows = this->collectOnIORank_.isParallel()
330 : this->outputModule_->getInterRegFlows();
331
332 this->evalSummary(reportStepNum,
333 curTime,
334 localWellData,
335 localWBP,
336 localGroupAndNetworkData,
337 localAquiferData,
338 blockData,
339 miscSummaryData,
340 regionData,
341 inplace,
342 this->outputModule_->initialInplace(),
343 interRegFlows,
344 this->summaryState(),
345 this->udqState());
346 }
347 }
348
351 {
352 const auto& gridView = simulator_.vanguard().gridView();
353 const int num_interior = detail::
355
356 this->outputModule_->
357 allocBuffers(num_interior, 0, false, false, /*isRestart*/ false);
358
359#ifdef _OPENMP
360#pragma omp parallel for
361#endif
362 for (int dofIdx = 0; dofIdx < num_interior; ++dofIdx) {
363 const auto& intQuants = *simulator_.model().cachedIntensiveQuantities(dofIdx, /*timeIdx=*/0);
364 const auto totVolume = simulator_.model().dofTotalVolume(dofIdx);
365
366 this->outputModule_->updateFluidInPlace(dofIdx, intQuants, totVolume);
367 }
368
369 // We always calculate the initial fip values as it may be used by various
370 // keywords in the Schedule, e.g. FIP=2 in RPTSCHED but no FIP in RPTSOL
371 outputModule_->calc_initial_inplace(simulator_.gridView().comm());
372
373 // check if RPTSOL entry has FIP output
374 const auto& fip = simulator_.vanguard().eclState().getEclipseConfig().fip();
375 if (fip.output(FIPConfig::OutputField::FIELD) ||
376 fip.output(FIPConfig::OutputField::RESV))
377 {
378 OPM_TIMEBLOCK(outputFipLogAndFipresvLog);
379 boost::posix_time::ptime start_time =
380 boost::posix_time::from_time_t(simulator_.vanguard().schedule().getStartTime());
381
382 if (this->collectOnIORank_.isIORank()) {
383 inplace_ = outputModule_->initialInplace().value();
384 outputModule_->outputFipAndResvLog(inplace_, 0, 0.0, start_time,
385 false, simulator_.gridView().comm());
386 }
387 }
388
389 outputModule_->outputFipAndResvLogToCSV(0, false, simulator_.gridView().comm());
390 }
391
392 void writeReports(const SimulatorTimer& timer)
393 {
394 if (! this->collectOnIORank_.isIORank()) {
395 return;
396 }
397
398 // SimulatorTimer::reportStepNum() is the simulator's zero-based
399 // "episode index". This is generally the index value needed to
400 // look up objects in the Schedule container. That said, function
401 // writeReports() is invoked at the *beginning* of a report
402 // step/episode which means we typically need the objects from the
403 // *previous* report step/episode. We therefore need special case
404 // handling for reportStepNum() == 0 in base runs and
405 // reportStepNum() <= restart step in restarted runs.
406 const auto firstStep = this->initialStep();
407 const auto simStep =
408 std::max(timer.reportStepNum() - 1, firstStep);
409
410 const auto& rpt = this->schedule_[simStep].rpt_config();
411
412 if (rpt.contains("WELSPECS") && (rpt.at("WELSPECS") > 0)) {
413 // Requesting a well specification report is valid at all times,
414 // including reportStepNum() == initialStep().
415 this->writeWellspecReport(timer);
416 }
417
418 if (timer.reportStepNum() == firstStep) {
419 // No dynamic flows at the beginning of the initialStep().
420 return;
421 }
422
423 if (rpt.contains("WELLS") && rpt.at("WELLS") > 0) {
424 this->writeWellflowReport(timer, simStep, rpt.at("WELLS"));
425 }
426
427 this->outputModule_->outputFipAndResvLog(this->inplace_,
428 timer.reportStepNum(),
429 timer.simulationTimeElapsed(),
430 timer.currentDateTime(),
431 /* isSubstep = */ false,
432 simulator_.gridView().comm());
433
434 OpmLog::note(""); // Blank line after all reports.
435 }
436
437 void writeOutput(data::Solution&& localCellData, bool isSubStep)
438 {
439 OPM_TIMEBLOCK(writeOutput);
440
441 const int reportStepNum = simulator_.episodeIndex() + 1;
442 this->prepareLocalCellData(isSubStep, reportStepNum);
443 this->outputModule_->outputErrorLog(simulator_.gridView().comm());
444
445 // output using eclWriter if enabled
446 auto localWellData = simulator_.problem().wellModel().wellData();
447 auto localGroupAndNetworkData = simulator_.problem().wellModel()
448 .groupAndNetworkData(reportStepNum);
449
450 auto localAquiferData = simulator_.problem().aquiferModel().aquiferData();
451 auto localWellTestState = simulator_.problem().wellModel().wellTestState();
452
453 const bool isFlowsn = this->outputModule_->getFlows().hasFlowsn();
454 auto flowsn = this->outputModule_->getFlows().getFlowsn();
455
456 const bool isFloresn = this->outputModule_->getFlows().hasFloresn();
457 auto floresn = this->outputModule_->getFlows().getFloresn();
458
459 if (! isSubStep || Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
460
461 if (localCellData.empty()) {
462 this->outputModule_->assignToSolution(localCellData);
463 }
464
465 // Add cell data to perforations for RFT output
466 this->outputModule_->addRftDataToWells(localWellData,
467 reportStepNum,
468 simulator_.gridView().comm());
469 }
470
471 if (this->collectOnIORank_.isParallel() ||
472 this->collectOnIORank_.doesNeedReordering())
473 {
474 // Note: We don't need WBP (well-block averaged pressures) or
475 // inter-region flow rate values in order to create restart file
476 // output. There's consequently no need to collect those
477 // properties on the I/O rank.
478
479 this->collectOnIORank_.collect(localCellData,
480 this->outputModule_->getBlockData(),
481 this->outputModule_->getExtraBlockData(),
482 localWellData,
483 /* wbpData = */ {},
484 localGroupAndNetworkData,
485 localAquiferData,
486 localWellTestState,
487 /* interRegFlows = */ {},
488 flowsn,
489 floresn);
490 if (this->collectOnIORank_.isIORank()) {
491 this->outputModule_->assignGlobalFieldsToSolution(this->collectOnIORank_.globalCellData());
492 }
493 } else {
494 this->outputModule_->assignGlobalFieldsToSolution(localCellData);
495 }
496
497 if (this->collectOnIORank_.isIORank()) {
498 const Scalar curTime = simulator_.time() + simulator_.timeStepSize();
499 const Scalar nextStepSize = simulator_.problem().nextTimeStepSize();
500 std::optional<int> timeStepIdx;
501 if (Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
502 timeStepIdx = simulator_.timeStepIndex();
503 }
504 this->doWriteOutput(reportStepNum, timeStepIdx, isSubStep,
505 std::move(localCellData),
506 std::move(localWellData),
507 std::move(localGroupAndNetworkData),
508 std::move(localAquiferData),
509 std::move(localWellTestState),
510 this->actionState(),
511 this->udqState(),
512 this->summaryState(),
513 this->simulator_.problem().thresholdPressure().getRestartVector(),
514 curTime, nextStepSize,
515 Parameters::Get<Parameters::EclOutputDoublePrecision>(),
516 isFlowsn, std::move(flowsn),
517 isFloresn, std::move(floresn));
518 }
519 }
520
522 {
523 const auto enablePCHysteresis = simulator_.problem().materialLawManager()->enablePCHysteresis();
524 const auto enableNonWettingHysteresis = simulator_.problem().materialLawManager()->enableNonWettingHysteresis();
525 const auto enableWettingHysteresis = simulator_.problem().materialLawManager()->enableWettingHysteresis();
526 const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
527 const auto gasActive = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
528 const auto waterActive = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
529 const auto enableSwatinit = simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT");
530
531 std::vector<RestartKey> solutionKeys {
532 {"PRESSURE", UnitSystem::measure::pressure},
533 {"SWAT", UnitSystem::measure::identity, waterActive},
534 {"SGAS", UnitSystem::measure::identity, gasActive},
535 {"TEMP", UnitSystem::measure::temperature, enableEnergy},
536 {"SSOLVENT", UnitSystem::measure::identity, enableSolvent},
537
538 {"RS", UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()},
539 {"RV", UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()},
540 {"RVW", UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedWater()},
541 {"RSW", UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGasInWater()},
542
543 {"SGMAX", UnitSystem::measure::identity, enableNonWettingHysteresis && oilActive && gasActive},
544 {"SHMAX", UnitSystem::measure::identity, enableWettingHysteresis && oilActive && gasActive},
545
546 {"SOMAX", UnitSystem::measure::identity,
547 (enableNonWettingHysteresis && oilActive && waterActive)
548 || simulator_.problem().vapparsActive(simulator_.episodeIndex())},
549
550 {"SOMIN", UnitSystem::measure::identity, enablePCHysteresis && oilActive && gasActive},
551 {"SWHY1", UnitSystem::measure::identity, enablePCHysteresis && oilActive && waterActive},
552 {"SWMAX", UnitSystem::measure::identity, enableWettingHysteresis && oilActive && waterActive},
553
554 {"PPCW", UnitSystem::measure::pressure, enableSwatinit},
555 };
556
557 {
558 const auto& tracers = simulator_.vanguard().eclState().tracer();
559
560 for (const auto& tracer : tracers) {
561 const auto enableSolTracer =
562 ((tracer.phase == Phase::GAS) && FluidSystem::enableDissolvedGas()) ||
563 ((tracer.phase == Phase::OIL) && FluidSystem::enableVaporizedOil());
564
565 solutionKeys.emplace_back(tracer.fname(), UnitSystem::measure::identity, true);
566 solutionKeys.emplace_back(tracer.sname(), UnitSystem::measure::identity, enableSolTracer);
567 }
568 }
569
570 const auto& inputThpres = eclState().getSimulationConfig().getThresholdPressure();
571 const std::vector<RestartKey> extraKeys {
572 {"OPMEXTRA", UnitSystem::measure::identity, false},
573 {"THRESHPR", UnitSystem::measure::pressure, inputThpres.active()},
574 };
575
576 const auto& gridView = this->simulator_.vanguard().gridView();
577 const auto numElements = gridView.size(/*codim=*/0);
578
579 // Try to load restart step 0 to calculate initial FIP
580 {
581 this->outputModule_->allocBuffers(numElements,
582 0,
583 /*isSubStep = */false,
584 /*log = */ false,
585 /*isRestart = */true);
586
587 const auto restartSolution =
589 solutionKeys, gridView.comm(), 0);
590
591 if (!restartSolution.empty()) {
592 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
593 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
594 this->outputModule_->setRestart(restartSolution, elemIdx, globalIdx);
595 }
596
597 this->simulator_.problem().readSolutionFromOutputModule(0, true);
598 ElementContext elemCtx(this->simulator_);
599 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
600 elemCtx.updatePrimaryStencil(elem);
601 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
602
603 this->outputModule_->updateFluidInPlace(elemCtx);
604 }
605
606 this->outputModule_->calc_initial_inplace(this->simulator_.gridView().comm());
607 }
608 }
609
610 {
611 // The episodeIndex is rewound one step back before calling
612 // beginRestart() and cannot be used here. We just ask the
613 // initconfig directly to be sure that we use the correct index.
614 const auto restartStepIdx = this->simulator_.vanguard()
615 .eclState().getInitConfig().getRestartStep();
616
617 this->outputModule_->allocBuffers(numElements,
618 restartStepIdx,
619 /*isSubStep = */false,
620 /*log = */ false,
621 /*isRestart = */true);
622 }
623
624 {
625 const auto restartValues =
626 loadParallelRestart(this->eclIO_.get(),
627 this->actionState(),
628 this->summaryState(),
629 solutionKeys, extraKeys, gridView.comm());
630
631 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
632 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
633 this->outputModule_->setRestart(restartValues.solution, elemIdx, globalIdx);
634 }
635
636 auto& tracer_model = simulator_.problem().tracerModel();
637 for (int tracer_index = 0; tracer_index < tracer_model.numTracers(); ++tracer_index) {
638 // Free tracers
639 {
640 const auto& free_tracer_name = tracer_model.fname(tracer_index);
641 const auto& free_tracer_solution = restartValues.solution
642 .template data<double>(free_tracer_name);
643
644 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
645 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
646 tracer_model.setFreeTracerConcentration
647 (tracer_index, elemIdx, free_tracer_solution[globalIdx]);
648 }
649 }
650
651 // Solution tracer (only if DISGAS/VAPOIL are active for gas/oil tracers)
652 if ((tracer_model.phase(tracer_index) == Phase::GAS && FluidSystem::enableDissolvedGas()) ||
653 (tracer_model.phase(tracer_index) == Phase::OIL && FluidSystem::enableVaporizedOil()))
654 {
655 tracer_model.setEnableSolTracers(tracer_index, true);
656
657 const auto& sol_tracer_name = tracer_model.sname(tracer_index);
658 const auto& sol_tracer_solution = restartValues.solution
659 .template data<double>(sol_tracer_name);
660
661 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
662 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
663 tracer_model.setSolTracerConcentration
664 (tracer_index, elemIdx, sol_tracer_solution[globalIdx]);
665 }
666 }
667 else {
668 tracer_model.setEnableSolTracers(tracer_index, false);
669
670 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
671 tracer_model.setSolTracerConcentration(tracer_index, elemIdx, 0.0);
672 }
673 }
674 }
675
676 if (inputThpres.active()) {
677 const_cast<Simulator&>(this->simulator_)
678 .problem().thresholdPressure()
679 .setFromRestart(restartValues.getExtra("THRESHPR"));
680 }
681
682 restartTimeStepSize_ = restartValues.getExtra("OPMEXTRA")[0];
683 if (restartTimeStepSize_ <= 0) {
684 restartTimeStepSize_ = std::numeric_limits<double>::max();
685 }
686
687 // Initialize the well model from restart values
688 this->simulator_.problem().wellModel()
689 .initFromRestartFile(restartValues);
690
691 if (!restartValues.aquifer.empty()) {
692 this->simulator_.problem().mutableAquiferModel()
693 .initFromRestart(restartValues.aquifer);
694 }
695 }
696 }
697
699 {
700 // Calculate initial in-place volumes.
701 // Does nothing if they have already been calculated,
702 // e.g. from restart data at T=0.
703 this->outputModule_->calc_initial_inplace(this->simulator_.gridView().comm());
704
705 if (this->collectOnIORank_.isIORank()) {
706 if (this->outputModule_->initialInplace().has_value()) {
707 this->inplace_ = this->outputModule_->initialInplace().value();
708 }
709 }
710 }
711
712 const OutputModule& outputModule() const
713 { return *outputModule_; }
714
715 OutputModule& mutableOutputModule() const
716 { return *outputModule_; }
717
718 Scalar restartTimeStepSize() const
719 { return restartTimeStepSize_; }
720
721 template <class Serializer>
722 void serializeOp(Serializer& serializer)
723 {
724 serializer(*outputModule_);
725 }
726
727private:
728 static bool enableEclOutput_()
729 {
730 static bool enable = Parameters::Get<Parameters::EnableEclOutput>();
731 return enable;
732 }
733
734 const EclipseState& eclState() const
735 { return simulator_.vanguard().eclState(); }
736
737 SummaryState& summaryState()
738 { return simulator_.vanguard().summaryState(); }
739
740 Action::State& actionState()
741 { return simulator_.vanguard().actionState(); }
742
743 UDQState& udqState()
744 { return simulator_.vanguard().udqState(); }
745
746 const Schedule& schedule() const
747 { return simulator_.vanguard().schedule(); }
748
749 void prepareLocalCellData(const bool isSubStep,
750 const int reportStepNum)
751 {
752 OPM_TIMEBLOCK(prepareLocalCellData);
753
754 if (this->outputModule_->localDataValid()) {
755 return;
756 }
757
758 const auto& gridView = simulator_.vanguard().gridView();
759 const bool log = this->collectOnIORank_.isIORank();
760
761 const int num_interior = detail::
763 this->outputModule_->
764 allocBuffers(num_interior, reportStepNum,
765 isSubStep && !Parameters::Get<Parameters::EnableWriteAllSolutions>(),
766 log, /*isRestart*/ false);
767
768 ElementContext elemCtx(simulator_);
769
771
772 {
773 OPM_TIMEBLOCK(prepareCellBasedData);
774
775 this->outputModule_->prepareDensityAccumulation();
776 this->outputModule_->setupExtractors(isSubStep, reportStepNum);
777 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
778 elemCtx.updatePrimaryStencil(elem);
779 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
780
781 this->outputModule_->processElement(elemCtx);
782 this->outputModule_->processElementBlockData(elemCtx);
783 }
784 this->outputModule_->clearExtractors();
785
786 this->outputModule_->accumulateDensityParallel();
787 }
788
789 {
790 OPM_TIMEBLOCK(prepareFluidInPlace);
791
792#ifdef _OPENMP
793#pragma omp parallel for
794#endif
795 for (int dofIdx = 0; dofIdx < num_interior; ++dofIdx) {
796 const auto& intQuants = *simulator_.model().cachedIntensiveQuantities(dofIdx, /*timeIdx=*/0);
797 const auto totVolume = simulator_.model().dofTotalVolume(dofIdx);
798
799 this->outputModule_->updateFluidInPlace(dofIdx, intQuants, totVolume);
800 }
801 }
802
803 this->outputModule_->validateLocalData();
804
805 OPM_END_PARALLEL_TRY_CATCH("EclWriter::prepareLocalCellData() failed: ",
806 this->simulator_.vanguard().grid().comm());
807 }
808
809 void captureLocalFluxData()
810 {
811 OPM_TIMEBLOCK(captureLocalData);
812
813 const auto& gridView = this->simulator_.vanguard().gridView();
814 const auto timeIdx = 0u;
815
816 auto elemCtx = ElementContext { this->simulator_ };
817
818 const auto elemMapper = ElementMapper { gridView, Dune::mcmgElementLayout() };
819 const auto activeIndex = [&elemMapper](const Element& e)
820 {
821 return elemMapper.index(e);
822 };
823
824 const auto cartesianIndex = [this](const int elemIndex)
825 {
826 return this->cartMapper_.cartesianIndex(elemIndex);
827 };
828
829 this->outputModule_->initializeFluxData();
830
832
833 for (const auto& elem : elements(gridView, Dune::Partitions::interiorBorder)) {
834 elemCtx.updateStencil(elem);
835 elemCtx.updateIntensiveQuantities(timeIdx);
836 elemCtx.updateExtensiveQuantities(timeIdx);
837
838 this->outputModule_->processFluxes(elemCtx, activeIndex, cartesianIndex);
839 }
840
841 OPM_END_PARALLEL_TRY_CATCH("EclWriter::captureLocalFluxData() failed: ",
842 this->simulator_.vanguard().grid().comm())
843
844 this->outputModule_->finalizeFluxData();
845 }
846
847 void writeWellspecReport(const SimulatorTimer& timer) const
848 {
849 const auto changedWells = this->schedule_
850 .changed_wells(timer.reportStepNum(), this->initialStep());
851
852 if (changedWells.empty()) {
853 return;
854 }
855
856 this->outputModule_->outputWellspecReport(changedWells,
857 timer.reportStepNum(),
858 timer.simulationTimeElapsed(),
859 timer.currentDateTime());
860 }
861
862 void writeWellflowReport(const SimulatorTimer& timer,
863 const int simStep,
864 const int wellsRequest) const
865 {
866 this->outputModule_->outputTimeStamp("WELLS",
867 timer.simulationTimeElapsed(),
868 timer.reportStepNum(),
869 timer.currentDateTime());
870
871 const auto wantConnData = wellsRequest > 1;
872
873 this->outputModule_->outputProdLog(simStep, wantConnData);
874 this->outputModule_->outputInjLog(simStep, wantConnData);
875 this->outputModule_->outputCumLog(simStep, wantConnData);
876 this->outputModule_->outputMSWLog(simStep);
877 }
878
879 int initialStep() const
880 {
881 const auto& initConfig = this->eclState().cfg().init();
882
883 return initConfig.restartRequested()
884 ? initConfig.getRestartStep()
885 : 0;
886 }
887
888 Simulator& simulator_;
889 std::unique_ptr<OutputModule> outputModule_;
890 Scalar restartTimeStepSize_;
891 int rank_ ;
892 Inplace inplace_;
893};
894
895} // namespace Opm
896
897#endif // OPM_ECL_WRITER_HPP
#define OPM_END_PARALLEL_TRY_CATCH(prefix, comm)
Catch exception and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:192
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:158
Declares the properties required by the black oil model.
int localIdxToGlobalIdx(unsigned localIdx) const
Definition: CollectDataOnIORank_impl.hpp:1102
InterRegFlowMap & globalInterRegFlows()
Definition: CollectDataOnIORank.hpp:114
bool isParallel() const
Definition: CollectDataOnIORank.hpp:129
bool isIORank() const
Definition: CollectDataOnIORank.hpp:126
const std::map< std::pair< std::string, int >, double > & globalBlockData() const
Definition: CollectDataOnIORank.hpp:90
void collect(const data::Solution &localCellData, const std::map< std::pair< std::string, int >, double > &localBlockData, std::map< std::pair< std::string, int >, double > &localExtraBlockData, const data::Wells &localWellData, const data::WellBlockAveragePressures &localWBPData, const data::GroupAndNetworkValues &localGroupAndNetworkData, const data::Aquifers &localAquiferData, const WellTestState &localWellTestState, const InterRegFlowMap &interRegFlows, const std::array< FlowsData< double >, 3 > &localFlowsn, const std::array< FlowsData< double >, 3 > &localFloresn)
Definition: CollectDataOnIORank_impl.hpp:969
const data::Solution & globalCellData() const
Definition: CollectDataOnIORank.hpp:93
Definition: EclGenericWriter.hpp:65
void doWriteOutput(const int reportStepNum, const std::optional< int > timeStepNum, 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< GetPropType< TypeTag, Properties::Scalar > > &thresholdPressure, GetPropType< TypeTag, Properties::Scalar > curTime, GetPropType< TypeTag, Properties::Scalar > nextStepSize, bool doublePrecision, bool isFlowsn, std::array< FlowsData< double >, 3 > &&flowsn, bool isFloresn, std::array< FlowsData< double >, 3 > &&floresn)
Definition: EclGenericWriter_impl.hpp:533
void evalSummary(int reportStepNum, GetPropType< TypeTag, Properties::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 std::optional< Inplace > &initialInPlace, const InterRegFlowMap &interRegFlows, SummaryState &summaryState, UDQState &udqState)
Definition: EclGenericWriter_impl.hpp:626
Collects necessary output values and pass it to opm-common's ECL output.
Definition: EclWriter.hpp:114
OutputModule & mutableOutputModule() const
Definition: EclWriter.hpp:715
const OutputModule & outputModule() const
Definition: EclWriter.hpp:712
void evalSummaryState(bool isSubStep)
collect and pass data and pass it to eclIO writer
Definition: EclWriter.hpp:202
static void registerParameters()
Definition: EclWriter.hpp:137
void serializeOp(Serializer &serializer)
Definition: EclWriter.hpp:722
void writeInitialFIPReport()
Writes the initial FIP report as configured in RPTSOL.
Definition: EclWriter.hpp:350
void beginRestart()
Definition: EclWriter.hpp:521
EclWriter(Simulator &simulator)
Definition: EclWriter.hpp:151
void writeOutput(data::Solution &&localCellData, bool isSubStep)
Definition: EclWriter.hpp:437
void writeReports(const SimulatorTimer &timer)
Definition: EclWriter.hpp:392
void endRestart()
Definition: EclWriter.hpp:698
Scalar restartTimeStepSize() const
Definition: EclWriter.hpp:718
~EclWriter()
Definition: EclWriter.hpp:191
const EquilGrid & globalGrid() const
Definition: EclWriter.hpp:194
virtual int reportStepNum() const
Current report step number. This might differ from currentStepNum in case of sub stepping.
Definition: SimulatorTimerInterface.hpp:109
Definition: SimulatorTimer.hpp:39
virtual boost::posix_time::ptime currentDateTime() const
Return the current time as a posix time object.
double simulationTimeElapsed() const override
Defines the common properties required by the porous medium multi-phase models.
Definition: ActionHandler.hpp:34
Definition: blackoilnewtonmethodparams.hpp:31
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition: parametersystem.hpp:185
std::size_t countLocalInteriorCellsGridView(const GridView &gridView)
Get the number of local interior cells in a grid view.
Definition: countGlobalCells.hpp:47
Definition: blackoilboundaryratevector.hh:39
data::Solution loadParallelRestartSolution(const EclipseIO *eclIO, const std::vector< RestartKey > &solutionKeys, Parallel::Communication comm, const int step)
void eclBroadcast(Parallel::Communication, T &)
RestartValue loadParallelRestart(const EclipseIO *eclIO, Action::State &actionState, SummaryState &summaryState, const std::vector< RestartKey > &solutionKeys, const std::vector< RestartKey > &extraKeys, Parallel::Communication comm)
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:233
Definition: EclWriter.hpp:72
static constexpr bool value
Definition: EclWriter.hpp:72
Definition: EclWriter.hpp:69
static constexpr bool value
Definition: EclWriter.hpp:69
Definition: EclWriter.hpp:79
static constexpr bool value
Definition: EclWriter.hpp:79
Definition: EclWriter.hpp:76
static constexpr bool value
Definition: EclWriter.hpp:76
SimulatorReportSingle success
Definition: SimulatorReport.hpp:123
unsigned int min_linear_iterations
Definition: SimulatorReport.hpp:52
unsigned int total_newton_iterations
Definition: SimulatorReport.hpp:50
unsigned int max_linear_iterations
Definition: SimulatorReport.hpp:53
unsigned int total_linear_iterations
Definition: SimulatorReport.hpp:51