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();
384
385 outputModule_->outputFipAndResvLog(inplace_, 0, 0.0, start_time,
386 false, simulator_.gridView().comm());
387 }
388 }
389 }
390
391 void writeReports(const SimulatorTimer& timer)
392 {
393 auto rstep = timer.reportStepNum();
394
395 if ((rstep > 0) && (this->collectOnIORank_.isIORank())) {
396 const auto& rpt = this->schedule_[rstep-1].rpt_config.get();
397 if (rpt.contains("WELLS") && rpt.at("WELLS") > 0) {
398 outputModule_->outputTimeStamp("WELLS",
399 timer.simulationTimeElapsed(),
400 rstep,
401 timer.currentDateTime());
402 outputModule_->outputProdLog(rstep - 1, rpt.at("WELLS") > 1);
403 outputModule_->outputInjLog(rstep - 1, rpt.at("WELLS") > 1);
404 outputModule_->outputCumLog(rstep - 1, rpt.at("WELLS") > 1);
405 outputModule_->outputMSWLog(rstep - 1);
406 }
407
408 outputModule_->outputFipAndResvLog(inplace_,
409 rstep,
410 timer.simulationTimeElapsed(),
411 timer.currentDateTime(),
412 false,
413 simulator_.gridView().comm());
414
415
416 OpmLog::note(""); // Blank line after all reports.
417 }
418 }
419
420 void writeOutput(data::Solution&& localCellData, bool isSubStep)
421 {
422 OPM_TIMEBLOCK(writeOutput);
423
424 const int reportStepNum = simulator_.episodeIndex() + 1;
425 this->prepareLocalCellData(isSubStep, reportStepNum);
426 this->outputModule_->outputErrorLog(simulator_.gridView().comm());
427
428 // output using eclWriter if enabled
429 auto localWellData = simulator_.problem().wellModel().wellData();
430 auto localGroupAndNetworkData = simulator_.problem().wellModel()
431 .groupAndNetworkData(reportStepNum);
432
433 auto localAquiferData = simulator_.problem().aquiferModel().aquiferData();
434 auto localWellTestState = simulator_.problem().wellModel().wellTestState();
435
436 const bool isFlowsn = this->outputModule_->getFlows().hasFlowsn();
437 auto flowsn = this->outputModule_->getFlows().getFlowsn();
438
439 const bool isFloresn = this->outputModule_->getFlows().hasFloresn();
440 auto floresn = this->outputModule_->getFlows().getFloresn();
441
442 if (! isSubStep || Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
443
444 if (localCellData.empty()) {
445 this->outputModule_->assignToSolution(localCellData);
446 }
447
448 // Add cell data to perforations for RFT output
449 this->outputModule_->addRftDataToWells(localWellData,
450 reportStepNum,
451 simulator_.gridView().comm());
452 }
453
454 if (this->collectOnIORank_.isParallel() ||
455 this->collectOnIORank_.doesNeedReordering())
456 {
457 // Note: We don't need WBP (well-block averaged pressures) or
458 // inter-region flow rate values in order to create restart file
459 // output. There's consequently no need to collect those
460 // properties on the I/O rank.
461
462 this->collectOnIORank_.collect(localCellData,
463 this->outputModule_->getBlockData(),
464 this->outputModule_->getExtraBlockData(),
466 /* wbpData = */ {},
470 /* interRegFlows = */ {},
471 flowsn,
472 floresn);
473 if (this->collectOnIORank_.isIORank()) {
474 this->outputModule_->assignGlobalFieldsToSolution(this->collectOnIORank_.globalCellData());
475 }
476 } else {
477 this->outputModule_->assignGlobalFieldsToSolution(localCellData);
478 }
479
480 if (this->collectOnIORank_.isIORank()) {
481 const Scalar curTime = simulator_.time() + simulator_.timeStepSize();
482 const Scalar nextStepSize = simulator_.problem().nextTimeStepSize();
483 std::optional<int> timeStepIdx;
484 if (Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
485 timeStepIdx = simulator_.timeStepIndex();
486 }
487 this->doWriteOutput(reportStepNum, timeStepIdx, isSubStep,
488 std::move(localCellData),
489 std::move(localWellData),
490 std::move(localGroupAndNetworkData),
491 std::move(localAquiferData),
492 std::move(localWellTestState),
493 this->actionState(),
494 this->udqState(),
495 this->summaryState(),
496 this->simulator_.problem().thresholdPressure().getRestartVector(),
498 Parameters::Get<Parameters::EclOutputDoublePrecision>(),
499 isFlowsn, std::move(flowsn),
500 isFloresn, std::move(floresn));
501 }
502 }
503
504 void beginRestart()
505 {
506 const auto enablePCHysteresis = simulator_.problem().materialLawManager()->enablePCHysteresis();
507 const auto enableNonWettingHysteresis = simulator_.problem().materialLawManager()->enableNonWettingHysteresis();
508 const auto enableWettingHysteresis = simulator_.problem().materialLawManager()->enableWettingHysteresis();
509 const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
510 const auto gasActive = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
511 const auto waterActive = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
512 const auto enableSwatinit = simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT");
513
514 std::vector<RestartKey> solutionKeys {
515 {"PRESSURE", UnitSystem::measure::pressure},
516 {"SWAT", UnitSystem::measure::identity, waterActive},
517 {"SGAS", UnitSystem::measure::identity, gasActive},
518 {"TEMP", UnitSystem::measure::temperature, enableEnergy},
519 {"SSOLVENT", UnitSystem::measure::identity, enableSolvent},
520
521 {"RS", UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()},
522 {"RV", UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()},
523 {"RVW", UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedWater()},
524 {"RSW", UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGasInWater()},
525
526 {"SGMAX", UnitSystem::measure::identity, enableNonWettingHysteresis && oilActive && gasActive},
527 {"SHMAX", UnitSystem::measure::identity, enableWettingHysteresis && oilActive && gasActive},
528
529 {"SOMAX", UnitSystem::measure::identity,
530 (enableNonWettingHysteresis && oilActive && waterActive)
531 || simulator_.problem().vapparsActive(simulator_.episodeIndex())},
532
533 {"SOMIN", UnitSystem::measure::identity, enablePCHysteresis && oilActive && gasActive},
534 {"SWHY1", UnitSystem::measure::identity, enablePCHysteresis && oilActive && waterActive},
535 {"SWMAX", UnitSystem::measure::identity, enableWettingHysteresis && oilActive && waterActive},
536
537 {"PPCW", UnitSystem::measure::pressure, enableSwatinit},
538 };
539
540 {
541 const auto& tracers = simulator_.vanguard().eclState().tracer();
542
543 for (const auto& tracer : tracers) {
544 const auto enableSolTracer =
545 ((tracer.phase == Phase::GAS) && FluidSystem::enableDissolvedGas()) ||
546 ((tracer.phase == Phase::OIL) && FluidSystem::enableVaporizedOil());
547
548 solutionKeys.emplace_back(tracer.fname(), UnitSystem::measure::identity, true);
549 solutionKeys.emplace_back(tracer.sname(), UnitSystem::measure::identity, enableSolTracer);
550 }
551 }
552
553 const auto& inputThpres = eclState().getSimulationConfig().getThresholdPressure();
554 const std::vector<RestartKey> extraKeys {
555 {"OPMEXTRA", UnitSystem::measure::identity, false},
556 {"THRESHPR", UnitSystem::measure::pressure, inputThpres.active()},
557 };
558
559 const auto& gridView = this->simulator_.vanguard().gridView();
560 const auto numElements = gridView.size(/*codim=*/0);
561
562 // Try to load restart step 0 to calculate initial FIP
563 {
564 this->outputModule_->allocBuffers(numElements,
565 0,
566 /*isSubStep = */false,
567 /*log = */ false,
568 /*isRestart = */true);
569
570 const auto restartSolution =
571 loadParallelRestartSolution(this->eclIO_.get(),
572 solutionKeys, gridView.comm(), 0);
573
574 if (!restartSolution.empty()) {
575 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
576 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
577 this->outputModule_->setRestart(restartSolution, elemIdx, globalIdx);
578 }
579
580 this->simulator_.problem().readSolutionFromOutputModule(0, true);
581 ElementContext elemCtx(this->simulator_);
582 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
583 elemCtx.updatePrimaryStencil(elem);
584 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
585
586 this->outputModule_->updateFluidInPlace(elemCtx);
587 }
588
589 this->outputModule_->calc_initial_inplace(this->simulator_.gridView().comm());
590 }
591 }
592
593 {
594 // The episodeIndex is rewound one step back before calling
595 // beginRestart() and cannot be used here. We just ask the
596 // initconfig directly to be sure that we use the correct index.
597 const auto restartStepIdx = this->simulator_.vanguard()
598 .eclState().getInitConfig().getRestartStep();
599
600 this->outputModule_->allocBuffers(numElements,
602 /*isSubStep = */false,
603 /*log = */ false,
604 /*isRestart = */true);
605 }
606
607 {
608 const auto restartValues =
609 loadParallelRestart(this->eclIO_.get(),
610 this->actionState(),
611 this->summaryState(),
612 solutionKeys, extraKeys, gridView.comm());
613
614 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
615 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
616 this->outputModule_->setRestart(restartValues.solution, elemIdx, globalIdx);
617 }
618
619 auto& tracer_model = simulator_.problem().tracerModel();
620 for (int tracer_index = 0; tracer_index < tracer_model.numTracers(); ++tracer_index) {
621 // Free tracers
622 {
623 const auto& free_tracer_name = tracer_model.fname(tracer_index);
624 const auto& free_tracer_solution = restartValues.solution
626
627 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
628 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
629 tracer_model.setFreeTracerConcentration
630 (tracer_index, elemIdx, free_tracer_solution[globalIdx]);
631 }
632 }
633
634 // Solution tracer (only if DISGAS/VAPOIL are active for gas/oil tracers)
635 if ((tracer_model.phase(tracer_index) == Phase::GAS && FluidSystem::enableDissolvedGas()) ||
636 (tracer_model.phase(tracer_index) == Phase::OIL && FluidSystem::enableVaporizedOil()))
637 {
638 tracer_model.setEnableSolTracers(tracer_index, true);
639
640 const auto& sol_tracer_name = tracer_model.sname(tracer_index);
641 const auto& sol_tracer_solution = restartValues.solution
642 .template data<double>(sol_tracer_name);
643
644 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
645 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
646 tracer_model.setSolTracerConcentration
647 (tracer_index, elemIdx, sol_tracer_solution[globalIdx]);
648 }
649 }
650 else {
651 tracer_model.setEnableSolTracers(tracer_index, false);
652
653 for (auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
654 tracer_model.setSolTracerConcentration(tracer_index, elemIdx, 0.0);
655 }
656 }
657 }
658
659 if (inputThpres.active()) {
660 const_cast<Simulator&>(this->simulator_)
661 .problem().thresholdPressure()
662 .setFromRestart(restartValues.getExtra("THRESHPR"));
663 }
664
665 restartTimeStepSize_ = restartValues.getExtra("OPMEXTRA")[0];
666 if (restartTimeStepSize_ <= 0) {
667 restartTimeStepSize_ = std::numeric_limits<double>::max();
668 }
669
670 // Initialize the well model from restart values
671 this->simulator_.problem().wellModel()
672 .initFromRestartFile(restartValues);
673
674 if (!restartValues.aquifer.empty()) {
675 this->simulator_.problem().mutableAquiferModel()
676 .initFromRestart(restartValues.aquifer);
677 }
678 }
679 }
680
681 void endRestart()
682 {
683 // Calculate initial in-place volumes.
684 // Does nothing if they have already been calculated,
685 // e.g. from restart data at T=0.
686 this->outputModule_->calc_initial_inplace(this->simulator_.gridView().comm());
687
688 if (this->collectOnIORank_.isIORank()) {
689 this->inplace_ = this->outputModule_->initialInplace();
690 }
691 }
692
693 const OutputModule& outputModule() const
694 { return *outputModule_; }
695
696 OutputModule& mutableOutputModule() const
697 { return *outputModule_; }
698
699 Scalar restartTimeStepSize() const
700 { return restartTimeStepSize_; }
701
702 template <class Serializer>
703 void serializeOp(Serializer& serializer)
704 {
705 serializer(*outputModule_);
706 }
707
708private:
709 static bool enableEclOutput_()
710 {
711 static bool enable = Parameters::Get<Parameters::EnableEclOutput>();
712 return enable;
713 }
714
715 const EclipseState& eclState() const
716 { return simulator_.vanguard().eclState(); }
717
718 SummaryState& summaryState()
719 { return simulator_.vanguard().summaryState(); }
720
721 Action::State& actionState()
722 { return simulator_.vanguard().actionState(); }
723
724 UDQState& udqState()
725 { return simulator_.vanguard().udqState(); }
726
727 const Schedule& schedule() const
728 { return simulator_.vanguard().schedule(); }
729
730 void prepareLocalCellData(const bool isSubStep,
731 const int reportStepNum)
732 {
733 OPM_TIMEBLOCK(prepareLocalCellData);
734
735 if (this->outputModule_->localDataValid()) {
736 return;
737 }
738
739 const auto& gridView = simulator_.vanguard().gridView();
740 const bool log = this->collectOnIORank_.isIORank();
741
742 const int num_interior = detail::
743 countLocalInteriorCellsGridView(gridView);
744 this->outputModule_->
745 allocBuffers(num_interior, reportStepNum,
746 isSubStep && !Parameters::Get<Parameters::EnableWriteAllSolutions>(),
747 log, /*isRestart*/ false);
748
749 ElementContext elemCtx(simulator_);
750
751 OPM_BEGIN_PARALLEL_TRY_CATCH();
752
753 {
755
756 this->outputModule_->prepareDensityAccumulation();
757 this->outputModule_->setupExtractors(isSubStep, reportStepNum);
758 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
759 elemCtx.updatePrimaryStencil(elem);
760 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
761
762 this->outputModule_->processElement(elemCtx);
763 this->outputModule_->processElementBlockData(elemCtx);
764 }
765 this->outputModule_->clearExtractors();
766
767 this->outputModule_->accumulateDensityParallel();
768 }
769
770 {
772
773#ifdef _OPENMP
774#pragma omp parallel for
775#endif
776 for (int dofIdx = 0; dofIdx < num_interior; ++dofIdx) {
777 const auto& intQuants = *simulator_.model().cachedIntensiveQuantities(dofIdx, /*timeIdx=*/0);
778 const auto totVolume = simulator_.model().dofTotalVolume(dofIdx);
779
780 this->outputModule_->updateFluidInPlace(dofIdx, intQuants, totVolume);
781 }
782 }
783
784 this->outputModule_->validateLocalData();
785
786 OPM_END_PARALLEL_TRY_CATCH("EclWriter::prepareLocalCellData() failed: ",
787 this->simulator_.vanguard().grid().comm());
788 }
789
790 void captureLocalFluxData()
791 {
793
794 const auto& gridView = this->simulator_.vanguard().gridView();
795 const auto timeIdx = 0u;
796
797 auto elemCtx = ElementContext { this->simulator_ };
798
799 const auto elemMapper = ElementMapper { gridView, Dune::mcmgElementLayout() };
800 const auto activeIndex = [&elemMapper](const Element& e)
801 {
802 return elemMapper.index(e);
803 };
804
805 const auto cartesianIndex = [this](const int elemIndex)
806 {
807 return this->cartMapper_.cartesianIndex(elemIndex);
808 };
809
810 this->outputModule_->initializeFluxData();
811
812 OPM_BEGIN_PARALLEL_TRY_CATCH();
813
814 for (const auto& elem : elements(gridView, Dune::Partitions::interiorBorder)) {
815 elemCtx.updateStencil(elem);
816 elemCtx.updateIntensiveQuantities(timeIdx);
817 elemCtx.updateExtensiveQuantities(timeIdx);
818
819 this->outputModule_->processFluxes(elemCtx, activeIndex, cartesianIndex);
820 }
821
822 OPM_END_PARALLEL_TRY_CATCH("EclWriter::captureLocalFluxData() failed: ",
823 this->simulator_.vanguard().grid().comm())
824
825 this->outputModule_->finalizeFluxData();
826 }
827
828 Simulator& simulator_;
829 std::unique_ptr<OutputModule> outputModule_;
830 Scalar restartTimeStepSize_;
831 int rank_ ;
832 Inplace inplace_;
833};
834
835} // namespace Opm
836
837#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