OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5268 from svenn-t/tracer_dis_vap

Browse Source 
Partitioning tracers
This commit is contained in:
Tor Harald Sandve
2024-06-14 09:08:12 +02:00
committed by GitHub
parent cc6acb7722 87362d5037
commit 151bba43a7
13 changed files with 765 additions and 185 deletions
Download Patch File Download Diff File Expand all files Collapse all files

31
opm/simulators/flow/EclWriter.hpp
View File

@@ -500,8 +500,12 @@ public:
{ {
const auto& tracers = simulator_.vanguard().eclState().tracer(); 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.fname(), UnitSystem::measure::identity, true);
solutionKeys.emplace_back(tracer.sname(), UnitSystem::measure::identity, enableSolTracer);
}
} }
// The episodeIndex is rewined one back before beginRestart is called // The episodeIndex is rewined one back before beginRestart is called
@@ -527,11 +531,30 @@ public:
auto& tracer_model = simulator_.problem().tracerModel(); auto& tracer_model = simulator_.problem().tracerModel();
for (int tracer_index = 0; tracer_index < tracer_model.numTracers(); tracer_index++) { for (int tracer_index = 0; tracer_index < tracer_model.numTracers(); tracer_index++) {
const auto& tracer_name = tracer_model.fname(tracer_index); // Free tracers
const auto& tracer_solution = restartValues.solution.template data<double>(tracer_name); const auto& free_tracer_name = tracer_model.fname(tracer_index);
const auto& free_tracer_solution = restartValues.solution.template data<double>(free_tracer_name);
for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
unsigned globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(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<double>(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);
}
} }
} }

38
opm/simulators/flow/GenericOutputBlackoilModule.cpp
View File

