diff --git a/opm/simulators/flow/EclWriter.hpp b/opm/simulators/flow/EclWriter.hpp index 92e01a6b6ff..6138b6d1863 100644 --- a/opm/simulators/flow/EclWriter.hpp +++ b/opm/simulators/flow/EclWriter.hpp @@ -501,8 +501,10 @@ class EclWriter : public EclGenericWriter { const auto& tracers = simulator_.vanguard().eclState().tracer(); 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, true); + solutionKeys.emplace_back(tracer.sname(), UnitSystem::measure::identity, enableSolTracer); } } @@ -529,14 +531,30 @@ class EclWriter : public EclGenericWriter auto& tracer_model = simulator_.problem().tracerModel(); for (int tracer_index = 0; tracer_index < tracer_model.numTracers(); tracer_index++) { + // Free tracers const auto& free_tracer_name = tracer_model.fname(tracer_index); const auto& free_tracer_solution = restartValues.solution.template data(free_tracer_name); - 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.setFreeTracerConcentration(tracer_index, globalIdx, free_tracer_solution[globalIdx]); - tracer_model.setSolTracerConcentration(tracer_index, globalIdx, sol_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 06c2069d3d2..1dcfd599f0d 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -721,6 +721,9 @@ assignToSolution(data::Solution& sol) 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]), @@ -854,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. @@ -1328,7 +1332,8 @@ doAllocBuffers(const unsigned bufferSize, solTracerConcentrations_.resize(numTracers); for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) { - solTracerConcentrations_[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 0575f920a42..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, 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 32d9a5c8544..312f7b06eb7 100644 --- a/opm/simulators/flow/GenericTracerModel.hpp +++ b/opm/simulators/flow/GenericTracerModel.hpp @@ -36,6 +36,8 @@ #include +#include + #include #include #include @@ -49,7 +51,7 @@ namespace Opm { class EclipseState; class Well; -template +template class GenericTracerModel { public: using TracerVectorSingle = Dune::BlockVector>; @@ -71,6 +73,8 @@ class GenericTracerModel { 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 @@ -79,6 +83,7 @@ class GenericTracerModel { 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 @@ -96,9 +101,12 @@ class GenericTracerModel { void serializeOp(Serializer& serializer) { serializer(tracerConcentration_); + serializer(freeTracerConcentration_); + serializer(solTracerConcentration_); serializer(wellTracerRate_); serializer(wellFreeTracerRate_); serializer(wellSolTracerRate_); + serializer(mSwTracerRate_); } protected: @@ -131,6 +139,7 @@ class GenericTracerModel { const DofMapper& dofMapper_; std::vector tracerPhaseIdx_; + std::vector enableSolTracers_; std::vector tracerConcentration_; std::unique_ptr tracerMatrix_; std::vector freeTracerConcentration_; diff --git a/opm/simulators/flow/GenericTracerModel_impl.hpp b/opm/simulators/flow/GenericTracerModel_impl.hpp index 5f646526a46..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,8 +111,8 @@ GenericTracerModel(const GridView& gridView, { } -template -Scalar GenericTracerModel:: +template +Scalar GenericTracerModel:: freeTracerConcentration(int tracerIdx, int globalDofIdx) const { if (freeTracerConcentration_.empty()) @@ -122,8 +121,8 @@ freeTracerConcentration(int tracerIdx, int globalDofIdx) const return freeTracerConcentration_[tracerIdx][globalDofIdx]; } -template -Scalar GenericTracerModel:: +template +Scalar GenericTracerModel:: solTracerConcentration(int tracerIdx, int globalDofIdx) const { if (solTracerConcentration_.empty()) @@ -132,71 +131,94 @@ solTracerConcentration(int tracerIdx, int globalDofIdx) const return solTracerConcentration_[tracerIdx][globalDofIdx]; } -template -void GenericTracerModel:: +template +void GenericTracerModel:: setFreeTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value) { this->freeTracerConcentration_[tracerIdx][globalDofIdx] = value; + this->tracerConcentration_[tracerIdx][globalDofIdx][0] = value; } -template -void GenericTracerModel:: +template +void GenericTracerModel:: setSolTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value) { this->solTracerConcentration_[tracerIdx][globalDofIdx] = value; + this->tracerConcentration_[tracerIdx][globalDofIdx][1] = value; } -template -int GenericTracerModel:: +template +void GenericTracerModel:: +setEnableSolTracers(int tracerIdx, bool enableSolTracer) +{ + this->enableSolTracers_[tracerIdx] = enableSolTracer; +} + +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 -std::string GenericTracerModel:: +template +std::string GenericTracerModel:: sname(int tracerIdx) const { return this->eclState_.tracer()[tracerIdx].sname(); } -template -std::string GenericTracerModel:: +template +std::string GenericTracerModel:: wellfname(int tracerIdx) const { return this->eclState_.tracer()[tracerIdx].wellfname(); } -template -std::string GenericTracerModel:: +template +std::string GenericTracerModel:: wellsname(int tracerIdx) const { return this->eclState_.tracer()[tracerIdx].wellsname(); } -template -Scalar GenericTracerModel:: +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) { @@ -207,6 +229,7 @@ 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); freeTracerConcentration_.resize(numTracers); solTracerConcentration_.resize(numTracers); @@ -223,7 +246,6 @@ doInit(bool rst, std::size_t numGridDof, else if (tracer.phase == Phase::GAS) tracerPhaseIdx_[tracerIdx] = gasPhaseIdx; - tracerConcentration_[tracerIdx].resize(numGridDof); tracerConcentration_[tracerIdx].resize(numGridDof); freeTracerConcentration_[tracerIdx].resize(numGridDof); solTracerConcentration_[tracerIdx].resize(numGridDof); @@ -258,36 +280,63 @@ doInit(bool rst, std::size_t numGridDof, } } else { - throw std::logic_error(fmt::format("Can not initialize free tracer concentration: {}", tracer.name)); + 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; + } } - // TBLKS keyword - if (tracer.solution_concentration.has_value()){ - 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); + // 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]; + } } - 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; + } } } - // TVDPS keyword - else if (tracer.solution_tvdp.has_value()) { - const auto& solution_tvdp = tracer.solution_tvdp.value(); + else { + // No solution tracers, default to zero + enableSolTracers_[tracerIdx] = false; 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]); + tracerConcentration_[tracerIdx][globalDofIdx][1] = 0.0; + solTracerConcentration_[tracerIdx][globalDofIdx] = 0.0; } - } - else { - throw std::logic_error(fmt::format("Can not initialize solution tracer concentration: {}", tracer.name)); } } @@ -331,8 +380,8 @@ doInit(bool rst, std::size_t numGridDof, tracerMatrix_->endindices(); } -template -bool GenericTracerModel:: +template +bool GenericTracerModel:: linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b) { x = 0.0; @@ -386,11 +435,11 @@ 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-6; + Scalar tolerance = 1e-2; int maxIter = 100; int verbosity = 0; diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 976616b6673..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()); } diff --git a/opm/simulators/flow/TracerModel.hpp b/opm/simulators/flow/TracerModel.hpp index 11a062e4774..30c414f778c 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; @@ -494,7 +496,7 @@ class TracerModel : public GenericTracerModelwellTracerRate_.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()) { @@ -506,39 +508,23 @@ class TracerModel : public GenericTracerModelwellTracerRate_.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]; - } } dfVol_[tr.phaseIdx_][I] -= rate_f * dt; - // Derivative matrix for producer + // Derivative matrix for free tracer producer (*tr.mat)[I][I][0][0] -= rate_f * variable(1.0, 0).derivative(0); } if (rate_s < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { + // Production of solution tracer tr.residual_[tIdx][I][1] -= rate_s * tr.concentration_[tIdx][I][1]; - - 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]; - } } 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); } } @@ -552,6 +538,13 @@ class TracerModel : public GenericTracerModelwellEcl(); + 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_); + + 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) { + // 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; } @@ -776,7 +812,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; } } } diff --git a/opm/simulators/flow/VtkTracerModule.hpp b/opm/simulators/flow/VtkTracerModule.hpp index 018ca785c09..c9df2b4fdeb 100644 --- a/opm/simulators/flow/VtkTracerModule.hpp +++ b/opm/simulators/flow/VtkTracerModule.hpp @@ -109,10 +109,12 @@ namespace Opm { const auto& tracerModel = this->simulator_.problem().tracerModel(); eclFreeTracerConcentration_.resize(tracerModel.numTracers()); eclSolTracerConcentration_.resize(tracerModel.numTracers()); - for (std::size_t tracerIdx = 0; tracerIdx < eclFreeTracerConcentration_.size(); ++tracerIdx) { + const auto& enableSolTracers = tracerModel.enableSolTracers(); + for (std::size_t tracerIdx = 0; tracerIdx < eclFreeTracerConcentration_.size(); ++tracerIdx) { this->resizeScalarBuffer_(eclFreeTracerConcentration_[tracerIdx]); - this->resizeScalarBuffer_(eclSolTracerConcentration_[tracerIdx]); + if (enableSolTracers[tracerIdx]) + this->resizeScalarBuffer_(eclSolTracerConcentration_[tracerIdx]); } } @@ -127,15 +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 < eclFreeTracerConcentration_.size(); ++tracerIdx) { + 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); - eclSolTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.solTracerConcentration(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); + } } } } @@ -152,11 +161,15 @@ namespace Opm { if (eclTracerConcentrationOutput_()){ const auto& tracerModel = this->simulator_.problem().tracerModel(); + 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]); - const std::string tmp2 = "solTracerConcentration_" + tracerModel.name(tracerIdx); - this->commitScalarBuffer_(baseWriter,tmp2.c_str(), eclSolTracerConcentration_[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]); + } } } diff --git a/tests/test_RestartSerialization.cpp b/tests/test_RestartSerialization.cpp index a53aa3f7911..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, @@ -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;