OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-01 21:39:09 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3708 from joakim-hove/tracer-rst-init

Browse Source 
Tracer rst init
This commit is contained in:
Joakim Hove 2021-11-30 14:08:46 +01:00 committed by GitHub
commit 90f2445557
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
7 changed files with 79 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
ebos/eclgenericoutputblackoilmodule.cc
View File

@ -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

53
ebos/eclgenerictracermodel.cc
View File

@ -106,15 +106,6 @@ EclGenericTracerModel(const GridView& gridView,
{
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
const std::string& EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
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<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
Scalar EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
@ -126,9 +117,39 @@ tracerConcentration(int tracerIdx, int globalDofIdx) const
return tracerConcentration_[tracerIdx][globalDofIdx];
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
void EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
doInit(bool enabled, size_t numGridDof,
setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value)
{
this->tracerConcentration_[tracerIdx][globalDofIdx] = value;
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
int EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
numTracers() const
{
return this->eclState_.tracer().size();
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
std::string EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
fname(int tracerIdx) const
{
return this->eclState_.tracer()[tracerIdx].fname();
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
const std::string& EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
name(int tracerIdx) const
{
return this->eclState_.tracer()[tracerIdx].name;
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
void EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
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

10
ebos/eclgenerictracermodel.hh
View File

@ -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<std::string> tracerNames_;
std::vector<int> tracerPhaseIdx_;
std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> tracerConcentration_;
std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> tracerConcentrationInitial_;

11
ebos/eclproblem.hh
View File

@ -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<TypeTag>& tracerModel() const
{ return tracerModel_; }
EclTracerModel<TypeTag>& 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.

24
ebos/ecltracermodel.hh
View File

@ -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<double> 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; tracerIdx<this->tracerPhaseIdx_.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]);

16
ebos/eclwriter.hh
View File

@ -305,6 +305,12 @@ public:
std::vector<RestartKey> 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&>(simulator_);
auto& thpres = mutableSimulator.problem().thresholdPressure();

2
ebos/vtkecltracermodule.hh
View File

@ -153,7 +153,7 @@ namespace Opm {
if (eclTracerConcentrationOutput_()){
const auto& tracerModel = this->simulator_.problem().tracerModel();
for(size_t tracerIdx=0; tracerIdx<eclTracerConcentration_.size();++tracerIdx){
const std::string tmp = "tracerConcentration_" + tracerModel.tracerName(tracerIdx);
const std::string tmp = "tracerConcentration_" + tracerModel.name(tracerIdx);
this->commitScalarBuffer_(baseWriter,tmp.c_str(), eclTracerConcentration_[tracerIdx]);
}
}