@@ -700,21 +700,40 @@ assignToSolution(data::Solution& sol)
} }
// Tracers // Tracers
if (! this->tracerConcentrations_.empty()) { if (! this->freeTracerConcentrations_.empty()) {
const auto& tracers = this->eclState_.tracer(); const auto& tracers = this->eclState_.tracer();
for (auto tracerIdx = 0*tracers.size(); for (auto tracerIdx = 0*tracers.size();
tracerIdx < tracers.size(); ++tracerIdx) tracerIdx < tracers.size(); ++tracerIdx)
{ {
sol.insert(tracers[tracerIdx].fname(), sol.insert(tracers[tracerIdx].fname(),
UnitSystem::measure::identity, UnitSystem::measure::identity,
std::move(tracerConcentrations_[tracerIdx]), std::move(freeTracerConcentrations_[tracerIdx]),
data::TargetType::RESTART_TRACER_SOLUTION); 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 // we'll move from vectors that have already been moved from if we
// get here and it's not a restart step. // 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 vapparsActive,
const bool enableHysteresis, const bool enableHysteresis,
const unsigned numTracers, const unsigned numTracers,
const std::vector<bool>& enableSolTracers,
const unsigned numOutputNnc) const unsigned numOutputNnc)
{ {
// Output RESTART_OPM_EXTENDED only when explicitly requested by user. // Output RESTART_OPM_EXTENDED only when explicitly requested by user.
@@ -1304,10 +1324,16 @@ doAllocBuffers(const unsigned bufferSize,
// tracers // tracers
if (numTracers > 0) { if (numTracers > 0) {
tracerConcentrations_.resize(numTracers); freeTracerConcentrations_.resize(numTracers);
for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx)
{ {
tracerConcentrations_[tracerIdx].resize(bufferSize, 0.0); freeTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
}
solTracerConcentrations_.resize(numTracers);
for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx)
{
if (enableSolTracers[tracerIdx])
solTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
} }
} }

4
opm/simulators/flow/GenericOutputBlackoilModule.hpp
View File

@@ -340,6 +340,7 @@ protected:
const bool vapparsActive, const bool vapparsActive,
const bool enableHysteresis, const bool enableHysteresis,
unsigned numTracers, unsigned numTracers,
const std::vector<bool>& enableSolTracers,
unsigned numOutputNnc); unsigned numOutputNnc);
void makeRegionSum(Inplace& inplace, void makeRegionSum(Inplace& inplace,
@@ -504,7 +505,8 @@ protected:
std::array<ScalarBuffer, numPhases> viscosity_; std::array<ScalarBuffer, numPhases> viscosity_;
std::array<ScalarBuffer, numPhases> relativePermeability_; std::array<ScalarBuffer, numPhases> relativePermeability_;
std::vector<ScalarBuffer> tracerConcentrations_; std::vector<ScalarBuffer> freeTracerConcentrations_;
std::vector<ScalarBuffer> solTracerConcentrations_;
std::array<ScalarBuffer, numPhases> residual_; std::array<ScalarBuffer, numPhases> residual_;

6
opm/simulators/flow/GenericTracerModel.cpp
View File

@@ -24,6 +24,8 @@
#include <opm/simulators/flow/GenericTracerModel_impl.hpp> #include <opm/simulators/flow/GenericTracerModel_impl.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#if HAVE_DUNE_FEM #if HAVE_DUNE_FEM
#include <dune/common/version.hh> #include <dune/common/version.hh>
#include <dune/fem/gridpart/adaptiveleafgridpart.hh> #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
@@ -39,6 +41,7 @@ template class GenericTracerModel<Dune::CpGrid,
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>, Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
Opm::EcfvStencil<double,Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,false,false>, Opm::EcfvStencil<double,Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,false,false>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>; double>;
#if HAVE_DUNE_FEM #if HAVE_DUNE_FEM
@@ -50,12 +53,14 @@ template class GenericTracerModel<Dune::CpGrid,
GV, GV,
Dune::MultipleCodimMultipleGeomTypeMapper<GV>, Dune::MultipleCodimMultipleGeomTypeMapper<GV>,
EcfvStencil<double, GV, false, false>, EcfvStencil<double, GV, false, false>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>; double>;
#else #else
template class GenericTracerModel<Dune::CpGrid, template class GenericTracerModel<Dune::CpGrid,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>, Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>>,
EcfvStencil<double,Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,false,false>, EcfvStencil<double,Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,false,false>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>; double>;
template class GenericTracerModel<Dune::CpGrid, template class GenericTracerModel<Dune::CpGrid,
Dune::Fem::GridPart2GridViewImpl<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, (Dune::PartitionIteratorType)4, false> >, Dune::Fem::GridPart2GridViewImpl<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, (Dune::PartitionIteratorType)4, false> >,
@@ -65,6 +70,7 @@ template class GenericTracerModel<Dune::CpGrid,
EcfvStencil<double, Dune::Fem::GridPart2GridViewImpl< EcfvStencil<double, Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false> >, Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false> >,
false, false>, false, false>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>; double>;
#endif #endif
#endif // HAVE_DUNE_FEM #endif // HAVE_DUNE_FEM

46
opm/simulators/flow/GenericTracerModel.hpp
View File

@@ -36,6 +36,8 @@
#include <opm/simulators/linalg/matrixblock.hh> #include <opm/simulators/linalg/matrixblock.hh>
#include <opm/input/eclipse/EclipseState/Phase.hpp>
#include <array> #include <array>
#include <cstddef> #include <cstddef>
#include <functional> #include <functional>
@@ -49,11 +51,12 @@ namespace Opm {
class EclipseState; class EclipseState;
class Well; class Well;
template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
class GenericTracerModel { class GenericTracerModel {
public: public:
using TracerMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 1, 1>>; using TracerVectorSingle = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
using TracerVector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>; using TracerMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 2, 2>>;
using TracerVector = Dune::BlockVector<Dune::FieldVector<Scalar, 2>>;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>; using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
static constexpr int dimWorld = Grid::dimensionworld; static constexpr int dimWorld = Grid::dimensionworld;
/*! /*!
@@ -66,25 +69,44 @@ public:
*/ */
const std::string& name(int tracerIdx) const; const std::string& name(int tracerIdx) const;
std::string fname(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<bool>& enableSolTracers() const;
/*! /*!
* \brief Return the tracer concentration for tracer index and global DofIdx * \brief Return the tracer concentration for tracer index and global DofIdx
*/ */
Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const; Scalar freeTracerConcentration(int tracerIdx, int globalDofIdx) const;
void setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value); 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 * \brief Return well tracer rates
*/ */
const std::map<std::pair<std::string, std::string>, Scalar>& const std::map<std::pair<std::string, std::string>, Scalar>&
getWellTracerRates() const {return wellTracerRate_;} getWellTracerRates() const {return wellTracerRate_;}
const std::map<std::pair<std::string, std::string>, Scalar>&
getWellFreeTracerRates() const {return wellFreeTracerRate_;}
const std::map<std::pair<std::string, std::string>, Scalar>&
getWellSolTracerRates() const {return wellSolTracerRate_;}
const std::map<std::tuple<std::string, std::string, std::size_t>, Scalar>&
getMswTracerRates() const {return mSwTracerRate_;}
template<class Serializer> template<class Serializer>
void serializeOp(Serializer& serializer) void serializeOp(Serializer& serializer)
{ {
serializer(tracerConcentration_); serializer(tracerConcentration_);
serializer(freeTracerConcentration_);
serializer(solTracerConcentration_);
serializer(wellTracerRate_); serializer(wellTracerRate_);
serializer(wellFreeTracerRate_);
serializer(wellSolTracerRate_);
serializer(mSwTracerRate_);
} }
protected: protected:
@@ -117,12 +139,20 @@ protected:
const DofMapper& dofMapper_; const DofMapper& dofMapper_;
std::vector<int> tracerPhaseIdx_; std::vector<int> tracerPhaseIdx_;
std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> tracerConcentration_; std::vector<bool> enableSolTracers_;
std::vector<TracerVector> tracerConcentration_;
std::unique_ptr<TracerMatrix> tracerMatrix_; std::unique_ptr<TracerMatrix> tracerMatrix_;
std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> storageOfTimeIndex1_; std::vector<TracerVectorSingle> freeTracerConcentration_;
std::vector<TracerVectorSingle> solTracerConcentration_;
// <wellName, tracerIdx> -> wellRate // <wellName, tracerName> -> wellRate
std::map<std::pair<std::string, std::string>, Scalar> wellTracerRate_; std::map<std::pair<std::string, std::string>, Scalar> wellTracerRate_;
std::map<std::pair<std::string, std::string>, Scalar> wellFreeTracerRate_;
std::map<std::pair<std::string, std::string>, Scalar> wellSolTracerRate_;
// <wellName, tracerName, segNum> -> wellRate
std::map<std::tuple<std::string, std::string, std::size_t>, Scalar> mSwTracerRate_;
/// \brief Function returning the cell centers /// \brief Function returning the cell centers
std::function<std::array<double,dimWorld>(int)> centroids_; std::function<std::array<double,dimWorld>(int)> centroids_;
}; };

201
opm/simulators/flow/GenericTracerModel_impl.hpp
View File

@@ -39,7 +39,6 @@
#include <opm/grid/CpGrid.hpp> #include <opm/grid/CpGrid.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp> #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Phase.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TracerVdTable.hpp> #include <opm/input/eclipse/EclipseState/Tables/TracerVdTable.hpp>
#include <opm/input/eclipse/Schedule/Well/Well.hpp> #include <opm/input/eclipse/Schedule/Well/Well.hpp>
#include <opm/input/eclipse/Schedule/Well/WellTracerProperties.hpp> #include <opm/input/eclipse/Schedule/Well/WellTracerProperties.hpp>
@@ -97,8 +96,8 @@ createParallelFlexibleSolver(const Dune::CpGrid& grid, const Matrix& M, const Pr
} }
#endif #endif
template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>:: GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
GenericTracerModel(const GridView& gridView, GenericTracerModel(const GridView& gridView,
const EclipseState& eclState, const EclipseState& eclState,
const CartesianIndexMapper& cartMapper, const CartesianIndexMapper& cartMapper,
@@ -112,53 +111,114 @@ GenericTracerModel(const GridView& gridView,
{ {
} }
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>:: Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
tracerConcentration(int tracerIdx, int globalDofIdx) const freeTracerConcentration(int tracerIdx, int globalDofIdx) const
{ {
if (tracerConcentration_.empty()) if (freeTracerConcentration_.empty())
return 0.0; return 0.0;
return tracerConcentration_[tracerIdx][globalDofIdx]; return freeTracerConcentration_[tracerIdx][globalDofIdx];
} }
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
void GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>:: Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value) solTracerConcentration(int tracerIdx, int globalDofIdx) const
{ {
this->tracerConcentration_[tracerIdx][globalDofIdx] = value; if (solTracerConcentration_.empty())
return 0.0;
return solTracerConcentration_[tracerIdx][globalDofIdx];
} }
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
int GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>:: void GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
setFreeTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value)
{
this->freeTracerConcentration_[tracerIdx][globalDofIdx] = value;
this->tracerConcentration_[tracerIdx][globalDofIdx][0] = value;
}
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
void GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
setSolTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value)
{
this->solTracerConcentration_[tracerIdx][globalDofIdx] = value;
this->tracerConcentration_[tracerIdx][globalDofIdx][1] = value;
}
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
void GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
setEnableSolTracers(int tracerIdx, bool enableSolTracer)
{
this->enableSolTracers_[tracerIdx] = enableSolTracer;
}
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
int GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
numTracers() const numTracers() const
{ {
return this->eclState_.tracer().size(); return this->eclState_.tracer().size();
} }
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>:: std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
fname(int tracerIdx) const fname(int tracerIdx) const
{ {
return this->eclState_.tracer()[tracerIdx].fname(); return this->eclState_.tracer()[tracerIdx].fname();
} }
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>:: std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
sname(int tracerIdx) const
{
return this->eclState_.tracer()[tracerIdx].sname();
}
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
wellfname(int tracerIdx) const
{
return this->eclState_.tracer()[tracerIdx].wellfname();
}
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
wellsname(int tracerIdx) const
{
return this->eclState_.tracer()[tracerIdx].wellsname();
}
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
Phase GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
phase(int tracerIdx) const
{
return this->eclState_.tracer()[tracerIdx].phase;
}
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
const std::vector<bool>& GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
enableSolTracers() const
{
return this->enableSolTracers_;
}
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
currentConcentration_(const Well& eclWell, const std::string& name) const currentConcentration_(const Well& eclWell, const std::string& name) const
{ {
return eclWell.getTracerProperties().getConcentration(name); return eclWell.getTracerProperties().getConcentration(name);
} }
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
const std::string& GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>:: const std::string& GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
name(int tracerIdx) const name(int tracerIdx) const
{ {
return this->eclState_.tracer()[tracerIdx].name; return this->eclState_.tracer()[tracerIdx].name;
} }
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
void GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>:: void GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
doInit(bool rst, std::size_t numGridDof, doInit(bool rst, std::size_t numGridDof,
std::size_t gasPhaseIdx, std::size_t oilPhaseIdx, std::size_t waterPhaseIdx) 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 // retrieve the number of tracers from the deck
const std::size_t numTracers = tracers.size(); const std::size_t numTracers = tracers.size();
enableSolTracers_.resize(numTracers);
tracerConcentration_.resize(numTracers); tracerConcentration_.resize(numTracers);
storageOfTimeIndex1_.resize(numTracers); freeTracerConcentration_.resize(numTracers);
solTracerConcentration_.resize(numTracers);
// the phase where the tracer is // the phase where the tracer is
tracerPhaseIdx_.resize(numTracers); tracerPhaseIdx_.resize(numTracers);
@@ -185,34 +247,97 @@ doInit(bool rst, std::size_t numGridDof,
tracerPhaseIdx_[tracerIdx] = gasPhaseIdx; tracerPhaseIdx_[tracerIdx] = gasPhaseIdx;
tracerConcentration_[tracerIdx].resize(numGridDof); tracerConcentration_[tracerIdx].resize(numGridDof);
storageOfTimeIndex1_[tracerIdx].resize(numGridDof); freeTracerConcentration_[tracerIdx].resize(numGridDof);
solTracerConcentration_[tracerIdx].resize(numGridDof);
if (rst) if (rst)
continue; continue;
//TBLK keyword // TBLKF keyword
if (tracer.free_concentration.has_value()){ if (tracer.free_concentration.has_value()){
const auto& free_concentration = tracer.free_concentration.value(); const auto& free_concentration = tracer.free_concentration.value();
int tblkDatasize = free_concentration.size(); int tblkDatasize = free_concentration.size();
if (tblkDatasize < cartMapper_.cartesianSize()){ 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) { for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
int cartDofIdx = cartMapper_.cartesianIndex(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()) { else if (tracer.free_tvdp.has_value()) {
const auto& free_tvdp = tracer.free_tvdp.value(); const auto& free_tvdp = tracer.free_tvdp.value();
for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) { 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", free_tvdp.evaluate("TRACER_CONCENTRATION",
centroids_(globalDofIdx)[2]); centroids_(globalDofIdx)[2]);
} }
} else }
throw std::logic_error(fmt::format("Can not initialize tracer: {}", tracer.name)); 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;
}
}
} }
// allocate matrix for storing the Jacobian of the tracer residual // 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 // 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_->setrowsize(dofIdx, neighbors[dofIdx].size());
}
tracerMatrix_->endrowsizes(); tracerMatrix_->endrowsizes();
// fill the rows with indices. each degree of freedom talks to // 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) { for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
typename NeighborSet::iterator nIt = neighbors[dofIdx].begin(); typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end(); typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
for (; nIt != nEndIt; ++nIt) for (; nIt != nEndIt; ++nIt) {
tracerMatrix_->addindex(dofIdx, *nIt); tracerMatrix_->addindex(dofIdx, *nIt);
}
} }
tracerMatrix_->endindices(); tracerMatrix_->endindices();
} }
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
bool GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>:: bool GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b) linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
{ {
x = 0.0; x = 0.0;
@@ -308,8 +435,8 @@ linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
#endif #endif
} }
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
bool GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>:: bool GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b) linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b)
{ {
Scalar tolerance = 1e-2; Scalar tolerance = 1e-2;

20
opm/simulators/flow/OutputBlackoilModule.hpp
View File

@@ -240,6 +240,7 @@ public:
problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)), problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)),
problem.materialLawManager()->enableHysteresis(), problem.materialLawManager()->enableHysteresis(),
problem.tracerModel().numTracers(), problem.tracerModel().numTracers(),
problem.tracerModel().enableSolTracers(),
problem.eclWriter()->getOutputNnc().size()); problem.eclWriter()->getOutputNnc().size());
} }
@@ -661,14 +662,23 @@ public:
// tracers // tracers
const auto& tracerModel = simulator_.problem().tracerModel(); const auto& tracerModel = simulator_.problem().tracerModel();
if (! this->tracerConcentrations_.empty()) { if (! this->freeTracerConcentrations_.empty()) {
for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) { for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
if (this->tracerConcentrations_[tracerIdx].empty()) { if (this->freeTracerConcentrations_[tracerIdx].empty()) {
continue; continue;
} }
this->freeTracerConcentrations_[tracerIdx][globalDofIdx] =
this->tracerConcentrations_[tracerIdx][globalDofIdx] = tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx);
tracerModel.tracerConcentration(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);
} }
} }

