diff --git a/ebos/eclgenericoutputblackoilmodule.cc b/ebos/eclgenericoutputblackoilmodule.cc index 820338a7b..e7f465a74 100644 --- a/ebos/eclgenericoutputblackoilmodule.cc +++ b/ebos/eclgenericoutputblackoilmodule.cc @@ -668,10 +668,9 @@ assignToSolution(data::Solution& sol) // tracers if (!tracerConcentrations_.empty()) { const auto& tracers = eclState_.tracer(); - size_t tracerIdx = 0; - for (const auto& tracer : tracers) { - std::string tmp = tracer.name + "F"; - sol.insert(tmp, UnitSystem::measure::identity, std::move(tracerConcentrations_[tracerIdx++]), data::TargetType::RESTART_TRACER_SOLUTION); + for (std::size_t tracerIdx = 0; tracerIdx < tracers.size(); tracerIdx++) { + const auto& tracer = tracers[tracerIdx]; + sol.insert(tracer.fname(), UnitSystem::measure::identity, std::move(tracerConcentrations_[tracerIdx]), data::TargetType::RESTART_TRACER_SOLUTION); } // We need put tracerConcentrations into a valid state. // Otherwise next time we end up here outside of a restart write we will diff --git a/ebos/eclgenerictracermodel.cc b/ebos/eclgenerictracermodel.cc index 61ae6fcc3..f7d5b6d35 100644 --- a/ebos/eclgenerictracermodel.cc +++ b/ebos/eclgenerictracermodel.cc @@ -106,15 +106,6 @@ EclGenericTracerModel(const GridView& gridView, { } -template -const std::string& EclGenericTracerModel:: -tracerName(int tracerIdx) const -{ - if (tracerNames_.empty()) - throw std::logic_error("This method should never be called when there are no tracers in the model"); - - return tracerNames_[tracerIdx]; -} template Scalar EclGenericTracerModel:: @@ -126,9 +117,39 @@ tracerConcentration(int tracerIdx, int globalDofIdx) const return tracerConcentration_[tracerIdx][globalDofIdx]; } + template void EclGenericTracerModel:: -doInit(bool enabled, size_t numGridDof, +setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value) +{ + this->tracerConcentration_[tracerIdx][globalDofIdx] = value; +} + +template +int EclGenericTracerModel:: +numTracers() const +{ + return this->eclState_.tracer().size(); +} + +template +std::string EclGenericTracerModel:: +fname(int tracerIdx) const +{ + return this->eclState_.tracer()[tracerIdx].fname(); +} + +template +const std::string& EclGenericTracerModel:: +name(int tracerIdx) const +{ + return this->eclState_.tracer()[tracerIdx].name; +} + + +template +void EclGenericTracerModel:: +doInit(bool enabled, bool rst, size_t numGridDof, size_t gasPhaseIdx, size_t oilPhaseIdx, size_t waterPhaseIdx) { const auto& tracers = eclState_.tracer(); @@ -148,15 +169,14 @@ doInit(bool enabled, size_t numGridDof, // retrieve the number of tracers from the deck const size_t numTracers = tracers.size(); - tracerNames_.resize(numTracers); tracerConcentration_.resize(numTracers); storageOfTimeIndex1_.resize(numTracers); // the phase where the tracer is tracerPhaseIdx_.resize(numTracers); - size_t tracerIdx = 0; - for (const auto& tracer : tracers) { - tracerNames_[tracerIdx] = tracer.name; + for (size_t tracerIdx = 0; tracerIdx < numTracers; tracerIdx++) { + const auto& tracer = tracers[tracerIdx]; + if (tracer.phase == Phase::WATER) tracerPhaseIdx_[tracerIdx] = waterPhaseIdx; else if (tracer.phase == Phase::OIL) @@ -167,6 +187,9 @@ doInit(bool enabled, size_t numGridDof, tracerConcentration_[tracerIdx].resize(numGridDof); storageOfTimeIndex1_[tracerIdx].resize(numGridDof); + if (rst) + continue; + //TBLK keyword if (tracer.free_concentration.has_value()){ @@ -190,8 +213,6 @@ doInit(bool enabled, size_t numGridDof, } } else throw std::logic_error(fmt::format("Can not initialize tracer: {}", tracer.name)); - - ++tracerIdx; } // initial tracer concentration diff --git a/ebos/eclgenerictracermodel.hh b/ebos/eclgenerictracermodel.hh index aa02163e5..90efc840e 100644 --- a/ebos/eclgenerictracermodel.hh +++ b/ebos/eclgenerictracermodel.hh @@ -55,18 +55,20 @@ public: /*! * \brief Return the number of tracers considered by the tracerModel. */ - int numTracers() const - { return tracerNames_.size(); } + int numTracers() const; /*! * \brief Return the tracer name */ - const std::string& tracerName(int tracerIdx) const; + const std::string& name(int tracerIdx) const; + std::string fname(int tracerIdx) const; + /*! * \brief Return the tracer concentration for tracer index and global DofIdx */ Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const; + void setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value); /*! * \brief Return well tracer rates @@ -85,6 +87,7 @@ protected: * \brief Initialize all internal data structures needed by the tracer module */ void doInit(bool enabled, + bool rst, size_t numGridDof, size_t gasPhaseIdx, size_t oilPhaseIdx, @@ -99,7 +102,6 @@ protected: const CartesianIndexMapper& cartMapper_; const DofMapper& dofMapper_; - std::vector tracerNames_; std::vector tracerPhaseIdx_; std::vector>> tracerConcentration_; std::vector>> tracerConcentrationInitial_; diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index a8da3ccd8..05e464def 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -875,6 +875,7 @@ public: transmissibilities_.finishInit(); const auto& initconfig = eclState.getInitConfig(); + tracerModel_.init(initconfig.restartRequested()); if (initconfig.restartRequested()) readEclRestartSolution_(); else @@ -889,8 +890,6 @@ public: this->maxPolymerAdsorption_.resize(numElements, 0.0); } - tracerModel_.init(); - readBoundaryConditions_(); if (enableDriftCompensation_) { @@ -1465,6 +1464,9 @@ public: const EclTracerModel& tracerModel() const { return tracerModel_; } + EclTracerModel& tracerModel() + { return tracerModel_; } + /*! * \copydoc FvBaseMultiPhaseProblem::porosity * @@ -2475,11 +2477,6 @@ private: for (size_t pvtRegionIdx = 0; pvtRegionIdx < this->maxDRv_.size(); ++pvtRegionIdx) this->maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*simulator.timeStepSize(); - if (tracerModel().numTracers() > 0 && this->gridView().comm().rank() == 0) - std::cout << "Warning: Restart is not implemented for the tracer model, it will initialize itself " - << "with the initial tracer concentration.\n" - << std::flush; - // assign the restart solution to the current solution. note that we still need // to compute real initial solution after this because the initial fluid states // need to be correct for stuff like boundary conditions. diff --git a/ebos/ecltracermodel.hh b/ebos/ecltracermodel.hh index 597c4705f..0718060b1 100644 --- a/ebos/ecltracermodel.hh +++ b/ebos/ecltracermodel.hh @@ -103,10 +103,10 @@ public: /*! * \brief Initialize all internal data structures needed by the tracer module */ - void init() + void init(bool rst) { bool enabled = EWOMS_GET_PARAM(TypeTag, bool, EnableTracerModel); - this->doInit(enabled, simulator_.model().numGridDof(), + this->doInit(enabled, rst, simulator_.model().numGridDof(), gasPhaseIdx, oilPhaseIdx, waterPhaseIdx); prepareTracerBatches(); @@ -306,12 +306,12 @@ protected: const auto& well = wellPtr->wellEcl(); for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) { - this->wellTracerRate_[std::make_pair(well.name(),this->tracerNames_[tr.idx_[tIdx]])] = 0.0; + this->wellTracerRate_[std::make_pair(well.name(), this->name(tr.idx_[tIdx]))] = 0.0; } std::vector wtracer(tr.numTracer()); for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) { - wtracer[tIdx] = well.getTracerProperties().getConcentration(this->tracerNames_[tr.idx_[tIdx]]); + wtracer[tIdx] = well.getTracerProperties().getConcentration(this->name(tr.idx_[tIdx])); } for (auto& perfData : wellPtr->perforationData()) { @@ -321,7 +321,7 @@ protected: for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) { tr.residual_[tIdx][I][0] -= rate*wtracer[tIdx]; // Store _injector_ tracer rate for reporting - this->wellTracerRate_.at(std::make_pair(well.name(),this->tracerNames_[tr.idx_[tIdx]])) += rate*wtracer[tIdx]; + this->wellTracerRate_.at(std::make_pair(well.name(),this->name(tr.idx_[tIdx]))) += rate*wtracer[tIdx]; } } else if (rate < 0) { @@ -402,7 +402,7 @@ protected: if (rate < 0) { rateWellNeg += rate; for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) { - this->wellTracerRate_.at(std::make_pair(well.name(),this->tracerNames_[tr.idx_[tIdx]])) += rate*tr.concentration_[tIdx][I]; + this->wellTracerRate_.at(std::make_pair(well.name(),this->name(tr.idx_[tIdx]))) += rate*tr.concentration_[tIdx][I]; } } else { @@ -424,7 +424,7 @@ protected: const Scalar bucketPrDay = 10.0/(1000.*3600.*24.); // ... keeps (some) trouble away const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal/rateWellNeg : 0.0; for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) { - this->wellTracerRate_.at(std::make_pair(well.name(),this->tracerNames_[tr.idx_[tIdx]])) *= factor; + this->wellTracerRate_.at(std::make_pair(well.name(),this->name(tr.idx_[tIdx]))) *= factor; } } } @@ -435,27 +435,27 @@ protected: for (size_t tracerIdx=0; tracerIdxtracerPhaseIdx_.size(); ++tracerIdx) { if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::waterPhaseIdx) { if (! FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)){ - throw std::runtime_error("Water tracer specified for non-water fluid system:" + this->tracerNames_[tracerIdx]); + throw std::runtime_error("Water tracer specified for non-water fluid system:" + this->name(tracerIdx)); } wat_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]); } else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::oilPhaseIdx) { if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){ - throw std::runtime_error("Oil tracer specified for non-oil fluid system:" + this->tracerNames_[tracerIdx]); + throw std::runtime_error("Oil tracer specified for non-oil fluid system:" + this->name(tracerIdx)); } if (FluidSystem::enableVaporizedOil()) { - throw std::runtime_error("Oil tracer in combination with kw VAPOIL is not supported: " + this->tracerNames_[tracerIdx]); + throw std::runtime_error("Oil tracer in combination with kw VAPOIL is not supported: " + this->name(tracerIdx)); } oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]); } else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::gasPhaseIdx) { if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){ - throw std::runtime_error("Gas tracer specified for non-gas fluid system:" + this->tracerNames_[tracerIdx]); + throw std::runtime_error("Gas tracer specified for non-gas fluid system:" + this->name(tracerIdx)); } if (FluidSystem::enableDissolvedGas()) { - throw std::runtime_error("Gas tracer in combination with kw DISGAS is not supported: " + this->tracerNames_[tracerIdx]); + throw std::runtime_error("Gas tracer in combination with kw DISGAS is not supported: " + this->name(tracerIdx)); } gas_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]); diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index b0cc96ae5..559e04a18 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -305,6 +305,12 @@ public: std::vector extraKeys = {{"OPMEXTRA", UnitSystem::measure::identity, false}, {"THRESHPR", UnitSystem::measure::pressure, inputThpres.active()}}; + { + const auto& tracers = simulator_.vanguard().eclState().tracer(); + for (const auto& tracer : tracers) + solutionKeys.emplace_back(tracer.fname(), UnitSystem::measure::identity, true); + } + // The episodeIndex is rewined one back before beginRestart is called // and can not be used here. // We just ask the initconfig directly to be sure that we use the correct @@ -326,6 +332,16 @@ public: eclOutputModule_->setRestart(restartValues.solution, elemIdx, globalIdx); } + auto& tracer_model = simulator_.problem().tracerModel(); + for (int tracer_index = 0; tracer_index < tracer_model.numTracers(); tracer_index++) { + const auto& tracer_name = tracer_model.fname(tracer_index); + const auto& tracer_solution = restartValues.solution.data(tracer_name); + for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { + unsigned globalIdx = this->collectToIORank_.localIdxToGlobalIdx(elemIdx); + tracer_model.setTracerConcentration(tracer_index, globalIdx, tracer_solution[globalIdx]); + } + } + if (inputThpres.active()) { Simulator& mutableSimulator = const_cast(simulator_); auto& thpres = mutableSimulator.problem().thresholdPressure(); diff --git a/ebos/vtkecltracermodule.hh b/ebos/vtkecltracermodule.hh index 78263b147..b1b2ec608 100644 --- a/ebos/vtkecltracermodule.hh +++ b/ebos/vtkecltracermodule.hh @@ -153,7 +153,7 @@ namespace Opm { if (eclTracerConcentrationOutput_()){ const auto& tracerModel = this->simulator_.problem().tracerModel(); for(size_t tracerIdx=0; tracerIdxcommitScalarBuffer_(baseWriter,tmp.c_str(), eclTracerConcentration_[tracerIdx]); } }