diff --git a/opm/simulators/flow/EclWriter.hpp b/opm/simulators/flow/EclWriter.hpp index 6ebb2109765..6138b6d1863 100644 --- a/opm/simulators/flow/EclWriter.hpp +++ b/opm/simulators/flow/EclWriter.hpp @@ -500,8 +500,12 @@ class EclWriter : public EclGenericWriter { const auto& tracers = simulator_.vanguard().eclState().tracer(); - for (const auto& tracer : tracers) + for (const auto& tracer : tracers) { + bool enableSolTracer = (tracer.phase == Phase::GAS && FluidSystem::enableDissolvedGas()) || + (tracer.phase == Phase::OIL && FluidSystem::enableVaporizedOil()); solutionKeys.emplace_back(tracer.fname(), UnitSystem::measure::identity, true); + solutionKeys.emplace_back(tracer.sname(), UnitSystem::measure::identity, enableSolTracer); + } } // The episodeIndex is rewined one back before beginRestart is called @@ -527,11 +531,30 @@ class EclWriter : public EclGenericWriter 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.template data(tracer_name); + // Free tracers + const auto& free_tracer_name = tracer_model.fname(tracer_index); + const auto& free_tracer_solution = restartValues.solution.template data(free_tracer_name); for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { unsigned globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx); - tracer_model.setTracerConcentration(tracer_index, globalIdx, tracer_solution[globalIdx]); + tracer_model.setFreeTracerConcentration(tracer_index, globalIdx, free_tracer_solution[globalIdx]); + } + // Solution tracer (only if DISGAS/VAPOIL are active for gas/oil tracers) + if ((tracer_model.phase(tracer_index) == Phase::GAS && FluidSystem::enableDissolvedGas()) || + (tracer_model.phase(tracer_index) == Phase::OIL && FluidSystem::enableVaporizedOil())) { + tracer_model.setEnableSolTracers(tracer_index, true); + const auto& sol_tracer_name = tracer_model.sname(tracer_index); + const auto& sol_tracer_solution = restartValues.solution.template data(sol_tracer_name); + for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { + unsigned globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx); + tracer_model.setSolTracerConcentration(tracer_index, globalIdx, sol_tracer_solution[globalIdx]); + } + } + else { + tracer_model.setEnableSolTracers(tracer_index, false); + for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { + unsigned globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx); + tracer_model.setSolTracerConcentration(tracer_index, globalIdx, 0.0); + } } } diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index c82582bcbf5..1dcfd599f0d 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -700,21 +700,40 @@ assignToSolution(data::Solution& sol) } // Tracers - if (! this->tracerConcentrations_.empty()) { + if (! this->freeTracerConcentrations_.empty()) { const auto& tracers = this->eclState_.tracer(); for (auto tracerIdx = 0*tracers.size(); tracerIdx < tracers.size(); ++tracerIdx) { sol.insert(tracers[tracerIdx].fname(), UnitSystem::measure::identity, - std::move(tracerConcentrations_[tracerIdx]), + std::move(freeTracerConcentrations_[tracerIdx]), data::TargetType::RESTART_TRACER_SOLUTION); } - // Put tracerConcentrations container into a valid state. Otherwise + // Put freeTracerConcentrations container into a valid state. Otherwise // we'll move from vectors that have already been moved from if we // get here and it's not a restart step. - this->tracerConcentrations_.clear(); + this->freeTracerConcentrations_.clear(); + } + if (! this->solTracerConcentrations_.empty()) { + const auto& tracers = this->eclState_.tracer(); + for (auto tracerIdx = 0*tracers.size(); + tracerIdx < tracers.size(); ++tracerIdx) + { + if (solTracerConcentrations_[tracerIdx].empty()) + continue; + + sol.insert(tracers[tracerIdx].sname(), + UnitSystem::measure::identity, + std::move(solTracerConcentrations_[tracerIdx]), + data::TargetType::RESTART_TRACER_SOLUTION); + } + + // Put solTracerConcentrations container into a valid state. Otherwise + // we'll move from vectors that have already been moved from if we + // get here and it's not a restart step. + this->solTracerConcentrations_.clear(); } } @@ -838,6 +857,7 @@ doAllocBuffers(const unsigned bufferSize, const bool vapparsActive, const bool enableHysteresis, const unsigned numTracers, + const std::vector& enableSolTracers, const unsigned numOutputNnc) { // Output RESTART_OPM_EXTENDED only when explicitly requested by user. @@ -1304,10 +1324,16 @@ doAllocBuffers(const unsigned bufferSize, // tracers if (numTracers > 0) { - tracerConcentrations_.resize(numTracers); + freeTracerConcentrations_.resize(numTracers); + for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) + { + freeTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0); + } + solTracerConcentrations_.resize(numTracers); for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) { - tracerConcentrations_[tracerIdx].resize(bufferSize, 0.0); + if (enableSolTracers[tracerIdx]) + solTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0); } } diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index 4eb6216944a..914c9a4a8d2 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -340,6 +340,7 @@ class GenericOutputBlackoilModule { const bool vapparsActive, const bool enableHysteresis, unsigned numTracers, + const std::vector& enableSolTracers, unsigned numOutputNnc); void makeRegionSum(Inplace& inplace, @@ -504,7 +505,8 @@ class GenericOutputBlackoilModule { std::array viscosity_; std::array relativePermeability_; - std::vector tracerConcentrations_; + std::vector freeTracerConcentrations_; + std::vector solTracerConcentrations_; std::array residual_; diff --git a/opm/simulators/flow/GenericTracerModel.cpp b/opm/simulators/flow/GenericTracerModel.cpp index 0076847920c..f58ed53d34a 100644 --- a/opm/simulators/flow/GenericTracerModel.cpp +++ b/opm/simulators/flow/GenericTracerModel.cpp @@ -24,6 +24,8 @@ #include +#include + #if HAVE_DUNE_FEM #include #include @@ -39,6 +41,7 @@ template class GenericTracerModel>, Dune::MultipleCodimMultipleGeomTypeMapper>>, Opm::EcfvStencil>,false,false>, + BlackOilFluidSystem, double>; #if HAVE_DUNE_FEM @@ -50,12 +53,14 @@ template class GenericTracerModel, EcfvStencil, + BlackOilFluidSystem, double>; #else template class GenericTracerModel>>, Dune::MultipleCodimMultipleGeomTypeMapper>>>, EcfvStencil>>,false,false>, + BlackOilFluidSystem, double>; template class GenericTracerModel >, @@ -65,6 +70,7 @@ template class GenericTracerModel >, false, false>, + BlackOilFluidSystem, double>; #endif #endif // HAVE_DUNE_FEM diff --git a/opm/simulators/flow/GenericTracerModel.hpp b/opm/simulators/flow/GenericTracerModel.hpp index af456dfc70e..312f7b06eb7 100644 --- a/opm/simulators/flow/GenericTracerModel.hpp +++ b/opm/simulators/flow/GenericTracerModel.hpp @@ -36,6 +36,8 @@ #include +#include + #include #include #include @@ -49,11 +51,12 @@ namespace Opm { class EclipseState; class Well; -template +template class GenericTracerModel { public: - using TracerMatrix = Dune::BCRSMatrix>; - using TracerVector = Dune::BlockVector>; + using TracerVectorSingle = Dune::BlockVector>; + using TracerMatrix = Dune::BCRSMatrix>; + using TracerVector = Dune::BlockVector>; using CartesianIndexMapper = Dune::CartesianIndexMapper; static constexpr int dimWorld = Grid::dimensionworld; /*! @@ -66,25 +69,44 @@ class GenericTracerModel { */ const std::string& name(int tracerIdx) const; std::string fname(int tracerIdx) const; + std::string sname(int tracerIdx) const; + std::string wellfname(int tracerIdx) const; + std::string wellsname(int tracerIdx) const; + Phase phase(int tracerIdx) const; + const std::vector& enableSolTracers() 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); + Scalar freeTracerConcentration(int tracerIdx, int globalDofIdx) const; + Scalar solTracerConcentration(int tracerIdx, int globalDofIdx) const; + void setFreeTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value); + void setSolTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value); + void setEnableSolTracers(int tracerIdx, bool enableSolTracer); /*! * \brief Return well tracer rates */ const std::map, Scalar>& getWellTracerRates() const {return wellTracerRate_;} + const std::map, Scalar>& + getWellFreeTracerRates() const {return wellFreeTracerRate_;} + const std::map, Scalar>& + getWellSolTracerRates() const {return wellSolTracerRate_;} + const std::map, Scalar>& + getMswTracerRates() const {return mSwTracerRate_;} template void serializeOp(Serializer& serializer) { serializer(tracerConcentration_); + serializer(freeTracerConcentration_); + serializer(solTracerConcentration_); serializer(wellTracerRate_); + serializer(wellFreeTracerRate_); + serializer(wellSolTracerRate_); + serializer(mSwTracerRate_); } protected: @@ -117,12 +139,20 @@ class GenericTracerModel { const DofMapper& dofMapper_; std::vector tracerPhaseIdx_; - std::vector>> tracerConcentration_; + std::vector enableSolTracers_; + std::vector tracerConcentration_; std::unique_ptr tracerMatrix_; - std::vector>> storageOfTimeIndex1_; + std::vector freeTracerConcentration_; + std::vector solTracerConcentration_; - // -> wellRate + // -> wellRate std::map, Scalar> wellTracerRate_; + std::map, Scalar> wellFreeTracerRate_; + std::map, Scalar> wellSolTracerRate_; + + // -> wellRate + std::map, Scalar> mSwTracerRate_; + /// \brief Function returning the cell centers std::function(int)> centroids_; }; diff --git a/opm/simulators/flow/GenericTracerModel_impl.hpp b/opm/simulators/flow/GenericTracerModel_impl.hpp index 4707c05d1f7..5d3d149d9bc 100644 --- a/opm/simulators/flow/GenericTracerModel_impl.hpp +++ b/opm/simulators/flow/GenericTracerModel_impl.hpp @@ -39,7 +39,6 @@ #include #include -#include #include #include #include @@ -97,8 +96,8 @@ createParallelFlexibleSolver(const Dune::CpGrid& grid, const Matrix& M, const Pr } #endif -template -GenericTracerModel:: +template +GenericTracerModel:: GenericTracerModel(const GridView& gridView, const EclipseState& eclState, const CartesianIndexMapper& cartMapper, @@ -112,53 +111,114 @@ GenericTracerModel(const GridView& gridView, { } -template -Scalar GenericTracerModel:: -tracerConcentration(int tracerIdx, int globalDofIdx) const +template +Scalar GenericTracerModel:: +freeTracerConcentration(int tracerIdx, int globalDofIdx) const { - if (tracerConcentration_.empty()) + if (freeTracerConcentration_.empty()) return 0.0; - return tracerConcentration_[tracerIdx][globalDofIdx]; + return freeTracerConcentration_[tracerIdx][globalDofIdx]; } -template -void GenericTracerModel:: -setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value) +template +Scalar GenericTracerModel:: +solTracerConcentration(int tracerIdx, int globalDofIdx) const { - this->tracerConcentration_[tracerIdx][globalDofIdx] = value; + if (solTracerConcentration_.empty()) + return 0.0; + + return solTracerConcentration_[tracerIdx][globalDofIdx]; +} + +template +void GenericTracerModel:: +setFreeTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value) +{ + this->freeTracerConcentration_[tracerIdx][globalDofIdx] = value; + this->tracerConcentration_[tracerIdx][globalDofIdx][0] = value; +} + +template +void GenericTracerModel:: +setSolTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value) +{ + this->solTracerConcentration_[tracerIdx][globalDofIdx] = value; + this->tracerConcentration_[tracerIdx][globalDofIdx][1] = value; +} + +template +void GenericTracerModel:: +setEnableSolTracers(int tracerIdx, bool enableSolTracer) +{ + this->enableSolTracers_[tracerIdx] = enableSolTracer; } -template -int GenericTracerModel:: +template +int GenericTracerModel:: numTracers() const { return this->eclState_.tracer().size(); } -template -std::string GenericTracerModel:: +template +std::string GenericTracerModel:: fname(int tracerIdx) const { return this->eclState_.tracer()[tracerIdx].fname(); } -template -Scalar GenericTracerModel:: +template +std::string GenericTracerModel:: +sname(int tracerIdx) const +{ + return this->eclState_.tracer()[tracerIdx].sname(); +} + +template +std::string GenericTracerModel:: +wellfname(int tracerIdx) const +{ + return this->eclState_.tracer()[tracerIdx].wellfname(); +} + +template +std::string GenericTracerModel:: +wellsname(int tracerIdx) const +{ + return this->eclState_.tracer()[tracerIdx].wellsname(); +} + +template +Phase GenericTracerModel:: +phase(int tracerIdx) const +{ + return this->eclState_.tracer()[tracerIdx].phase; +} + +template +const std::vector& GenericTracerModel:: +enableSolTracers() const +{ + return this->enableSolTracers_; +} + +template +Scalar GenericTracerModel:: currentConcentration_(const Well& eclWell, const std::string& name) const { return eclWell.getTracerProperties().getConcentration(name); } -template -const std::string& GenericTracerModel:: +template +const std::string& GenericTracerModel:: name(int tracerIdx) const { return this->eclState_.tracer()[tracerIdx].name; } -template -void GenericTracerModel:: +template +void GenericTracerModel:: doInit(bool rst, std::size_t numGridDof, std::size_t gasPhaseIdx, std::size_t oilPhaseIdx, std::size_t waterPhaseIdx) { @@ -169,8 +229,10 @@ doInit(bool rst, std::size_t numGridDof, // retrieve the number of tracers from the deck const std::size_t numTracers = tracers.size(); + enableSolTracers_.resize(numTracers); tracerConcentration_.resize(numTracers); - storageOfTimeIndex1_.resize(numTracers); + freeTracerConcentration_.resize(numTracers); + solTracerConcentration_.resize(numTracers); // the phase where the tracer is tracerPhaseIdx_.resize(numTracers); @@ -185,34 +247,97 @@ doInit(bool rst, std::size_t numGridDof, tracerPhaseIdx_[tracerIdx] = gasPhaseIdx; tracerConcentration_[tracerIdx].resize(numGridDof); - storageOfTimeIndex1_[tracerIdx].resize(numGridDof); + freeTracerConcentration_[tracerIdx].resize(numGridDof); + solTracerConcentration_[tracerIdx].resize(numGridDof); if (rst) continue; - //TBLK keyword + // TBLKF keyword if (tracer.free_concentration.has_value()){ const auto& free_concentration = tracer.free_concentration.value(); int tblkDatasize = free_concentration.size(); if (tblkDatasize < cartMapper_.cartesianSize()){ - throw std::runtime_error("Wrong size of TBLK for" + tracer.name); + throw std::runtime_error("Wrong size of TBLKF for" + tracer.name); } for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) { int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx); - tracerConcentration_[tracerIdx][globalDofIdx] = free_concentration[cartDofIdx]; + tracerConcentration_[tracerIdx][globalDofIdx][0] = free_concentration[cartDofIdx]; + freeTracerConcentration_[tracerIdx][globalDofIdx] = free_concentration[cartDofIdx]; } } - //TVDPF keyword + // TVDPF keyword else if (tracer.free_tvdp.has_value()) { const auto& free_tvdp = tracer.free_tvdp.value(); for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) { - tracerConcentration_[tracerIdx][globalDofIdx] = + tracerConcentration_[tracerIdx][globalDofIdx][0] = free_tvdp.evaluate("TRACER_CONCENTRATION", centroids_(globalDofIdx)[2]); + freeTracerConcentration_[tracerIdx][globalDofIdx] = + free_tvdp.evaluate("TRACER_CONCENTRATION", + centroids_(globalDofIdx)[2]); + } + } + else { + OpmLog::warning(fmt::format("No TBLKF or TVDPF given for free tracer {}. " + "Initial values set to zero. ", tracer.name)); + for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) { + tracerConcentration_[tracerIdx][globalDofIdx][0] = 0.0; + freeTracerConcentration_[tracerIdx][globalDofIdx] = 0.0; + } + } + + // Solution tracer initialization only needed for gas/oil tracers with DISGAS/VAPOIL active + if (tracer.phase != Phase::WATER && + ((tracer.phase == Phase::GAS && FluidSystem::enableDissolvedGas()) || + (tracer.phase == Phase::OIL && FluidSystem::enableVaporizedOil()))) { + // TBLKS keyword + if (tracer.solution_concentration.has_value()){ + enableSolTracers_[tracerIdx] = true; + const auto& solution_concentration = tracer.solution_concentration.value(); + int tblkDatasize = solution_concentration.size(); + if (tblkDatasize < cartMapper_.cartesianSize()){ + throw std::runtime_error("Wrong size of TBLKS for" + tracer.name); + } + for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) { + int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx); + tracerConcentration_[tracerIdx][globalDofIdx][1] = solution_concentration[cartDofIdx]; + solTracerConcentration_[tracerIdx][globalDofIdx] = solution_concentration[cartDofIdx]; + } + } + // TVDPS keyword + else if (tracer.solution_tvdp.has_value()) { + enableSolTracers_[tracerIdx] = true; + const auto& solution_tvdp = tracer.solution_tvdp.value(); + for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) { + tracerConcentration_[tracerIdx][globalDofIdx][1] = + solution_tvdp.evaluate("TRACER_CONCENTRATION", + centroids_(globalDofIdx)[2]); + solTracerConcentration_[tracerIdx][globalDofIdx] = + solution_tvdp.evaluate("TRACER_CONCENTRATION", + centroids_(globalDofIdx)[2]); + } + } + else { + // No solution tracers, default to zero + enableSolTracers_[tracerIdx] = false; + OpmLog::warning(fmt::format("No TBLKS or TVDPS given for solution tracer {}. " + "Initial values set to zero. ", tracer.name)); + for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) { + tracerConcentration_[tracerIdx][globalDofIdx][1] = 0.0; + solTracerConcentration_[tracerIdx][globalDofIdx] = 0.0; + } + } + } + else { + // No solution tracers, default to zero + enableSolTracers_[tracerIdx] = false; + for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) { + tracerConcentration_[tracerIdx][globalDofIdx][1] = 0.0; + solTracerConcentration_[tracerIdx][globalDofIdx] = 0.0; } - } else - throw std::logic_error(fmt::format("Can not initialize tracer: {}", tracer.name)); + } } // allocate matrix for storing the Jacobian of the tracer residual @@ -237,8 +362,9 @@ doInit(bool rst, std::size_t numGridDof, } // allocate space for the rows of the matrix - for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) + for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) { tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size()); + } tracerMatrix_->endrowsizes(); // fill the rows with indices. each degree of freedom talks to @@ -247,14 +373,15 @@ doInit(bool rst, std::size_t numGridDof, for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) { typename NeighborSet::iterator nIt = neighbors[dofIdx].begin(); typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end(); - for (; nIt != nEndIt; ++nIt) + for (; nIt != nEndIt; ++nIt) { tracerMatrix_->addindex(dofIdx, *nIt); + } } tracerMatrix_->endindices(); } -template -bool GenericTracerModel:: +template +bool GenericTracerModel:: linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b) { x = 0.0; @@ -308,8 +435,8 @@ linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b) #endif } -template -bool GenericTracerModel:: +template +bool GenericTracerModel:: linearSolveBatchwise_(const TracerMatrix& M, std::vector& x, std::vector& b) { Scalar tolerance = 1e-2; diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 9ba8f35ae6d..8a87873e1a1 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -240,6 +240,7 @@ class OutputBlackOilModule : public GenericOutputBlackoilModuleenableHysteresis(), problem.tracerModel().numTracers(), + problem.tracerModel().enableSolTracers(), problem.eclWriter()->getOutputNnc().size()); } @@ -661,14 +662,23 @@ class OutputBlackOilModule : public GenericOutputBlackoilModuletracerConcentrations_.empty()) { + if (! this->freeTracerConcentrations_.empty()) { for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) { - if (this->tracerConcentrations_[tracerIdx].empty()) { + if (this->freeTracerConcentrations_[tracerIdx].empty()) { continue; } - - this->tracerConcentrations_[tracerIdx][globalDofIdx] = - tracerModel.tracerConcentration(tracerIdx, globalDofIdx); + this->freeTracerConcentrations_[tracerIdx][globalDofIdx] = + tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx); + } + } + if (! this->solTracerConcentrations_.empty()) { + for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) { + if (this->solTracerConcentrations_[tracerIdx].empty()) { + continue; + } + this->solTracerConcentrations_[tracerIdx][globalDofIdx] = + tracerModel.solTracerConcentration(tracerIdx, globalDofIdx); + } } diff --git a/opm/simulators/flow/TracerModel.hpp b/opm/simulators/flow/TracerModel.hpp index df53a119411..516d222083d 100644 --- a/opm/simulators/flow/TracerModel.hpp +++ b/opm/simulators/flow/TracerModel.hpp @@ -63,12 +63,14 @@ class TracerModel : public GenericTracerModel, GetPropType, GetPropType, + GetPropType, GetPropType> { using BaseType = GenericTracerModel, GetPropType, GetPropType, GetPropType, + GetPropType, GetPropType>; using Simulator = GetPropType; using GridView = GetPropType; @@ -113,7 +115,7 @@ class TracerModel : public GenericTracerModelname(tracerIdx)); } - if (FluidSystem::enableVaporizedOil()) { - throw std::runtime_error("Oil tracer in combination with kw VAPOIL is not supported: " + this->name(tracerIdx)); - } oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]); } @@ -154,12 +153,15 @@ class TracerModel : public GenericTracerModelname(tracerIdx)); } - if (FluidSystem::enableDissolvedGas()) { - throw std::runtime_error("Gas tracer in combination with kw DISGAS is not supported: " + this->name(tracerIdx)); - } gas_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]); } + + // resize free and solution volume storages + fVol1_[this->tracerPhaseIdx_[tracerIdx]].resize(this->freeTracerConcentration_[tracerIdx].size()); + sVol1_[this->tracerPhaseIdx_[tracerIdx]].resize(this->solTracerConcentration_[tracerIdx].size()); + dsVol_[this->tracerPhaseIdx_[tracerIdx]].resize(this->solTracerConcentration_[tracerIdx].size()); + dfVol_[this->tracerPhaseIdx_[tracerIdx]].resize(this->solTracerConcentration_[tracerIdx].size()); } // will be valid after we move out of tracerMatrix_ @@ -220,39 +222,63 @@ class TracerModel : public GenericTracerModel - void computeVolume_(LhsEval& freeVolume, - const int tracerPhaseIdx, - const ElementContext& elemCtx, - unsigned scvIdx, - unsigned timeIdx) + // compute volume associated with free concentration + Scalar computeFreeVolume_(const int tracerPhaseIdx, + unsigned globalDofIdx, + unsigned timeIdx) { - const auto& intQuants = elemCtx.intensiveQuantities(scvIdx, timeIdx); + const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx); const auto& fs = intQuants.fluidState(); + Scalar phaseVolume = decay(fs.saturation(tracerPhaseIdx)) *decay(fs.invB(tracerPhaseIdx)) *decay(intQuants.porosity()); - // avoid singular matrix if no water is present. - phaseVolume = max(phaseVolume, 1e-10); + return max(phaseVolume, 1e-10); + } + + // compute volume associated with solution concentration + Scalar computeSolutionVolume_(const int tracerPhaseIdx, + unsigned globalDofIdx, + unsigned timeIdx) + { + const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx); + const auto& fs = intQuants.fluidState(); + + Scalar phaseVolume; - if (std::is_same::value) - freeVolume = phaseVolume; - else - freeVolume = phaseVolume * variable(1.0, 0); + // vaporized oil + if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) { + phaseVolume = + decay(fs.saturation(FluidSystem::gasPhaseIdx)) + * decay(fs.invB(FluidSystem::gasPhaseIdx)) + * decay(fs.Rv()) + * decay(intQuants.porosity()); + } - } + // dissolved gas + else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) { + phaseVolume = + decay(fs.saturation(FluidSystem::oilPhaseIdx)) + * decay(fs.invB(FluidSystem::oilPhaseIdx)) + * decay(fs.Rs()) + * decay(intQuants.porosity()); + } + else { + phaseVolume = 0.0; + } - // evaluate the flux(es) over one face - void computeFlux_(TracerEvaluation & freeFlux, - bool & isUpFree, - const int tracerPhaseIdx, - const ElementContext& elemCtx, - unsigned scvfIdx, - unsigned timeIdx) + return max(phaseVolume, 1e-10); + } + + void computeFreeFlux_(TracerEvaluation & freeFlux, + bool & isUp, + const int tracerPhaseIdx, + const ElementContext& elemCtx, + unsigned scvfIdx, + unsigned timeIdx) { const auto& stencil = elemCtx.stencil(timeIdx); const auto& scvf = stencil.interiorFace(scvfIdx); @@ -265,17 +291,72 @@ class TracerModel : public GenericTracerModel(extQuants.volumeFlux(tracerPhaseIdx)) + * decay(fs.invB(tracerPhaseIdx)); + Scalar A = scvf.area(); - Scalar v = decay(extQuants.volumeFlux(tracerPhaseIdx)); - Scalar b = decay(fs.invB(tracerPhaseIdx)); + if (inIdx == upIdx) { + freeFlux = A*v*variable(1.0, 0); + isUp = true; + } + else { + freeFlux = A*v; + isUp = false; + } + } + void computeSolFlux_(TracerEvaluation & solFlux, + bool & isUp, + const int tracerPhaseIdx, + const ElementContext& elemCtx, + unsigned scvfIdx, + unsigned timeIdx) + { + const auto& stencil = elemCtx.stencil(timeIdx); + const auto& scvf = stencil.interiorFace(scvfIdx); + + const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx); + unsigned inIdx = extQuants.interiorIndex(); + + Scalar v; + unsigned upIdx; + + // vaporized oil + if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) { + upIdx = extQuants.upstreamIndex(FluidSystem::gasPhaseIdx); + + const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx); + const auto& fs = intQuants.fluidState(); + v = + decay(fs.invB(FluidSystem::gasPhaseIdx)) + * decay(extQuants.volumeFlux(FluidSystem::gasPhaseIdx)) + * decay(fs.Rv()); + } + // dissolved gas + else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) { + upIdx = extQuants.upstreamIndex(FluidSystem::oilPhaseIdx); + + const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx); + const auto& fs = intQuants.fluidState(); + v = + decay(fs.invB(FluidSystem::oilPhaseIdx)) + * decay(extQuants.volumeFlux(FluidSystem::oilPhaseIdx)) + * decay(fs.Rs()); + } + else { + upIdx = 0; + v = 0.0; + } + + Scalar A = scvf.area(); if (inIdx == upIdx) { - freeFlux = A*v*b*variable(1.0, 0); - isUpFree = true; + solFlux = A*v*variable(1.0, 0); + isUp = true; } else { - freeFlux = A*v*b*1.0; - isUpFree = false; + solFlux = A*v; + isUp = false; } } @@ -290,29 +371,46 @@ class TracerModel : public GenericTracerModel storageOfTimeIndex1(tr.numTracer()); + std::vector fStorageOfTimeIndex1(tr.numTracer()); + std::vector sStorageOfTimeIndex1(tr.numTracer()); if (elemCtx.enableStorageCache()) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { - storageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I]; + fStorageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I][0]; + sStorageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I][1]; } } else { - Scalar fVolume1; - computeVolume_(fVolume1, tr.phaseIdx_, elemCtx, 0, /*timeIdx=*/1); + Scalar fVolume1 = computeFreeVolume_(tr.phaseIdx_, I1, 1); + Scalar sVolume1 = computeSolutionVolume_(tr.phaseIdx_, I1, 1); for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { - storageOfTimeIndex1[tIdx] = fVolume1*tr.concentrationInitial_[tIdx][I1]; + fStorageOfTimeIndex1[tIdx] = fVolume1 * tr.concentration_[tIdx][I1][0]; + sStorageOfTimeIndex1[tIdx] = sVolume1 * tr.concentration_[tIdx][I1][1]; } } - TracerEvaluation fVolume; - computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, /*timeIdx=*/0); + TracerEvaluation fVol = computeFreeVolume_(tr.phaseIdx_, I, 0) * variable(1.0, 0); + TracerEvaluation sVol = computeSolutionVolume_(tr.phaseIdx_, I, 0) * variable(1.0, 0); + dsVol_[tr.phaseIdx_][I] += sVol.value() * scvVolume - sVol1_[tr.phaseIdx_][I]; + dfVol_[tr.phaseIdx_][I] += fVol.value() * scvVolume - fVol1_[tr.phaseIdx_][I]; for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { - Scalar storageOfTimeIndex0 = fVolume.value()*tr.concentration_[tIdx][I]; - Scalar localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1[tIdx]) * scvVolume/dt; - tr.residual_[tIdx][I][0] += localStorage; //residual + flux + // Free part + Scalar fStorageOfTimeIndex0 = fVol.value() * tr.concentration_[tIdx][I][0]; + Scalar fLocalStorage = (fStorageOfTimeIndex0 - fStorageOfTimeIndex1[tIdx]) * scvVolume/dt; + tr.residual_[tIdx][I][0] += fLocalStorage; //residual + flux + + // Solution part + Scalar sStorageOfTimeIndex0 = sVol.value() * tr.concentration_[tIdx][I][1]; + Scalar sLocalStorage = (sStorageOfTimeIndex0 - sStorageOfTimeIndex1[tIdx]) * scvVolume/dt; + tr.residual_[tIdx][I][1] += sLocalStorage; //residual + flux } - (*tr.mat)[I][I][0][0] += fVolume.derivative(0) * scvVolume/dt; + + // Derivative matrix + (*tr.mat)[I][I][0][0] += fVol.derivative(0) * scvVolume/dt; + (*tr.mat)[I][I][1][1] += sVol.derivative(0) * scvVolume/dt; + } template @@ -320,21 +418,36 @@ class TracerModel : public GenericTracerModelwellTracerRate_[std::make_pair(eclWell.name(), this->name(tr.idx_[tIdx]))] = 0.0; + this->wellFreeTracerRate_[std::make_pair(eclWell.name(), this->wellfname(tr.idx_[tIdx]))] = 0.0; + this->wellSolTracerRate_[std::make_pair(eclWell.name(), this->wellsname(tr.idx_[tIdx]))] = 0.0; + if (eclWell.isMultiSegment()) { + for (std::size_t i = 0; i < eclWell.getConnections().size(); ++i) { + this->mSwTracerRate_[std::make_tuple(eclWell.name(), + this->name(tr.idx_[tIdx]), + eclWell.getConnections().get(i).segment())] = 0.0; + } + } } - std::vector wtracer(tr.numTracer()); + std::vector wtracer(tr.numTracer()); for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { wtracer[tIdx] = this->currentConcentration_(eclWell, this->name(tr.idx_[tIdx])); } - for (auto& perfData : well.perforationData()) { - const auto I = perfData.cell_index; + Scalar dt = simulator_.timeStepSize(); + std::size_t well_index = simulator_.problem().wellModel().wellState().index(well.name()).value(); + const auto& ws = simulator_.problem().wellModel().wellState().well(well_index); + for (std::size_t i = 0; i < ws.perf_data.size(); ++i) { + const auto I = ws.perf_data.cell_index[i]; Scalar rate = well.volumetricSurfaceRateForConnection(I, tr.phaseIdx_); - if (rate > 0) { + Scalar rate_s; + if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) { + rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil]; + } + else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) { + rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas]; + } + else { + rate_s = 0.0; + } + + Scalar rate_f = rate - rate_s; + if (rate_f > 0) { 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(eclWell.name(),this->name(tr.idx_[tIdx]))) += rate*wtracer[tIdx]; + // Injection of free tracer only + tr.residual_[tIdx][I][0] -= rate_f*wtracer[tIdx]; + + // Store _injector_ tracer rate for reporting (can be done here since WTRACER is constant) + this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) += rate_f*wtracer[tIdx]; + this->wellFreeTracerRate_.at(std::make_pair(eclWell.name(),this->wellfname(tr.idx_[tIdx]))) += rate_f*wtracer[tIdx]; + if (eclWell.isMultiSegment()) { + this->mSwTracerRate_[std::make_tuple(eclWell.name(), + this->name(tr.idx_[tIdx]), + eclWell.getConnections().get(i).segment())] += rate_f*wtracer[tIdx]; + } + } + dfVol_[tr.phaseIdx_][I] -= rate_f * dt; + } + else if (rate_f < 0) { + for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { + // Production of free tracer + tr.residual_[tIdx][I][0] -= rate_f * tr.concentration_[tIdx][I][0]; + } + dfVol_[tr.phaseIdx_][I] -= rate_f * dt; + + // Derivative matrix for free tracer producer + (*tr.mat)[I][I][0][0] -= rate_f * variable(1.0, 0).derivative(0); } - else if (rate < 0) { + if (rate_s < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { - tr.residual_[tIdx][I][0] -= rate*tr.concentration_[tIdx][I]; + // Production of solution tracer + tr.residual_[tIdx][I][1] -= rate_s * tr.concentration_[tIdx][I][1]; } - (*tr.mat)[I][I][0][0] -= rate*variable(1.0, 0).derivative(0); + dsVol_[tr.phaseIdx_][I] -= rate_s * dt; + + // Derivative matrix for solution tracer producer + (*tr.mat)[I][I][1][1] -= rate_s * variable(1.0, 0).derivative(0); } } } + + template + void assembleTracerEquationSource(TrRe& tr, + const Scalar dt, + unsigned I) + { + if (tr.numTracer() == 0) + return; + + // Skip if solution tracers do not exist + if (tr.phaseIdx_ == FluidSystem::waterPhaseIdx || + (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && !FluidSystem::enableDissolvedGas()) || + (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && !FluidSystem::enableVaporizedOil())) { + return; + } + + const Scalar& dsVol = dsVol_[tr.phaseIdx_][I]; + const Scalar& dfVol = dfVol_[tr.phaseIdx_][I]; + + // Source term determined by sign of dsVol: if dsVol > 0 then ms -> mf, else mf -> ms + for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { + if (dsVol >= 0) { + tr.residual_[tIdx][I][0] -= (dfVol / dt) * tr.concentration_[tIdx][I][0]; + tr.residual_[tIdx][I][1] += (dfVol / dt) * tr.concentration_[tIdx][I][0]; + } + else { + tr.residual_[tIdx][I][0] += (dsVol / dt) * tr.concentration_[tIdx][I][1]; + tr.residual_[tIdx][I][1] -= (dsVol / dt) * tr.concentration_[tIdx][I][1]; + } + } + + // Derivative matrix + if (dsVol >= 0) { + (*tr.mat)[I][I][0][0] -= (dfVol / dt) * variable(1.0, 0).derivative(0); + (*tr.mat)[I][I][1][0] += (dfVol / dt) * variable(1.0, 0).derivative(0); + } + else { + (*tr.mat)[I][I][0][1] += (dsVol / dt) * variable(1.0, 0).derivative(0); + (*tr.mat)[I][I][1][1] -= (dsVol / dt) * variable(1.0, 0).derivative(0); + } + } void assembleTracerEquations_() { @@ -385,8 +589,17 @@ class TracerModel : public GenericTracerModelassembleTracerEquationWell(tr, *wellPtr); } } @@ -400,8 +613,10 @@ class TracerModel : public GenericTracerModelassembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J); + this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J, dt); } } - } - - // Well terms - const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells(); - for (const auto& wellPtr : wellPtrs) { + + // Source terms (mass transfer between free and solution tracer) for (auto& tr : tbatch) { - this->assembleTracerEquationWell(tr, *wellPtr); + this->assembleTracerEquationSource(tr, dt, I); } + } // Communicate overlap using grid Communication @@ -457,22 +670,31 @@ class TracerModel : public GenericTracerModel dx(tr.concentration_); - for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) + for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { dx[tIdx] = 0.0; + } bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_); if (!converged) { @@ -498,30 +721,99 @@ class TracerModel : public GenericTracerModeltracerConcentration_[tr.idx_[tIdx]] = tr.concentration_[tIdx]; + for (std::size_t globalDofIdx = 0; globalDofIdx < tr.concentration_[tIdx].size(); ++globalDofIdx) { + // New concetration. Concentrations that are negative or where free/solution phase is not + // present are set to zero + const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, 0); + const auto& fs = intQuants.fluidState(); + Scalar Sf = decay(fs.saturation(tr.phaseIdx_)); + Scalar Ss = 0.0; + if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) + Ss = decay(fs.saturation(FluidSystem::oilPhaseIdx)); + else if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) + Ss = decay(fs.saturation(FluidSystem::gasPhaseIdx)); + + const Scalar tol_gas_sat = 1e-6; + if ((tr.concentration_[tIdx][globalDofIdx][0] - dx[tIdx][globalDofIdx][0] < 0.0) + || Sf < tol_gas_sat) + tr.concentration_[tIdx][globalDofIdx][0] = 0.0; + else + tr.concentration_[tIdx][globalDofIdx][0] -= dx[tIdx][globalDofIdx][0]; + if (tr.concentration_[tIdx][globalDofIdx][1] - dx[tIdx][globalDofIdx][1] < 0.0 + || Ss < tol_gas_sat) + tr.concentration_[tIdx][globalDofIdx][1] = 0.0; + else + tr.concentration_[tIdx][globalDofIdx][1] -= dx[tIdx][globalDofIdx][1]; + + // Partition concentration into free and solution tracers for output + this->freeTracerConcentration_[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][0]; + this->solTracerConcentration_[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][1]; + } } // Store _producer_ tracer rate for reporting const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells(); for (const auto& wellPtr : wellPtrs) { - const auto& well = wellPtr->wellEcl(); + const auto& eclWell = wellPtr->wellEcl(); - if (!well.isProducer()) //Injection rates already reported during assembly + // Injection rates already reported during assembly + if (!eclWell.isProducer()) continue; Scalar rateWellPos = 0.0; Scalar rateWellNeg = 0.0; - for (auto& perfData : wellPtr->perforationData()) { - const int I = perfData.cell_index; + std::size_t well_index = simulator_.problem().wellModel().wellState().index(eclWell.name()).value(); + const auto& ws = simulator_.problem().wellModel().wellState().well(well_index); + for (std::size_t i = 0; i < ws.perf_data.size(); ++i) { + const auto I = ws.perf_data.cell_index[i]; Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_); - if (rate < 0) { - rateWellNeg += rate; + + Scalar rate_s; + if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) { + rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil]; + } + else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) { + rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas]; + } + else { + rate_s = 0.0; + } + + Scalar rate_f = rate - rate_s; + if (rate_f < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { - this->wellTracerRate_.at(std::make_pair(well.name(),this->name(tr.idx_[tIdx]))) += rate*tr.concentration_[tIdx][I]; + // Store _producer_ free tracer rate for reporting + this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) += + rate_f * tr.concentration_[tIdx][I][0]; + this->wellFreeTracerRate_.at(std::make_pair(eclWell.name(),this->wellfname(tr.idx_[tIdx]))) += + rate_f * tr.concentration_[tIdx][I][0]; + if (eclWell.isMultiSegment()) { + this->mSwTracerRate_[std::make_tuple(eclWell.name(), + this->name(tr.idx_[tIdx]), + eclWell.getConnections().get(i).segment())] += + rate_f * tr.concentration_[tIdx][I][0]; + } } } + if (rate_s < 0) { + for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { + // Store _producer_ solution tracer rate for reporting + this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) += + rate_s * tr.concentration_[tIdx][I][1]; + this->wellSolTracerRate_.at(std::make_pair(eclWell.name(),this->wellsname(tr.idx_[tIdx]))) += + rate_s * tr.concentration_[tIdx][I][1]; + if (eclWell.isMultiSegment()) { + this->mSwTracerRate_[std::make_tuple(eclWell.name(), + this->name(tr.idx_[tIdx]), + eclWell.getConnections().get(i).segment())] += + rate_s * tr.concentration_[tIdx][I][1]; + } + } + } + + if (rate < 0) { + rateWellNeg += rate; + } else { rateWellPos += rate; } @@ -532,7 +824,6 @@ class TracerModel : public GenericTracerModelwellTracerRate_.at(std::make_pair(well.name(),this->name(tr.idx_[tIdx]))) *= factor; + this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) *= factor; } } } } } + Simulator& simulator_; @@ -561,13 +853,13 @@ class TracerModel : public GenericTracerModel struct TracerBatch { - std::vector idx_; - const int phaseIdx_; - std::vector concentrationInitial_; - std::vector concentration_; - std::vector storageOfTimeIndex1_; - std::vector residual_; - std::unique_ptr mat; + std::vector idx_; + const int phaseIdx_; + std::vector concentrationInitial_; + std::vector concentration_; + std::vector storageOfTimeIndex1_; + std::vector residual_; + std::unique_ptr mat; bool operator==(const TracerBatch& rhs) const { @@ -600,12 +892,12 @@ class TracerModel : public GenericTracerModel& wat_; TracerBatch& oil_; TracerBatch& gas_; + std::array, 3> fVol1_; + std::array, 3> sVol1_; + std::array, 3> dsVol_; + std::array, 3> dfVol_; }; } // namespace Opm diff --git a/opm/simulators/flow/VtkTracerModule.hpp b/opm/simulators/flow/VtkTracerModule.hpp index 3d45089e0b6..c9df2b4fdeb 100644 --- a/opm/simulators/flow/VtkTracerModule.hpp +++ b/opm/simulators/flow/VtkTracerModule.hpp @@ -107,10 +107,14 @@ namespace Opm { { if (eclTracerConcentrationOutput_()){ const auto& tracerModel = this->simulator_.problem().tracerModel(); - eclTracerConcentration_.resize(tracerModel.numTracers()); - for (std::size_t tracerIdx = 0; tracerIdx < eclTracerConcentration_.size(); ++tracerIdx) { - - this->resizeScalarBuffer_(eclTracerConcentration_[tracerIdx]); + eclFreeTracerConcentration_.resize(tracerModel.numTracers()); + eclSolTracerConcentration_.resize(tracerModel.numTracers()); + const auto& enableSolTracers = tracerModel.enableSolTracers(); + + for (std::size_t tracerIdx = 0; tracerIdx < eclFreeTracerConcentration_.size(); ++tracerIdx) { + this->resizeScalarBuffer_(eclFreeTracerConcentration_[tracerIdx]); + if (enableSolTracers[tracerIdx]) + this->resizeScalarBuffer_(eclSolTracerConcentration_[tracerIdx]); } } @@ -125,14 +129,22 @@ namespace Opm { if (!Parameters::get()) return; - const auto& tracerModel = elemCtx.problem().tracerModel(); - - for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { - unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + if (eclTracerConcentrationOutput_()) { + const auto& tracerModel = elemCtx.problem().tracerModel(); + const auto& enableSolTracers = tracerModel.enableSolTracers(); - if (eclTracerConcentrationOutput_()){ - for (std::size_t tracerIdx = 0; tracerIdx < eclTracerConcentration_.size(); ++tracerIdx) { - eclTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.tracerConcentration(tracerIdx, globalDofIdx); + for (std::size_t tracerIdx = 0; tracerIdx < eclFreeTracerConcentration_.size(); ++tracerIdx) { + // free tracer + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + eclFreeTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx); + } + // solution tracer (only if it exist) + if (enableSolTracers[tracerIdx]) { + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + eclSolTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.solTracerConcentration(tracerIdx, globalDofIdx); + } } } } @@ -149,9 +161,15 @@ namespace Opm { if (eclTracerConcentrationOutput_()){ const auto& tracerModel = this->simulator_.problem().tracerModel(); - for (std::size_t tracerIdx = 0; tracerIdx < eclTracerConcentration_.size(); ++tracerIdx) { - const std::string tmp = "tracerConcentration_" + tracerModel.name(tracerIdx); - this->commitScalarBuffer_(baseWriter,tmp.c_str(), eclTracerConcentration_[tracerIdx]); + const auto& enableSolTracers = tracerModel.enableSolTracers(); + + for (std::size_t tracerIdx = 0; tracerIdx < eclFreeTracerConcentration_.size(); ++tracerIdx) { + const std::string tmp = "freeTracerConcentration_" + tracerModel.name(tracerIdx); + this->commitScalarBuffer_(baseWriter, tmp.c_str(), eclFreeTracerConcentration_[tracerIdx]); + if (enableSolTracers[tracerIdx]) { + const std::string tmp2 = "solTracerConcentration_" + tracerModel.name(tracerIdx); + this->commitScalarBuffer_(baseWriter, tmp2.c_str(), eclSolTracerConcentration_[tracerIdx]); + } } } @@ -167,7 +185,8 @@ namespace Opm { } - std::vector eclTracerConcentration_; + std::vector eclFreeTracerConcentration_; + std::vector eclSolTracerConcentration_; }; } // namespace Opm diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 2b4cb1fd954..0fe1ef12148 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -267,8 +267,17 @@ template class WellContributions; { const auto& tracerRates = this->simulator_.problem() .tracerModel().getWellTracerRates(); + const auto& freeTracerRates = simulator_.problem() + .tracerModel().getWellFreeTracerRates(); + const auto& solTracerRates = simulator_.problem() + .tracerModel().getWellSolTracerRates(); + const auto& mswTracerRates = simulator_.problem() + .tracerModel().getMswTracerRates(); this->assignWellTracerRates(wsrpt, tracerRates); + this->assignWellTracerRates(wsrpt, freeTracerRates); + this->assignWellTracerRates(wsrpt, solTracerRates); + this->assignMswTracerRates(wsrpt, mswTracerRates); } BlackoilWellModelGuideRates(*this) diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 3fae38c5476..eb4dc2d1160 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1698,6 +1698,33 @@ assignWellTracerRates(data::Wells& wsrpt, } } +template +void BlackoilWellModelGeneric:: +assignMswTracerRates(data::Wells& wsrpt, + const MswTracerRates& mswTracerRates) const +{ + if (mswTracerRates.empty()) + return; + + for (const auto& mswTR : mswTracerRates) { + std::string wellName = std::get<0>(mswTR.first); + auto xwPos = wsrpt.find(wellName); + if (xwPos == wsrpt.end()) { // No well results. + continue; + } + std::string tracerName = std::get<1>(mswTR.first); + std::size_t segNumber = std::get<2>(mswTR.first); + Scalar rate = mswTR.second; + + auto& wData = xwPos->second; + auto segPos = wData.segments.find(segNumber); + if (segPos != wData.segments.end()) { + auto& segment = segPos->second; + segment.rates.set(data::Rates::opt::tracer, rate, tracerName); + } + } +} + template std::vector> BlackoilWellModelGeneric:: getMaxWellConnections() const diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp index 1fb7f09a341..98f512622aa 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -437,6 +437,9 @@ class BlackoilWellModelGeneric using WellTracerRates = std::map, Scalar>; void assignWellTracerRates(data::Wells& wsrpt, const WellTracerRates& wellTracerRates) const; + using MswTracerRates = std::map, Scalar>; + void assignMswTracerRates(data::Wells& wsrpt, + const MswTracerRates& mswTracerRates) const; Schedule& schedule_; const SummaryState& summaryState_; diff --git a/tests/test_RestartSerialization.cpp b/tests/test_RestartSerialization.cpp index 1cc23282f98..6c71c70d81a 100644 --- a/tests/test_RestartSerialization.cpp +++ b/tests/test_RestartSerialization.cpp @@ -381,10 +381,10 @@ BOOST_AUTO_TEST_CASE(BlackoilWellModelGeneric) BOOST_CHECK_MESSAGE(data_out == data_in, "Deserialized BlackoilWellModelGeneric differ"); } -template -class GenericTracerModelTest : public Opm::GenericTracerModel +template +class GenericTracerModelTest : public Opm::GenericTracerModel { - using Base = Opm::GenericTracerModel; + using Base = Opm::GenericTracerModel; public: GenericTracerModelTest(const GridView& gridView, const Opm::EclipseState& eclState, @@ -402,7 +402,7 @@ class GenericTracerModelTest : public Opm::GenericTracerModel(int)> centroids) { GenericTracerModelTest result(gridView, eclState, cartMapper, dofMapper, centroids); - result.tracerConcentration_ = {{1.0}, {2.0}, {3.0}}; + result.tracerConcentration_ = {{{1.0, 2.0}}, {{3.0, 4.0}}, {{5.0, 6.0}}}; result.wellTracerRate_.insert({{"foo", "bar"}, 4.0}); return result; @@ -438,6 +438,7 @@ BOOST_AUTO_TEST_CASE(FlowGenericTracerModel) GridView, Dune::MultipleCodimMultipleGeomTypeMapper, Opm::EcfvStencil, + Opm::BlackOilFluidSystem, double> ::serializationTestObject(gridView, eclState, mapper, dofMapper, centroids); Opm::Serialization::MemPacker packer; @@ -474,6 +475,7 @@ BOOST_AUTO_TEST_CASE(FlowGenericTracerModelFem) GridView, Dune::MultipleCodimMultipleGeomTypeMapper, Opm::EcfvStencil, + Opm::BlackOilFluidSystem, double> ::serializationTestObject(gridView, eclState, mapper, dofMapper, centroids); Opm::Serialization::MemPacker packer;