508
opm/simulators/flow/TracerModel.hpp
View File

@@ -63,12 +63,14 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
GetPropType<TypeTag, Properties::GridView>, GetPropType<TypeTag, Properties::GridView>,
GetPropType<TypeTag, Properties::DofMapper>, GetPropType<TypeTag, Properties::DofMapper>,
GetPropType<TypeTag, Properties::Stencil>, GetPropType<TypeTag, Properties::Stencil>,
GetPropType<TypeTag, Properties::FluidSystem>,
GetPropType<TypeTag, Properties::Scalar>> GetPropType<TypeTag, Properties::Scalar>>
{ {
using BaseType = GenericTracerModel<GetPropType<TypeTag, Properties::Grid>, using BaseType = GenericTracerModel<GetPropType<TypeTag, Properties::Grid>,
GetPropType<TypeTag, Properties::GridView>, GetPropType<TypeTag, Properties::GridView>,
GetPropType<TypeTag, Properties::DofMapper>, GetPropType<TypeTag, Properties::DofMapper>,
GetPropType<TypeTag, Properties::Stencil>, GetPropType<TypeTag, Properties::Stencil>,
GetPropType<TypeTag, Properties::FluidSystem>,
GetPropType<TypeTag, Properties::Scalar>>; GetPropType<TypeTag, Properties::Scalar>>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>; using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using GridView = GetPropType<TypeTag, Properties::GridView>; using GridView = GetPropType<TypeTag, Properties::GridView>;
@@ -113,7 +115,7 @@ public:
some phase index stuff. If this is a normal run the initial tracer some phase index stuff. If this is a normal run the initial tracer
concentrations will be assigned from the TBLK or TVDPF keywords. concentrations will be assigned from the TBLK or TVDPF keywords.
2. [Restart only:] The tracer concenntration are read from the restart 2. [Restart only:] The tracer concentration are read from the restart
file and the concentrations are applied with repeated calls to the file and the concentrations are applied with repeated calls to the
setTracerConcentration() method. This is currently done in the setTracerConcentration() method. This is currently done in the
eclwriter::beginRestart() method. eclwriter::beginRestart() method.
@@ -144,9 +146,6 @@ public:
if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){ if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
throw std::runtime_error("Oil tracer specified for non-oil fluid system:" + this->name(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->name(tracerIdx));
}
oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]); oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
} }
@@ -154,12 +153,15 @@ public:
if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){ if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
throw std::runtime_error("Gas tracer specified for non-gas fluid system:" + this->name(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->name(tracerIdx));
}
gas_.addTracer(tracerIdx, this->tracerConcentration_[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_ // will be valid after we move out of tracerMatrix_
@@ -220,39 +222,63 @@ public:
protected: protected:
// evaluate water storage volume(s) in a single cell // compute volume associated with free concentration
template <class LhsEval> Scalar computeFreeVolume_(const int tracerPhaseIdx,
void computeVolume_(LhsEval& freeVolume, unsigned globalDofIdx,
const int tracerPhaseIdx, unsigned timeIdx)
const ElementContext& elemCtx,
unsigned scvIdx,
unsigned timeIdx)
{ {
const auto& intQuants = elemCtx.intensiveQuantities(scvIdx, timeIdx); const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx);
const auto& fs = intQuants.fluidState(); const auto& fs = intQuants.fluidState();
Scalar phaseVolume = Scalar phaseVolume =
decay<Scalar>(fs.saturation(tracerPhaseIdx)) decay<Scalar>(fs.saturation(tracerPhaseIdx))
*decay<Scalar>(fs.invB(tracerPhaseIdx)) *decay<Scalar>(fs.invB(tracerPhaseIdx))
*decay<Scalar>(intQuants.porosity()); *decay<Scalar>(intQuants.porosity());
// avoid singular matrix if no water is present. return max(phaseVolume, 1e-10);
phaseVolume = max(phaseVolume, 1e-10); }
if (std::is_same<LhsEval, Scalar>::value) // compute volume associated with solution concentration
freeVolume = phaseVolume; Scalar computeSolutionVolume_(const int tracerPhaseIdx,
else unsigned globalDofIdx,
freeVolume = phaseVolume * variable<LhsEval>(1.0, 0); unsigned timeIdx)
{
const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx);
const auto& fs = intQuants.fluidState();
Scalar phaseVolume;
// vaporized oil
if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
phaseVolume =
decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx))
* decay<Scalar>(fs.invB(FluidSystem::gasPhaseIdx))
* decay<Scalar>(fs.Rv())
* decay<Scalar>(intQuants.porosity());
}
// dissolved gas
else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
phaseVolume =
decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx))
* decay<Scalar>(fs.invB(FluidSystem::oilPhaseIdx))
* decay<Scalar>(fs.Rs())
* decay<Scalar>(intQuants.porosity());
}
else {
phaseVolume = 0.0;
}
return max(phaseVolume, 1e-10);
} }
// evaluate the flux(es) over one face void computeFreeFlux_(TracerEvaluation & freeFlux,
void computeFlux_(TracerEvaluation & freeFlux, bool & isUp,
bool & isUpFree, const int tracerPhaseIdx,
const int tracerPhaseIdx, const ElementContext& elemCtx,
const ElementContext& elemCtx, unsigned scvfIdx,
unsigned scvfIdx, unsigned timeIdx)
unsigned timeIdx)
{ {
const auto& stencil = elemCtx.stencil(timeIdx); const auto& stencil = elemCtx.stencil(timeIdx);
const auto& scvf = stencil.interiorFace(scvfIdx); const auto& scvf = stencil.interiorFace(scvfIdx);
@@ -265,17 +291,72 @@ protected:
const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx); const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
const auto& fs = intQuants.fluidState(); const auto& fs = intQuants.fluidState();
Scalar A = scvf.area(); Scalar v =
Scalar v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx)); decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx))
Scalar b = decay<Scalar>(fs.invB(tracerPhaseIdx)); * decay<Scalar>(fs.invB(tracerPhaseIdx));
Scalar A = scvf.area();
if (inIdx == upIdx) { if (inIdx == upIdx) {
freeFlux = A*v*b*variable<TracerEvaluation>(1.0, 0); freeFlux = A*v*variable<TracerEvaluation>(1.0, 0);
isUpFree = true; isUp = true;
} }
else { else {
freeFlux = A*v*b*1.0; freeFlux = A*v;
isUpFree = false; 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<Scalar>(fs.invB(FluidSystem::gasPhaseIdx))
* decay<Scalar>(extQuants.volumeFlux(FluidSystem::gasPhaseIdx))
* decay<Scalar>(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<Scalar>(fs.invB(FluidSystem::oilPhaseIdx))
* decay<Scalar>(extQuants.volumeFlux(FluidSystem::oilPhaseIdx))
* decay<Scalar>(fs.Rs());
}
else {
upIdx = 0;
v = 0.0;
}
Scalar A = scvf.area();
if (inIdx == upIdx) {
solFlux = A*v*variable<TracerEvaluation>(1.0, 0);
isUp = true;
}
else {
solFlux = A*v;
isUp = false;
} }
} }
@@ -290,29 +371,46 @@ protected:
{ {
if (tr.numTracer() == 0) if (tr.numTracer() == 0)
return; return;
// Storage terms at previous time step (timeIdx = 1)
std::vector<Scalar> storageOfTimeIndex1(tr.numTracer()); std::vector<Scalar> storageOfTimeIndex1(tr.numTracer());
std::vector<Scalar> fStorageOfTimeIndex1(tr.numTracer());
std::vector<Scalar> sStorageOfTimeIndex1(tr.numTracer());
if (elemCtx.enableStorageCache()) { if (elemCtx.enableStorageCache()) {
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { 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 { else {
Scalar fVolume1; Scalar fVolume1 = computeFreeVolume_(tr.phaseIdx_, I1, 1);
computeVolume_(fVolume1, tr.phaseIdx_, elemCtx, 0, /*timeIdx=*/1); Scalar sVolume1 = computeSolutionVolume_(tr.phaseIdx_, I1, 1);
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { 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; TracerEvaluation fVol = computeFreeVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, /*timeIdx=*/0); TracerEvaluation sVol = computeSolutionVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(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) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
Scalar storageOfTimeIndex0 = fVolume.value()*tr.concentration_[tIdx][I]; // Free part
Scalar localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1[tIdx]) * scvVolume/dt; Scalar fStorageOfTimeIndex0 = fVol.value() * tr.concentration_[tIdx][I][0];
tr.residual_[tIdx][I][0] += localStorage; //residual + flux 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<class TrRe> template<class TrRe>
@@ -320,21 +418,36 @@ protected:
const ElementContext& elemCtx, const ElementContext& elemCtx,
unsigned scvfIdx, unsigned scvfIdx,
unsigned I, unsigned I,
unsigned J) unsigned J,
const Scalar dt)
{ {
if (tr.numTracer() == 0) if (tr.numTracer() == 0)
return; return;
TracerEvaluation flux; TracerEvaluation fFlux;
TracerEvaluation sFlux;
bool isUpF; bool isUpF;
computeFlux_(flux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0); bool isUpS;
int globalUpIdx = isUpF ? I : J; computeFreeFlux_(fFlux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0);
computeSolFlux_(sFlux, isUpS, tr.phaseIdx_, elemCtx, scvfIdx, 0);
dsVol_[tr.phaseIdx_][I] += sFlux.value() * dt;
dfVol_[tr.phaseIdx_][I] += fFlux.value() * dt;
int fGlobalUpIdx = isUpF ? I : J;
int sGlobalUpIdx = isUpS ? I : J;
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) { for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
tr.residual_[tIdx][I][0] += flux.value()*tr.concentration_[tIdx][globalUpIdx]; //residual + flux // Free and solution fluxes
tr.residual_[tIdx][I][0] += fFlux.value()*tr.concentration_[tIdx][fGlobalUpIdx][0]; //residual + flux
tr.residual_[tIdx][I][1] += sFlux.value()*tr.concentration_[tIdx][sGlobalUpIdx][1]; //residual + flux
} }
if (isUpF) {
(*tr.mat)[J][I][0][0] = -flux.derivative(0); // Derivative matrix
(*tr.mat)[I][I][0][0] += flux.derivative(0); if (isUpF){
(*tr.mat)[J][I][0][0] = -fFlux.derivative(0);
(*tr.mat)[I][I][0][0] += fFlux.derivative(0);
}
if (isUpS) {
(*tr.mat)[J][I][1][1] = -sFlux.derivative(0);
(*tr.mat)[I][I][1][1] += sFlux.derivative(0);
} }
} }
@@ -346,33 +459,124 @@ protected:
return; return;
const auto& eclWell = well.wellEcl(); const auto& eclWell = well.wellEcl();
// Init. well output to zero
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
this->wellTracerRate_[std::make_pair(eclWell.name(), this->name(tr.idx_[tIdx]))] = 0.0; this->wellTracerRate_[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<double> wtracer(tr.numTracer()); std::vector<Scalar> wtracer(tr.numTracer());
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
wtracer[tIdx] = this->currentConcentration_(eclWell, this->name(tr.idx_[tIdx])); wtracer[tIdx] = this->currentConcentration_(eclWell, this->name(tr.idx_[tIdx]));
} }
for (auto& perfData : well.perforationData()) { Scalar dt = simulator_.timeStepSize();
const auto I = perfData.cell_index; 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_); Scalar rate = well.volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
if (rate > 0) { Scalar rate_s;
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
tr.residual_[tIdx][I][0] -= rate*wtracer[tIdx]; rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
// Store _injector_ tracer rate for reporting
this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) += rate*wtracer[tIdx];
}
} }
else if (rate < 0) { 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) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
tr.residual_[tIdx][I][0] -= rate*tr.concentration_[tIdx][I]; // 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];
}
} }
(*tr.mat)[I][I][0][0] -= rate*variable<TracerEvaluation>(1.0, 0).derivative(0); 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<TracerEvaluation>(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];
}
dsVol_[tr.phaseIdx_][I] -= rate_s * dt;
// Derivative matrix for solution tracer producer
(*tr.mat)[I][I][1][1] -= rate_s * variable<TracerEvaluation>(1.0, 0).derivative(0);
} }
} }
} }
template<class TrRe>
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<TracerEvaluation>(1.0, 0).derivative(0);
(*tr.mat)[I][I][1][0] += (dfVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
}
else {
(*tr.mat)[I][I][0][1] += (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
(*tr.mat)[I][I][1][1] -= (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
}
}
void assembleTracerEquations_() void assembleTracerEquations_()
{ {
@@ -385,8 +589,17 @@ protected:
for (auto& tr : tbatch) { for (auto& tr : tbatch) {
if (tr.numTracer() != 0) { if (tr.numTracer() != 0) {
(*tr.mat) = 0.0; (*tr.mat) = 0.0;
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
tr.residual_[tIdx] = 0.0; tr.residual_[tIdx] = 0.0;
}
}
}
// Well terms
const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
for (const auto& wellPtr : wellPtrs) {
for (auto& tr : tbatch) {
this->assembleTracerEquationWell(tr, *wellPtr);
} }
} }
@@ -400,8 +613,10 @@ protected:
{ {
// Dirichlet boundary conditions needed for the parallel matrix // Dirichlet boundary conditions needed for the parallel matrix
for (auto& tr : tbatch) { for (auto& tr : tbatch) {
if (tr.numTracer() != 0) if (tr.numTracer() != 0) {
(*tr.mat)[I][I][0][0] = 1.; (*tr.mat)[I][I][0][0] = 1.;
(*tr.mat)[I][I][1][1] = 1.;
}
} }
continue; continue;
} }
@@ -430,17 +645,15 @@ protected:
unsigned j = face.exteriorIndex(); unsigned j = face.exteriorIndex();
unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0); unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
for (auto& tr : tbatch) { for (auto& tr : tbatch) {
this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J); this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J, dt);
} }
} }
}
// Source terms (mass transfer between free and solution tracer)
// Well terms
const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
for (const auto& wellPtr : wellPtrs) {
for (auto& tr : tbatch) { for (auto& tr : tbatch) {
this->assembleTracerEquationWell(tr, *wellPtr); this->assembleTracerEquationSource(tr, dt, I);
} }
} }
// Communicate overlap using grid Communication // Communicate overlap using grid Communication
@@ -457,22 +670,31 @@ protected:
void updateStorageCache() void updateStorageCache()
{ {
for (auto& tr : tbatch) { for (auto& tr : tbatch) {
if (tr.numTracer() != 0) if (tr.numTracer() != 0) {
tr.concentrationInitial_ = tr.concentration_; tr.concentrationInitial_ = tr.concentration_;
}
} }
ElementContext elemCtx(simulator_); ElementContext elemCtx(simulator_);
for (const auto& elem : elements(simulator_.gridView())) { for (const auto& elem : elements(simulator_.gridView())) {
elemCtx.updatePrimaryStencil(elem); elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
Scalar extrusionFactor = elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume() * extrusionFactor;
int globalDofIdx = elemCtx.globalSpaceIndex(0, /*timeIdx=*/0); int globalDofIdx = elemCtx.globalSpaceIndex(0, /*timeIdx=*/0);
for (auto& tr : tbatch) { for (auto& tr : tbatch) {
if (tr.numTracer() == 0) if (tr.numTracer() == 0)
continue; continue;
Scalar fVolume;
computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, /*timeIdx=*/0); Scalar fVol1 = computeFreeVolume_(tr.phaseIdx_, globalDofIdx, 0);
Scalar sVol1 = computeSolutionVolume_(tr.phaseIdx_, globalDofIdx, 0);
fVol1_[tr.phaseIdx_][globalDofIdx] = fVol1 * scvVolume;
sVol1_[tr.phaseIdx_][globalDofIdx] = sVol1 * scvVolume;
dsVol_[tr.phaseIdx_][globalDofIdx] = 0.0;
dfVol_[tr.phaseIdx_][globalDofIdx] = 0.0;
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
tr.storageOfTimeIndex1_[tIdx][globalDofIdx] = fVolume*tr.concentrationInitial_[tIdx][globalDofIdx]; tr.storageOfTimeIndex1_[tIdx][globalDofIdx][0] = fVol1 * tr.concentrationInitial_[tIdx][globalDofIdx][0];
tr.storageOfTimeIndex1_[tIdx][globalDofIdx][1] = sVol1 * tr.concentrationInitial_[tIdx][globalDofIdx][1];
} }
} }
} }
@@ -489,8 +711,9 @@ protected:
// Note that we solve for a concentration update (compared to previous time step) // Note that we solve for a concentration update (compared to previous time step)
// Confer also assembleTracerEquations_(...) above. // Confer also assembleTracerEquations_(...) above.
std::vector<TracerVector> dx(tr.concentration_); std::vector<TracerVector> dx(tr.concentration_);
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
dx[tIdx] = 0.0; dx[tIdx] = 0.0;
}
bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_); bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_);
if (!converged) { if (!converged) {
@@ -498,29 +721,98 @@ protected:
} }
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
tr.concentration_[tIdx] -= dx[tIdx]; for (std::size_t globalDofIdx = 0; globalDofIdx < tr.concentration_[tIdx].size(); ++globalDofIdx) {
// Tracer concentrations for restart report // New concetration. Concentrations that are negative or where free/solution phase is not
this->tracerConcentration_[tr.idx_[tIdx]] = tr.concentration_[tIdx]; // present are set to zero
const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, 0);
const auto& fs = intQuants.fluidState();
Scalar Sf = decay<Scalar>(fs.saturation(tr.phaseIdx_));
Scalar Ss = 0.0;
if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas())
Ss = decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx));
else if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil())
Ss = decay<Scalar>(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 // Store _producer_ tracer rate for reporting
const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells(); const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
for (const auto& wellPtr : wellPtrs) { 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; continue;
Scalar rateWellPos = 0.0; Scalar rateWellPos = 0.0;
Scalar rateWellNeg = 0.0; Scalar rateWellNeg = 0.0;
for (auto& perfData : wellPtr->perforationData()) { std::size_t well_index = simulator_.problem().wellModel().wellState().index(eclWell.name()).value();
const int I = perfData.cell_index; 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 = 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) { if (rate < 0) {
rateWellNeg += rate; rateWellNeg += rate;
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];
}
} }
else { else {
rateWellPos += rate; rateWellPos += rate;
@@ -532,7 +824,6 @@ protected:
// TODO: Some inconsistencies here that perhaps should be clarified. The "offical" rate as reported below is // TODO: Some inconsistencies here that perhaps should be clarified. The "offical" rate as reported below is
// occasionally significant different from the sum over connections (as calculated above). Only observed // occasionally significant different from the sum over connections (as calculated above). Only observed
// for small values, neglible for the rate itself, but matters when used to calculate tracer concentrations. // for small values, neglible for the rate itself, but matters when used to calculate tracer concentrations.
std::size_t well_index = simulator_.problem().wellModel().wellState().index(well.name()).value();
Scalar official_well_rate_total = simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_]; Scalar official_well_rate_total = simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_];
rateWellTotal = official_well_rate_total; rateWellTotal = official_well_rate_total;
@@ -541,13 +832,14 @@ protected:
const Scalar bucketPrDay = 10.0/(1000.*3600.*24.); // ... keeps (some) trouble away const Scalar bucketPrDay = 10.0/(1000.*3600.*24.); // ... keeps (some) trouble away
const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal/rateWellNeg : 0.0; const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal/rateWellNeg : 0.0;
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
this->wellTracerRate_.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_; Simulator& simulator_;
@@ -561,13 +853,13 @@ protected:
template <typename TV> template <typename TV>
struct TracerBatch { struct TracerBatch {
std::vector<int> idx_; std::vector<int> idx_;
const int phaseIdx_; const int phaseIdx_;
std::vector<TV> concentrationInitial_; std::vector<TV> concentrationInitial_;
std::vector<TV> concentration_; std::vector<TV> concentration_;
std::vector<TV> storageOfTimeIndex1_; std::vector<TV> storageOfTimeIndex1_;
std::vector<TV> residual_; std::vector<TV> residual_;
std::unique_ptr<TracerMatrix> mat; std::unique_ptr<TracerMatrix> mat;
bool operator==(const TracerBatch& rhs) const bool operator==(const TracerBatch& rhs) const
{ {
@@ -600,12 +892,12 @@ protected:
void addTracer(const int idx, const TV & concentration) void addTracer(const int idx, const TV & concentration)
{ {
int numGridDof = concentration.size(); int numGridDof = concentration.size();
idx_.emplace_back(idx); idx_.emplace_back(idx);
concentrationInitial_.emplace_back(concentration); concentrationInitial_.emplace_back(concentration);
concentration_.emplace_back(concentration); concentration_.emplace_back(concentration);
storageOfTimeIndex1_.emplace_back(numGridDof); residual_.emplace_back(numGridDof);
residual_.emplace_back(numGridDof); storageOfTimeIndex1_.emplace_back(numGridDof);
} }
}; };
@@ -613,6 +905,10 @@ protected:
TracerBatch<TracerVector>& wat_; TracerBatch<TracerVector>& wat_;
TracerBatch<TracerVector>& oil_; TracerBatch<TracerVector>& oil_;
TracerBatch<TracerVector>& gas_; TracerBatch<TracerVector>& gas_;
std::array<std::vector<Scalar>, 3> fVol1_;
std::array<std::vector<Scalar>, 3> sVol1_;
std::array<std::vector<Scalar>, 3> dsVol_;
std::array<std::vector<Scalar>, 3> dfVol_;
}; };
} // namespace Opm } // namespace Opm

