My Project
Loading...
Searching...
No Matches
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
46#include <opm/simulators/flow/CollectDataOnIORank.hpp>
47#include <opm/simulators/flow/countGlobalCells.hpp>
50#include <opm/simulators/timestepping/SimulatorTimer.hpp>
51#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
52#include <opm/simulators/utils/ParallelRestart.hpp>
53#include <opm/simulators/utils/ParallelSerialization.hpp>
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
191 ~EclWriter()
192 {}
193
194 const EquilGrid& globalGrid() const
195 {
196 return simulator_.vanguard().equilGrid();
197 }
198
203 {
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()) {
248 OPM_BEGIN_PARALLEL_TRY_CATCH()
249
250 std::map<std::pair<std::string,int>,double> dummy;
251 this->collectOnIORank_.collect({},
252 outputModule_->getBlockData(),
253 dummy,
255 localWBP,
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
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;
284
285 {
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) {
298 }
299 if (this->sub_step_report_.total_newton_iterations != 0) {
300 miscSummaryData["NEWTON"] = this->sub_step_report_.total_newton_iterations;
301 }
302 if (this->sub_step_report_.total_linear_iterations != 0) {
303 miscSummaryData["MLINEARS"] = this->sub_step_report_.total_linear_iterations;
304 }
305 if (this->sub_step_report_.total_newton_iterations != 0) {
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 }
311 if (this->sub_step_report_.max_linear_iterations != 0) {
312 miscSummaryData["NLINSMAX"] = this->sub_step_report_.max_linear_iterations;
313 }
314 if (this->simulation_report_.success.total_linear_iterations != 0) {
315 miscSummaryData["MSUMLINS"] = this->simulation_report_.success.total_linear_iterations;
316 }
317 if (this->simulation_report_.success.total_newton_iterations != 0) {
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()
325 ? this->collectOnIORank_.globalBlockData()
326 : this->outputModule_->getBlockData();
327
328 const auto& interRegFlows = this->collectOnIORank_.isParallel()
329 ? this->collectOnIORank_.globalInterRegFlows()
330 : this->outputModule_->getInterRegFlows();
331
332 this->evalSummary(reportStepNum,
333 curTime,
335 localWBP,
338 blockData,
341 inplace,
342 this->outputModule_->initialInplace(),
344 this->summaryState(),
345 this->udqState());
346 }
347 }
348
351 {
352 const auto& gridView = simulator_.vanguard().gridView();
353 const int num_interior = detail::
354 countLocalInteriorCellsGridView(gridView);
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 {
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
390 void writeReports(const SimulatorTimer& timer)
391 {
392 if (! this->collectOnIORank_.isIORank()) {
393 return;
394 }
395
396 // SimulatorTimer::reportStepNum() is the simulator's zero-based
397 // "episode index". This is generally the index value needed to
398 // look up objects in the Schedule container. That said, function
399 // writeReports() is invoked at the *beginning* of a report
400 // step/episode which means we typically need the objects from the
401 // *previous* report step/episode. We therefore need special case
402 // handling for reportStepNum() == 0 in base runs and
403 // reportStepNum() <= restart step in restarted runs.
404 const auto firstStep = this->initialStep();
405 const auto simStep =
406 std::max(timer.reportStepNum() - 1, firstStep);
407
408 const auto& rpt = this->schedule_[simStep].rpt_config();
409
410 if (rpt.contains("WELSPECS") && (rpt.at("WELSPECS") > 0)) {
411 // Requesting a well specification report is valid at all times,
412 // including reportStepNum() == initialStep().
413 this->writeWellspecReport(timer);
414 }
415
416 if (timer.reportStepNum() == firstStep) {
417 // No dynamic flows at the beginning of the initialStep().
418 return;
419 }
420
421 if (rpt.contains("WELLS") && rpt.at("WELLS") > 0) {
422 this->writeWellflowReport(timer, simStep, rpt.at("WELLS"));
423 }
424
425 this->outputModule_->outputFipAndResvLog(this->inplace_,
426 timer.reportStepNum(),
427 timer.simulationTimeElapsed(),
428 timer.currentDateTime(),
429 /* isSubstep = */ false,
430 simulator_.gridView().comm());
431
432 OpmLog::note(""); // Blank line after all reports.
433 }
434
435 void writeOutput(data::Solution&& localCellData, bool isSubStep)
436 {
437 OPM_TIMEBLOCK(writeOutput);
438
439 const int reportStepNum = simulator_.episodeIndex() + 1;
440 this->prepareLocalCellData(isSubStep, reportStepNum);
441 this->outputModule_->outputErrorLog(simulator_.gridView().comm());
442
443 // output using eclWriter if enabled
444 auto localWellData = simulator_.problem().wellModel().wellData();
445 auto localGroupAndNetworkData = simulator_.problem().wellModel()
446 .groupAndNetworkData(reportStepNum);
447
448 auto localAquiferData = simulator_.problem().aquiferModel().aquiferData();
449 auto localWellTestState = simulator_.problem().wellModel().wellTestState();
450
451 const bool isFlowsn = this->outputModule_->getFlows().hasFlowsn();
452 auto flowsn = this->outputModule_->getFlows().getFlowsn();
453
454 const bool isFloresn = this->outputModule_->getFlows().hasFloresn();
455 auto floresn = this->outputModule_->getFlows().getFloresn();
456
457 if (! isSubStep || Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
458
459 if (localCellData.empty()) {
460 this->outputModule_->assignToSolution(localCellData);
461 }
462
463 // Add cell data to perforations for RFT output
464 this->outputModule_->addRftDataToWells(localWellData,
465 reportStepNum,
466 simulator_.gridView().comm());
467 }
468
469 if (this->collectOnIORank_.isParallel() ||
470 this->collectOnIORank_.doesNeedReordering())
471 {
472 // Note: We don't need WBP (well-block averaged pressures) or
473 // inter-region flow rate values in order to create restart file
474 // output. There's consequently no need to collect those
475 // properties on the I/O rank.
476
477 this->collectOnIORank_.collect(localCellData,
478 this->outputModule_->getBlockData(),
479 this->outputModule_->getExtraBlockData(),
481 /* wbpData = */ {},
485 /* interRegFlows = */ {},
486 flowsn,
487 floresn);
488 if (this->collectOnIORank_.isIORank()) {
489 this->outputModule_->assignGlobalFieldsToSolution(this->collectOnIORank_.globalCellData());
490 }
491 } else {
492 this->outputModule_->assignGlobalFieldsToSolution(localCellData);
493 }
494
495 if (this->collectOnIORank_.isIORank()) {
496 const Scalar curTime = simulator_.time() + simulator_.timeStepSize();
497 const Scalar nextStepSize = simulator_.problem().nextTimeStepSize();
498 std::optional<int> timeStepIdx;
499 if (Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
500 timeStepIdx = simulator_.timeStepIndex();
501 }
502 this->doWriteOutput(reportStepNum, timeStepIdx, isSubStep,
503 std::move(localCellData),
504 std::move(localWellData),
505 std::move(localGroupAndNetworkData),
506 std::move(localAquiferData),
507 std::move(localWellTestState),
508 this->actionState(),
509 this->udqState(),
510 this->summaryState(),
511 this->simulator_.problem().thresholdPressure().getRestartVector(),
513 Parameters::Get<Parameters::EclOutputDoublePrecision>(),
514 isFlowsn, std::move(flowsn),
515 isFloresn, std::move(floresn));
516 }
517 }
518
519 void beginRestart()
520 {
521 const auto enablePCHysteresis = simulator_.problem().materialLawManager()->enablePCHysteresis();
522 const auto enableNonWettingHysteresis = simulator_.problem().materialLawManager()->enableNonWettingHysteresis();
523 const auto enableWettingHysteresis = simulator_.problem().materialLawManager()->enableWettingHysteresis();
524 const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
525 const auto gasActive = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
526 const auto waterActive = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
527 const auto enableSwatinit = simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT");
528
529 std::vector<RestartKey> solutionKeys {
530 {"PRESSURE", UnitSystem::measure::pressure},
531 {"SWAT", UnitSystem::measure::identity, waterActive},
532 {"SGAS", UnitSystem::measure::identity, gasActive},
533 {"TEMP", UnitSystem::measure::temperature, enableEnergy},
534 {"SSOLVENT", UnitSystem::measure::identity, enableSolvent},
535
536 {"RS", UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()},
537 {"RV", UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()},
538 {"RVW", UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedWater()},
539 {"RSW", UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGasInWater()},
540
541 {"SGMAX", UnitSystem::measure::identity, enableNonWettingHysteresis && oilActive && gasActive},
542 {"SHMAX", UnitSystem::measure::identity, enableWettingHysteresis && oilActive && gasActive},
543
544 {"SOMAX", UnitSystem::measure::identity,
545 (enableNonWettingHysteresis && oilActive && waterActive)
546 || simulator_.problem().vapparsActive(simulator_.episodeIndex())},
547
548 {"SOMIN", UnitSystem::measure::identity, enablePCHysteresis && oilActive && gasActive},
549 {"SWHY1", UnitSystem::measure::identity, enablePCHysteresis && oilActive && waterActive},
550 {"SWMAX", UnitSystem::measure::identity, enableWettingHysteresis && oilActive && waterActive},
551
552 {"PPCW", UnitSystem::measure::pressure, enableSwatinit},
553 };
554
555 {
556 const auto& tracers = simulator_.vanguard().eclState().tracer();
557
558 for (const auto& tracer : tracers) {
559 const auto enableSolTracer =
560 ((tracer.phase == Phase::GAS) && FluidSystem::enableDissolvedGas()) ||
561 ((tracer.phase == Phase::OIL) && FluidSystem::enableVaporizedOil());
562
563 solutionKeys.emplace_back(tracer.fname(), UnitSystem::measure::identity, true);
564 solutionKeys.emplace_back(tracer.sname(), UnitSystem::measure::identity, enableSolTracer);
565 }
566 }
567
568 const auto& inputThpres = eclState().getSimulationConfig().getThresholdPressure();
569 const std::vector<RestartKey> extraKeys {
570 {"OPMEXTRA", UnitSystem::measure::identity, false},
571 {"THRESHPR", UnitSystem::measure::pressure, inputThpres.active()},
572 };
573
574 const auto& gridView = this->simulator_.vanguard().gridView();
575 const auto numElements = gridView.size(/*codim=*/0);
576
577 // Try to load restart step 0 to calculate initial FIP
578 {
579 this->outputModule_->allocBuffers(numElements,
580 0,
581 /*isSubStep = */false,
582 /*log = */ false,
583 /*isRestart = */true);
584
585 const auto restartSolution =
586 loadParallelRestartSolution(this->eclIO_.get(),
587 solutionKeys, gridView.comm(), 0);
588
589 if (!restartSolution.empty()) {
590 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
591 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
592 this->outputModule_->setRestart(restartSolution, elemIdx, globalIdx);
593 }
594
595 this->simulator_.problem().readSolutionFromOutputModule(0, true);
596 ElementContext elemCtx(this->simulator_);
597 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
598 elemCtx.updatePrimaryStencil(elem);
599 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
600
601 this->outputModule_->updateFluidInPlace(elemCtx);
602 }
603
604 this->outputModule_->calc_initial_inplace(this->simulator_.gridView().comm());
605 }
606 }
607
608 {
609 // The episodeIndex is rewound one step back before calling
610 // beginRestart() and cannot be used here. We just ask the
611 // initconfig directly to be sure that we use the correct index.
612 const auto restartStepIdx = this->simulator_.vanguard()
613 .eclState().getInitConfig().getRestartStep();
614
615 this->outputModule_->allocBuffers(numElements,
617 /*isSubStep = */false,
618 /*log = */ false,
619 /*isRestart = */true);
620 }
621
622 {
623 const auto restartValues =
624 loadParallelRestart(this->eclIO_.get(),
625 this->actionState(),
626 this->summaryState(),
627 solutionKeys, extraKeys, gridView.comm());
628
629 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
630 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
631 this->outputModule_->setRestart(restartValues.solution, elemIdx, globalIdx);
632 }
633
634 auto& tracer_model = simulator_.problem().tracerModel();
635 for (int tracer_index = 0; tracer_index < tracer_model.numTracers(); ++tracer_index) {
636 // Free tracers
637 {
638 const auto& free_tracer_name = tracer_model.fname(tracer_index);
639 const auto& free_tracer_solution = restartValues.solution
641
642 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
643 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
644 tracer_model.setFreeTracerConcentration
645 (tracer_index, elemIdx, free_tracer_solution[globalIdx]);
646 }
647 }
648
649 // Solution tracer (only if DISGAS/VAPOIL are active for gas/oil tracers)
650 if ((tracer_model.phase(tracer_index) == Phase::GAS && FluidSystem::enableDissolvedGas()) ||
651 (tracer_model.phase(tracer_index) == Phase::OIL && FluidSystem::enableVaporizedOil()))
652 {
653 tracer_model.setEnableSolTracers(tracer_index, true);
654
655 const auto& sol_tracer_name = tracer_model.sname(tracer_index);
656 const auto& sol_tracer_solution = restartValues.solution
657 .template data<double>(sol_tracer_name);
658
659 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
660 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
661 tracer_model.setSolTracerConcentration
662 (tracer_index, elemIdx, sol_tracer_solution[globalIdx]);
663 }
664 }
665 else {
666 tracer_model.setEnableSolTracers(tracer_index, false);
667
668 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
669 tracer_model.setSolTracerConcentration(tracer_index, elemIdx, 0.0);
670 }
671 }
672 }
673
674 if (inputThpres.active()) {
675 const_cast<Simulator&>(this->simulator_)
676 .problem().thresholdPressure()
677 .setFromRestart(restartValues.getExtra("THRESHPR"));
678 }
679
680 restartTimeStepSize_ = restartValues.getExtra("OPMEXTRA")[0];
681 if (restartTimeStepSize_ <= 0) {
682 restartTimeStepSize_ = std::numeric_limits<double>::max();
683 }
684
685 // Initialize the well model from restart values
686 this->simulator_.problem().wellModel()
687 .initFromRestartFile(restartValues);
688
689 if (!restartValues.aquifer.empty()) {
690 this->simulator_.problem().mutableAquiferModel()
691 .initFromRestart(restartValues.aquifer);
692 }
693 }
694 }
695
696 void endRestart()
697 {
698 // Calculate initial in-place volumes.
699 // Does nothing if they have already been calculated,
700 // e.g. from restart data at T=0.
701 this->outputModule_->calc_initial_inplace(this->simulator_.gridView().comm());
702
703 if (this->collectOnIORank_.isIORank()) {
704 if (this->outputModule_->initialInplace().has_value()) {
705 this->inplace_ = this->outputModule_->initialInplace().value();
706 }
707 }
708 }
709
710 const OutputModule& outputModule() const
711 { return *outputModule_; }
712
713 OutputModule& mutableOutputModule() const
714 { return *outputModule_; }
715
716 Scalar restartTimeStepSize() const
717 { return restartTimeStepSize_; }
718
719 template <class Serializer>
720 void serializeOp(Serializer& serializer)
721 {
722 serializer(*outputModule_);
723 }
724
725private:
726 static bool enableEclOutput_()
727 {
728 static bool enable = Parameters::Get<Parameters::EnableEclOutput>();
729 return enable;
730 }
731
732 const EclipseState& eclState() const
733 { return simulator_.vanguard().eclState(); }
734
735 SummaryState& summaryState()
736 { return simulator_.vanguard().summaryState(); }
737
738 Action::State& actionState()
739 { return simulator_.vanguard().actionState(); }
740
741 UDQState& udqState()
742 { return simulator_.vanguard().udqState(); }
743
744 const Schedule& schedule() const
745 { return simulator_.vanguard().schedule(); }
746
747 void prepareLocalCellData(const bool isSubStep,
748 const int reportStepNum)
749 {
750 OPM_TIMEBLOCK(prepareLocalCellData);
751
752 if (this->outputModule_->localDataValid()) {
753 return;
754 }
755
756 const auto& gridView = simulator_.vanguard().gridView();
757 const bool log = this->collectOnIORank_.isIORank();
758
759 const int num_interior = detail::
760 countLocalInteriorCellsGridView(gridView);
761 this->outputModule_->
762 allocBuffers(num_interior, reportStepNum,
763 isSubStep && !Parameters::Get<Parameters::EnableWriteAllSolutions>(),
764 log, /*isRestart*/ false);
765
766 ElementContext elemCtx(simulator_);
767
768 OPM_BEGIN_PARALLEL_TRY_CATCH();
769
770 {
772
773 this->outputModule_->prepareDensityAccumulation();
774 this->outputModule_->setupExtractors(isSubStep, reportStepNum);
775 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
776 elemCtx.updatePrimaryStencil(elem);
777 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
778
779 this->outputModule_->processElement(elemCtx);
780 this->outputModule_->processElementBlockData(elemCtx);
781 }
782 this->outputModule_->clearExtractors();
783
784 this->outputModule_->accumulateDensityParallel();
785 }
786
787 {
789
790#ifdef _OPENMP
791#pragma omp parallel for
792#endif
793 for (int dofIdx = 0; dofIdx < num_interior; ++dofIdx) {
794 const auto& intQuants = *simulator_.model().cachedIntensiveQuantities(dofIdx, /*timeIdx=*/0);
795 const auto totVolume = simulator_.model().dofTotalVolume(dofIdx);
796
797 this->outputModule_->updateFluidInPlace(dofIdx, intQuants, totVolume);
798 }
799 }
800
801 this->outputModule_->validateLocalData();
802
803 OPM_END_PARALLEL_TRY_CATCH("EclWriter::prepareLocalCellData() failed: ",
804 this->simulator_.vanguard().grid().comm());
805 }
806
807 void captureLocalFluxData()
808 {
810
811 const auto& gridView = this->simulator_.vanguard().gridView();
812 const auto timeIdx = 0u;
813
814 auto elemCtx = ElementContext { this->simulator_ };
815
816 const auto elemMapper = ElementMapper { gridView, Dune::mcmgElementLayout() };
817 const auto activeIndex = [&elemMapper](const Element& e)
818 {
819 return elemMapper.index(e);
820 };
821
822 const auto cartesianIndex = [this](const int elemIndex)
823 {
824 return this->cartMapper_.cartesianIndex(elemIndex);
825 };
826
827 this->outputModule_->initializeFluxData();
828
829 OPM_BEGIN_PARALLEL_TRY_CATCH();
830
831 for (const auto& elem : elements(gridView, Dune::Partitions::interiorBorder)) {
832 elemCtx.updateStencil(elem);
833 elemCtx.updateIntensiveQuantities(timeIdx);
834 elemCtx.updateExtensiveQuantities(timeIdx);
835
836 this->outputModule_->processFluxes(elemCtx, activeIndex, cartesianIndex);
837 }
838
839 OPM_END_PARALLEL_TRY_CATCH("EclWriter::captureLocalFluxData() failed: ",
840 this->simulator_.vanguard().grid().comm())
841
842 this->outputModule_->finalizeFluxData();
843 }
844
845 void writeWellspecReport(const SimulatorTimer& timer) const
846 {
847 const auto changedWells = this->schedule_
848 .changed_wells(timer.reportStepNum(), this->initialStep());
849
850 if (changedWells.empty()) {
851 return;
852 }
853
854 this->outputModule_->outputWellspecReport(changedWells,
855 timer.reportStepNum(),
856 timer.simulationTimeElapsed(),
857 timer.currentDateTime());
858 }
859
860 void writeWellflowReport(const SimulatorTimer& timer,
861 const int simStep,
862 const int wellsRequest) const
863 {
864 this->outputModule_->outputTimeStamp("WELLS",
865 timer.simulationTimeElapsed(),
866 timer.reportStepNum(),
867 timer.currentDateTime());
868
869 const auto wantConnData = wellsRequest > 1;
870
871 this->outputModule_->outputProdLog(simStep, wantConnData);
872 this->outputModule_->outputInjLog(simStep, wantConnData);
873 this->outputModule_->outputCumLog(simStep, wantConnData);
874 this->outputModule_->outputMSWLog(simStep);
875 }
876
877 int initialStep() const
878 {
879 const auto& initConfig = this->eclState().cfg().init();
880
881 return initConfig.restartRequested()
882 ? initConfig.getRestartStep()
883 : 0;
884 }
885
886 Simulator& simulator_;
887 std::unique_ptr<OutputModule> outputModule_;
888 Scalar restartTimeStepSize_;
889 int rank_ ;
890 Inplace inplace_;
891};
892
893} // namespace Opm
894
895#endif // OPM_ECL_WRITER_HPP
Collects necessary output values and pass it to opm-common's ECL output.
Helper class for grid instantiation of ECL file-format using problems.
Declares the properties required by the black oil model.
Definition EclGenericWriter.hpp:65
Collects necessary output values and pass it to opm-common's ECL output.
Definition EclWriter.hpp:114
void evalSummaryState(bool isSubStep)
collect and pass data and pass it to eclIO writer
Definition EclWriter.hpp:202
void writeInitialFIPReport()
Writes the initial FIP report as configured in RPTSOL.
Definition EclWriter.hpp:350
void compress()
Form CSR adjacency matrix representation of input graph from connections established in previous call...
Definition InterRegFlows.cpp:165
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
double simulationTimeElapsed() const override
Time elapsed since the start of the simulation until the beginning of the current time step [s].
Definition SimulatorTimer.cpp:122
virtual boost::posix_time::ptime currentDateTime() const
Return the current time as a posix time object.
Definition SimulatorTimerInterface.cpp:28
Defines the common properties required by the porous medium multi-phase models.
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
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:235
Definition EclWriter.hpp:69
Definition EclWriter.hpp:79
Definition EclWriter.hpp:76