OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-07-04 11:33:06 -05:00
Code Issues Packages Projects Releases Wiki Activity

Fix solution tracers and well output.

Browse Source 
-Only output or restart solution tracers for gas/oil tracers with DISGAS/VAPOIL enabled (no solution tracers in water phase!).
-Initial tracers (free/solution) will be set to zero initially if TBLK/TVDP is not given.
- Do not calculate mass transfer between free and solution tracers if it is not necessary.
-Calculate well rates using updated tracer concentrations
This commit is contained in:
Svenn Tveit 2024-06-10 14:38:06 +02:00
parent 61cfcfa523
commit 84cdef1135
10 changed files with 244 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -501,8 +501,10 @@ public:
{
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 @@ public:
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<double>(free_tracer_name);
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.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<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);
}
}
}

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

@ -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<bool>& 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);
}
}

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

@ -340,6 +340,7 @@ protected:
const bool vapparsActive,
const bool enableHysteresis,
unsigned numTracers,
const std::vector<bool>& enableSolTracers,
unsigned numOutputNnc);
void makeRegionSum(Inplace& inplace,

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

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

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

@ -36,6 +36,8 @@
#include <opm/simulators/linalg/matrixblock.hh>
#include <opm/input/eclipse/EclipseState/Phase.hpp>
#include <array>
#include <cstddef>
#include <functional>
@ -49,7 +51,7 @@ namespace Opm {
class EclipseState;
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 {
public:
using TracerVectorSingle = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
@ -71,6 +73,8 @@ public:
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
@ -79,6 +83,7 @@ public:
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 @@ public:
void serializeOp(Serializer& serializer)
{
serializer(tracerConcentration_);
serializer(freeTracerConcentration_);
serializer(solTracerConcentration_);
serializer(wellTracerRate_);
serializer(wellFreeTracerRate_);
serializer(wellSolTracerRate_);
serializer(mSwTracerRate_);
}
protected:
@ -131,6 +139,7 @@ protected:
const DofMapper& dofMapper_;
std::vector<int> tracerPhaseIdx_;
std::vector<bool> enableSolTracers_;
std::vector<TracerVector> tracerConcentration_;
std::unique_ptr<TracerMatrix> tracerMatrix_;
std::vector<TracerVectorSingle> freeTracerConcentration_;

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

@ -39,7 +39,6 @@
#include <opm/grid/CpGrid.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/Schedule/Well/Well.hpp>
#include <opm/input/eclipse/Schedule/Well/WellTracerProperties.hpp>
@ -97,8 +96,8 @@ createParallelFlexibleSolver(const Dune::CpGrid& grid, const Matrix& M, const Pr
}
#endif
template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar>
GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
GenericTracerModel(const GridView& gridView,
const EclipseState& eclState,
const CartesianIndexMapper& cartMapper,
@ -112,8 +111,8 @@ GenericTracerModel(const GridView& gridView,
{
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
freeTracerConcentration(int tracerIdx, int globalDofIdx) const
{
if (freeTracerConcentration_.empty())
@ -122,8 +121,8 @@ freeTracerConcentration(int tracerIdx, int globalDofIdx) const
return freeTracerConcentration_[tracerIdx][globalDofIdx];
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
solTracerConcentration(int tracerIdx, int globalDofIdx) const
{
if (solTracerConcentration_.empty())
@ -132,71 +131,94 @@ solTracerConcentration(int tracerIdx, int globalDofIdx) const
return solTracerConcentration_[tracerIdx][globalDofIdx];
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
void GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class 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 Scalar>
void GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
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 Scalar>
int GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
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
{
return this->eclState_.tracer().size();
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
fname(int tracerIdx) const
{
return this->eclState_.tracer()[tracerIdx].fname();
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class 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 Scalar>
std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
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 Scalar>
std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
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 Scalar>
Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
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
{
return eclWell.getTracerProperties().getConcentration(name);
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
const std::string& GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
const std::string& GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
name(int tracerIdx) const
{
return this->eclState_.tracer()[tracerIdx].name;
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
void GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
void GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
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();
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 {
throw std::logic_error(fmt::format("Can not initialize solution tracer concentration: {}", tracer.name));
// 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;
}
}
}
@ -331,8 +380,8 @@ doInit(bool rst, std::size_t numGridDof,
tracerMatrix_->endindices();
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
bool GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
bool GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
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<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
bool GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
bool GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b)
{
Scalar tolerance = 1e-6;
Scalar tolerance = 1e-2;
int maxIter = 100;
int verbosity = 0;

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

@ -240,6 +240,7 @@ public:
problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)),
problem.materialLawManager()->enableHysteresis(),
problem.tracerModel().numTracers(),
problem.tracerModel().enableSolTracers(),
problem.eclWriter()->getOutputNnc().size());
}

89
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::DofMapper>,
GetPropType<TypeTag, Properties::Stencil>,
GetPropType<TypeTag, Properties::FluidSystem>,
GetPropType<TypeTag, Properties::Scalar>>
{
using BaseType = GenericTracerModel<GetPropType<TypeTag, Properties::Grid>,
GetPropType<TypeTag, Properties::GridView>,
GetPropType<TypeTag, Properties::DofMapper>,
GetPropType<TypeTag, Properties::Stencil>,
GetPropType<TypeTag, Properties::FluidSystem>,
GetPropType<TypeTag, Properties::Scalar>>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
@ -494,7 +496,7 @@ protected:
// Injection of free tracer only
tr.residual_[tIdx][I][0] -= rate_f*wtracer[tIdx];
// Store _injector_ tracer rate for reporting
// 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()) {
@ -506,39 +508,23 @@ protected:
}
else if (rate_f < 0) {
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
// Free and solution tracer production
// Production of free tracer
tr.residual_[tIdx][I][0] -= rate_f * tr.concentration_[tIdx][I][0];
// Store _producer_ 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];
}
}
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<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];
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<TracerEvaluation>(1.0, 0).derivative(0);
}
}
@ -552,6 +538,13 @@ protected:
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];
@ -753,16 +746,59 @@ protected:
// Store _producer_ tracer rate for reporting
const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
for (const auto& wellPtr : wellPtrs) {
const auto& well = wellPtr->wellEcl();
const auto& eclWell = wellPtr->wellEcl();
if (!well.isProducer()) //Injection rates already reported during assembly
// Injection rates already reported during assembly
if (!eclWell.isProducer())
continue;
Scalar rateWellPos = 0.0;
Scalar rateWellNeg = 0.0;
for (auto& perfData : wellPtr->perforationData()) {
const int I = perfData.cell_index;
std::size_t well_index = simulator_.problem().wellModel().wellState().index(eclWell.name()).value();
const auto& ws = simulator_.problem().wellModel().wellState().well(well_index);
for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
const auto I = ws.perf_data.cell_index[i];
Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
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 @@ protected:
// 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
// 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_];
rateWellTotal = official_well_rate_total;
@ -785,7 +820,7 @@ protected:
const Scalar bucketPrDay = 10.0/(1000.*3600.*24.); // ... keeps (some) trouble away
const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal/rateWellNeg : 0.0;
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
this->wellTracerRate_.at(std::make_pair(well.name(),this->name(tr.idx_[tIdx]))) *= factor;
this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) *= factor;
}
}
}

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

@ -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<TypeTag, Properties::EnableVtkOutput>())
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) {
unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
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]);
}
}
}

8
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");
}
template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar>
class GenericTracerModelTest : public Opm::GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>
template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class 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:
GenericTracerModelTest(const GridView& gridView,
const Opm::EclipseState& eclState,
@ -438,6 +438,7 @@ BOOST_AUTO_TEST_CASE(FlowGenericTracerModel)
GridView,
Dune::MultipleCodimMultipleGeomTypeMapper<GridView>,
Opm::EcfvStencil<double, GridView, false, false>,
Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits>,
double>
::serializationTestObject(gridView, eclState, mapper, dofMapper, centroids);
Opm::Serialization::MemPacker packer;
@ -474,6 +475,7 @@ BOOST_AUTO_TEST_CASE(FlowGenericTracerModelFem)
GridView,
Dune::MultipleCodimMultipleGeomTypeMapper<GridView>,
Opm::EcfvStencil<double, GridView, false, false>,
Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits>,
double>
::serializationTestObject(gridView, eclState, mapper, dofMapper, centroids);
Opm::Serialization::MemPacker packer;