47
opm/simulators/flow/VtkTracerModule.hpp
View File

@@ -107,10 +107,14 @@ namespace Opm {
{ {
if (eclTracerConcentrationOutput_()){ if (eclTracerConcentrationOutput_()){
const auto& tracerModel = this->simulator_.problem().tracerModel(); const auto& tracerModel = this->simulator_.problem().tracerModel();
eclTracerConcentration_.resize(tracerModel.numTracers()); eclFreeTracerConcentration_.resize(tracerModel.numTracers());
for (std::size_t tracerIdx = 0; tracerIdx < eclTracerConcentration_.size(); ++tracerIdx) { eclSolTracerConcentration_.resize(tracerModel.numTracers());
const auto& enableSolTracers = tracerModel.enableSolTracers();
this->resizeScalarBuffer_(eclTracerConcentration_[tracerIdx]); 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<TypeTag, Properties::EnableVtkOutput>()) if (!Parameters::get<TypeTag, Properties::EnableVtkOutput>())
return; return;
const auto& tracerModel = elemCtx.problem().tracerModel(); if (eclTracerConcentrationOutput_()) {
const auto& tracerModel = elemCtx.problem().tracerModel();
const auto& enableSolTracers = tracerModel.enableSolTracers();
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { for (std::size_t tracerIdx = 0; tracerIdx < eclFreeTracerConcentration_.size(); ++tracerIdx) {
unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); // free tracer
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
if (eclTracerConcentrationOutput_()){ unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
for (std::size_t tracerIdx = 0; tracerIdx < eclTracerConcentration_.size(); ++tracerIdx) { eclFreeTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx);
eclTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.tracerConcentration(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_()){ if (eclTracerConcentrationOutput_()){
const auto& tracerModel = this->simulator_.problem().tracerModel(); const auto& tracerModel = this->simulator_.problem().tracerModel();
for (std::size_t tracerIdx = 0; tracerIdx < eclTracerConcentration_.size(); ++tracerIdx) { const auto& enableSolTracers = tracerModel.enableSolTracers();
const std::string tmp = "tracerConcentration_" + tracerModel.name(tracerIdx);
this->commitScalarBuffer_(baseWriter,tmp.c_str(), eclTracerConcentration_[tracerIdx]); 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<ScalarBuffer> eclTracerConcentration_; std::vector<ScalarBuffer> eclFreeTracerConcentration_;
std::vector<ScalarBuffer> eclSolTracerConcentration_;
}; };
} // namespace Opm } // namespace Opm

9
opm/simulators/wells/BlackoilWellModel.hpp
View File

@@ -267,8 +267,17 @@ template<class Scalar> class WellContributions;
{ {
const auto& tracerRates = this->simulator_.problem() const auto& tracerRates = this->simulator_.problem()
.tracerModel().getWellTracerRates(); .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, tracerRates);
this->assignWellTracerRates(wsrpt, freeTracerRates);
this->assignWellTracerRates(wsrpt, solTracerRates);
this->assignMswTracerRates(wsrpt, mswTracerRates);
} }
BlackoilWellModelGuideRates(*this) BlackoilWellModelGuideRates(*this)

27
opm/simulators/wells/BlackoilWellModelGeneric.cpp
View File

@@ -1698,6 +1698,33 @@ assignWellTracerRates(data::Wells& wsrpt,
} }
} }
template<class Scalar>
void BlackoilWellModelGeneric<Scalar>::
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<class Scalar> template<class Scalar>
std::vector<std::vector<int>> BlackoilWellModelGeneric<Scalar>:: std::vector<std::vector<int>> BlackoilWellModelGeneric<Scalar>::
getMaxWellConnections() const getMaxWellConnections() const

3
opm/simulators/wells/BlackoilWellModelGeneric.hpp
View File

@@ -437,6 +437,9 @@ protected:
using WellTracerRates = std::map<std::pair<std::string, std::string>, Scalar>; using WellTracerRates = std::map<std::pair<std::string, std::string>, Scalar>;
void assignWellTracerRates(data::Wells& wsrpt, void assignWellTracerRates(data::Wells& wsrpt,
const WellTracerRates& wellTracerRates) const; const WellTracerRates& wellTracerRates) const;
using MswTracerRates = std::map<std::tuple<std::string, std::string, std::size_t>, Scalar>;
void assignMswTracerRates(data::Wells& wsrpt,
const MswTracerRates& mswTracerRates) const;
Schedule& schedule_; Schedule& schedule_;
const SummaryState& summaryState_; const SummaryState& summaryState_;

10
tests/test_RestartSerialization.cpp
View File

@@ -381,10 +381,10 @@ BOOST_AUTO_TEST_CASE(BlackoilWellModelGeneric)
BOOST_CHECK_MESSAGE(data_out == data_in, "Deserialized BlackoilWellModelGeneric differ"); BOOST_CHECK_MESSAGE(data_out == data_in, "Deserialized BlackoilWellModelGeneric differ");
} }
template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar> template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
class GenericTracerModelTest : public Opm::GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar> class GenericTracerModelTest : public Opm::GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>
{ {
using Base = Opm::GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>; using Base = Opm::GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>;
public: public:
GenericTracerModelTest(const GridView& gridView, GenericTracerModelTest(const GridView& gridView,
const Opm::EclipseState& eclState, const Opm::EclipseState& eclState,
@@ -402,7 +402,7 @@ public:
const std::function<std::array<double,Grid::dimensionworld>(int)> centroids) const std::function<std::array<double,Grid::dimensionworld>(int)> centroids)
{ {
GenericTracerModelTest result(gridView, eclState, cartMapper, dofMapper, 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}); result.wellTracerRate_.insert({{"foo", "bar"}, 4.0});
return result; return result;
@@ -438,6 +438,7 @@ BOOST_AUTO_TEST_CASE(FlowGenericTracerModel)
GridView, GridView,
Dune::MultipleCodimMultipleGeomTypeMapper<GridView>, Dune::MultipleCodimMultipleGeomTypeMapper<GridView>,
Opm::EcfvStencil<double, GridView, false, false>, Opm::EcfvStencil<double, GridView, false, false>,
Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits>,
double> double>
::serializationTestObject(gridView, eclState, mapper, dofMapper, centroids); ::serializationTestObject(gridView, eclState, mapper, dofMapper, centroids);
Opm::Serialization::MemPacker packer; Opm::Serialization::MemPacker packer;
@@ -474,6 +475,7 @@ BOOST_AUTO_TEST_CASE(FlowGenericTracerModelFem)
GridView, GridView,
Dune::MultipleCodimMultipleGeomTypeMapper<GridView>, Dune::MultipleCodimMultipleGeomTypeMapper<GridView>,
Opm::EcfvStencil<double, GridView, false, false>, Opm::EcfvStencil<double, GridView, false, false>,
Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits>,
double> double>
::serializationTestObject(gridView, eclState, mapper, dofMapper, centroids); ::serializationTestObject(gridView, eclState, mapper, dofMapper, centroids);
Opm::Serialization::MemPacker packer; Opm::Serialization::MemPacker packer;