diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 0308917ea..0289aebd1 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -26,6 +26,7 @@ list (APPEND MAIN_SOURCE_FILES
ebos/collecttoiorank.cc
ebos/eclgenericcpgridvanguard.cc
ebos/eclgenericthresholdpressure.cc
+ ebos/eclgenerictracermodel.cc
ebos/eclgenericvanguard.cc
ebos/ecltransmissibility.cc
opm/core/props/phaseUsageFromDeck.cpp
diff --git a/ebos/eclgenerictracermodel.cc b/ebos/eclgenerictracermodel.cc
new file mode 100644
index 000000000..63f269118
--- /dev/null
+++ b/ebos/eclgenerictracermodel.cc
@@ -0,0 +1,265 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 2 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+
+ Consult the COPYING file in the top-level source directory of this
+ module for the precise wording of the license and the list of
+ copyright holders.
+*/
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#if HAVE_DUNE_FEM
+#include
+#include
+#include
+#endif
+
+#include
+#include
+#include
+
+namespace Opm {
+
+template
+EclGenericTracerModel::
+EclGenericTracerModel(const GridView& gridView,
+ const EclipseState& eclState,
+ const CartesianIndexMapper& cartMapper,
+ const DofMapper& dofMapper)
+ : gridView_(gridView)
+ , eclState_(eclState)
+ , cartMapper_(cartMapper)
+ , dofMapper_(dofMapper)
+{
+}
+
+template
+const std::string& EclGenericTracerModel::
+tracerName(int tracerIdx) const
+{
+ if (tracerNames_.empty())
+ throw std::logic_error("This method should never be called when there are no tracers in the model");
+
+ return tracerNames_[tracerIdx];
+}
+
+template
+Scalar EclGenericTracerModel::
+tracerConcentration(int tracerIdx, int globalDofIdx) const
+{
+ if (tracerConcentration_.empty())
+ return 0.0;
+
+ return tracerConcentration_[tracerIdx][globalDofIdx];
+}
+
+template
+void EclGenericTracerModel::
+doInit(bool enabled, size_t numGridDof,
+ size_t gasPhaseIdx, size_t oilPhaseIdx, size_t waterPhaseIdx)
+{
+ const auto& tracers = eclState_.tracer();
+ const auto& comm = gridView_.comm();
+
+ if (tracers.size() == 0)
+ return; // tracer treatment is supposed to be disabled
+
+ if (!enabled) {
+ if (gridView_.comm().rank() == 0) {
+ OpmLog::warning("Keyword TRACERS has only experimental support, and is hence ignored.\n"
+ "The experimental tracer model can still be used, but must be set explicitely.\n"
+ "To use tracers, set the command line option: --enable-tracer-model=true"
+ "\n");
+ }
+ return; // Tracer transport must be enabled by the user
+ }
+
+ if (comm.size() > 1) {
+ tracerNames_.resize(0);
+ if (comm.rank() == 0)
+ std::cout << "Warning: The tracer model currently does not work for parallel runs\n"
+ << std::flush;
+ return;
+ }
+
+ // retrieve the number of tracers from the deck
+ const size_t numTracers = tracers.size();
+ tracerNames_.resize(numTracers);
+ tracerConcentration_.resize(numTracers);
+ storageOfTimeIndex1_.resize(numTracers);
+
+ // the phase where the tracer is
+ tracerPhaseIdx_.resize(numTracers);
+ size_t tracerIdx = 0;
+ for (const auto& tracer : tracers) {
+ tracerNames_[tracerIdx] = tracer.name;
+ if (tracer.phase == Phase::WATER)
+ tracerPhaseIdx_[tracerIdx] = waterPhaseIdx;
+ else if (tracer.phase == Phase::OIL)
+ tracerPhaseIdx_[tracerIdx] = oilPhaseIdx;
+ else if (tracer.phase == Phase::GAS)
+ tracerPhaseIdx_[tracerIdx] = gasPhaseIdx;
+
+ tracerConcentration_[tracerIdx].resize(numGridDof);
+ storageOfTimeIndex1_[tracerIdx].resize(numGridDof);
+
+
+ //TBLK keyword
+ if (!tracer.concentration.empty()){
+ int tblkDatasize = tracer.concentration.size();
+ if (tblkDatasize < cartMapper_.cartesianSize()){
+ throw std::runtime_error("Wrong size of TBLK for" + tracer.name);
+ }
+ for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
+ int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
+ tracerConcentration_[tracerIdx][globalDofIdx] = tracer.concentration[cartDofIdx];
+ }
+ }
+ //TVDPF keyword
+ else {
+ const auto& eclGrid = eclState_.getInputGrid();
+
+ for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
+ int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
+ const auto& center = eclGrid.getCellCenter(cartDofIdx);
+ tracerConcentration_[tracerIdx][globalDofIdx] = tracer.tvdpf.evaluate("TRACER_CONCENTRATION", center[2]);
+ }
+ }
+ ++tracerIdx;
+ }
+
+ // initial tracer concentration
+ tracerConcentrationInitial_ = tracerConcentration_;
+
+ // residual of tracers
+ tracerResidual_.resize(numGridDof);
+
+ // allocate matrix for storing the Jacobian of the tracer residual
+ tracerMatrix_ = new TracerMatrix(numGridDof, numGridDof, TracerMatrix::random);
+
+ // find the sparsity pattern of the tracer matrix
+ using NeighborSet = std::set;
+ std::vector neighbors(numGridDof);
+
+ Stencil stencil(gridView_, dofMapper_);
+ auto elemIt = gridView_.template begin<0>();
+ const auto elemEndIt = gridView_.template end<0>();
+ for (; elemIt != elemEndIt; ++elemIt) {
+ const auto& elem = *elemIt;
+ stencil.update(elem);
+
+ for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
+ unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
+
+ for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
+ unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
+ neighbors[myIdx].insert(neighborIdx);
+ }
+ }
+ }
+
+ // allocate space for the rows of the matrix
+ for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx)
+ tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
+ tracerMatrix_->endrowsizes();
+
+ // fill the rows with indices. each degree of freedom talks to
+ // all of its neighbors. (it also talks to itself since
+ // degrees of freedom are sometimes quite egocentric.)
+ for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
+ typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
+ typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
+ for (; nIt != nEndIt; ++nIt)
+ tracerMatrix_->addindex(dofIdx, *nIt);
+ }
+ tracerMatrix_->endindices();
+
+ const int sizeCartGrid = cartMapper_.cartesianSize();
+ cartToGlobal_.resize(sizeCartGrid);
+ for (unsigned i = 0; i < numGridDof; ++i) {
+ int cartIdx = cartMapper_.cartesianIndex(i);
+ cartToGlobal_[cartIdx] = i;
+ }
+}
+
+template
+bool EclGenericTracerModel::
+linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
+{
+#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
+ Dune::FMatrixPrecision::set_singular_limit(1.e-30);
+ Dune::FMatrixPrecision::set_absolute_limit(1.e-30);
+#endif
+ x = 0.0;
+ Scalar tolerance = 1e-2;
+ int maxIter = 100;
+
+ int verbosity = 0;
+ using TracerSolver = Dune::BiCGSTABSolver;
+ using TracerOperator = Dune::MatrixAdapter;
+ using TracerScalarProduct = Dune::SeqScalarProduct;
+ using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
+
+ TracerOperator tracerOperator(M);
+ TracerScalarProduct tracerScalarProduct;
+ TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
+
+ TracerSolver solver (tracerOperator, tracerScalarProduct,
+ tracerPreconditioner, tolerance, maxIter,
+ verbosity);
+
+ Dune::InverseOperatorResult result;
+ solver.apply(x, b, result);
+
+ // return the result of the solver
+ return result.converged;
+}
+
+#if HAVE_DUNE_FEM
+template class EclGenericTracerModel>>,
+ Dune::MultipleCodimMultipleGeomTypeMapper>>,
+ Opm::EcfvStencil>>,false,false>,
+ double>;
+#else
+template class EclGenericTracerModel>,
+ Dune::MultipleCodimMultipleGeomTypeMapper>,Dune::Impl::MCMGFailLayout>,
+ Opm::EcfvStencil>,false,false>,
+ double>;
+#endif
+
+template class EclGenericTracerModel,
+ Dune::GridView>,Dune::MultipleCodimMultipleGeomTypeMapper>,Dune::Impl::MCMGFailLayout>,
+ Opm::EcfvStencil>,false,false>,
+ double>;
+
+} // namespace Opm
diff --git a/ebos/eclgenerictracermodel.hh b/ebos/eclgenerictracermodel.hh
new file mode 100644
index 000000000..35446051a
--- /dev/null
+++ b/ebos/eclgenerictracermodel.hh
@@ -0,0 +1,105 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 2 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+
+ Consult the COPYING file in the top-level source directory of this
+ module for the precise wording of the license and the list of
+ copyright holders.
+*/
+/**
+ * \file
+ *
+ * \copydoc Opm::EclTracerModel
+ */
+#ifndef EWOMS_ECL_GENERIC_TRACER_MODEL_HH
+#define EWOMS_ECL_GENERIC_TRACER_MODEL_HH
+
+#include
+
+#include
+#include
+
+#include
+
+#include
+
+#include
+#include
+#include
+
+namespace Opm {
+
+class EclipseState;
+
+template
+class EclGenericTracerModel {
+public:
+ using TracerMatrix = Dune::BCRSMatrix>;
+ using TracerVector = Dune::BlockVector>;
+ using CartesianIndexMapper = Dune::CartesianIndexMapper;
+
+ /*!
+ * \brief Return the number of tracers considered by the tracerModel.
+ */
+ int numTracers() const
+ { return tracerNames_.size(); }
+
+ /*!
+ * \brief Return the tracer name
+ */
+ const std::string& tracerName(int tracerIdx) const;
+
+ /*!
+ * \brief Return the tracer concentration for tracer index and global DofIdx
+ */
+ Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const;
+
+protected:
+ EclGenericTracerModel(const GridView& gridView,
+ const EclipseState& eclState,
+ const CartesianIndexMapper& cartMapper,
+ const DofMapper& dofMapper);
+
+ /*!
+ * \brief Initialize all internal data structures needed by the tracer module
+ */
+ void doInit(bool enabled,
+ size_t numGridDof,
+ size_t gasPhaseIdx,
+ size_t oilPhaseIdx,
+ size_t waterPhaseIdx);
+
+ bool linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b);
+
+ const GridView& gridView_;
+ const EclipseState& eclState_;
+ const CartesianIndexMapper& cartMapper_;
+ const DofMapper& dofMapper_;
+
+ std::vector tracerNames_;
+ std::vector tracerPhaseIdx_;
+ std::vector>> tracerConcentration_;
+ std::vector>> tracerConcentrationInitial_;
+ TracerMatrix *tracerMatrix_;
+ TracerVector tracerResidual_;
+ std::vector cartToGlobal_;
+ std::vector>> storageOfTimeIndex1_;
+};
+
+} // namespace Opm
+
+#endif
diff --git a/ebos/ecltracermodel.hh b/ebos/ecltracermodel.hh
index 0b9cd9a39..c77f280bc 100644
--- a/ebos/ecltracermodel.hh
+++ b/ebos/ecltracermodel.hh
@@ -28,20 +28,12 @@
#ifndef EWOMS_ECL_TRACER_MODEL_HH
#define EWOMS_ECL_TRACER_MODEL_HH
-#include
+#include
-#include
-#include
-
-#include
-#include
-#include
-
-#include
+#include
#include
#include
-#include
namespace Opm::Properties {
@@ -62,8 +54,17 @@ namespace Opm {
* TODO: MPI parallelism.
*/
template
-class EclTracerModel
+class EclTracerModel : public EclGenericTracerModel,
+ GetPropType,
+ GetPropType,
+ GetPropType,
+ GetPropType>
{
+ using BaseType = EclGenericTracerModel,
+ GetPropType,
+ GetPropType,
+ GetPropType,
+ GetPropType>;
using Simulator = GetPropType;
using GridView = GetPropType;
using Grid = GetPropType;
@@ -74,7 +75,7 @@ class EclTracerModel
using RateVector = GetPropType;
using Indices = GetPropType;
- typedef DenseAd::Evaluation TracerEvaluation;
+ using TracerEvaluation = DenseAd::Evaluation;
enum { numEq = getPropValue() };
enum { numPhases = FluidSystem::numPhases };
@@ -82,15 +83,13 @@ class EclTracerModel
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
-
- typedef Dune::BCRSMatrix> TracerMatrix;
- typedef Dune::BlockVector> TracerVector;
-
public:
EclTracerModel(Simulator& simulator)
- : simulator_(simulator)
+ : BaseType(simulator.vanguard().gridView(),
+ simulator.vanguard().eclState(),
+ simulator.vanguard().cartesianIndexMapper(),
+ simulator.model().dofMapper())
+ , simulator_(simulator)
{ }
@@ -99,168 +98,17 @@ public:
*/
void init()
{
- const auto& tracers = simulator_.vanguard().eclState().tracer();
- const auto& comm = simulator_.gridView().comm();
-
- if (tracers.size() == 0)
- return; // tracer treatment is supposed to be disabled
-
- if (!EWOMS_GET_PARAM(TypeTag, bool, EnableTracerModel)) {
- if (simulator_.gridView().comm().rank() == 0) {
- OpmLog::warning("Keyword TRACERS has only experimental support, and is hence ignored.\n"
- "The experimental tracer model can still be used, but must be set explicitely.\n"
- "To use tracers, set the command line option: --enable-tracer-model=true"
- "\n");
- }
- return; // Tracer transport must be enabled by the user
- }
-
- if (comm.size() > 1) {
- tracerNames_.resize(0);
- if (comm.rank() == 0)
- std::cout << "Warning: The tracer model currently does not work for parallel runs\n"
- << std::flush;
- return;
- }
-
- // retrieve the number of tracers from the deck
- const size_t numTracers = tracers.size();
- tracerNames_.resize(numTracers);
- tracerConcentration_.resize(numTracers);
- storageOfTimeIndex1_.resize(numTracers);
-
- // the phase where the tracer is
- tracerPhaseIdx_.resize(numTracers);
- size_t numGridDof = simulator_.model().numGridDof();
- size_t tracerIdx = 0;
- for (const auto& tracer : tracers) {
- tracerNames_[tracerIdx] = tracer.name;
- if (tracer.phase == Phase::WATER)
- tracerPhaseIdx_[tracerIdx] = waterPhaseIdx;
- else if (tracer.phase == Phase::OIL)
- tracerPhaseIdx_[tracerIdx] = oilPhaseIdx;
- else if (tracer.phase == Phase::GAS)
- tracerPhaseIdx_[tracerIdx] = gasPhaseIdx;
-
- tracerConcentration_[tracerIdx].resize(numGridDof);
- storageOfTimeIndex1_[tracerIdx].resize(numGridDof);
-
-
- //TBLK keyword
- if (!tracer.concentration.empty()){
- const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper();
- int tblkDatasize = tracer.concentration.size();
- if (tblkDatasize < simulator_.vanguard().cartesianSize()){
- throw std::runtime_error("Wrong size of TBLK for" + tracer.name);
- }
- for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
- int cartDofIdx = cartMapper.cartesianIndex(globalDofIdx);
- tracerConcentration_[tracerIdx][globalDofIdx] = tracer.concentration[cartDofIdx];
- }
- }
- //TVDPF keyword
- else {
- const auto& eclGrid = simulator_.vanguard().eclState().getInputGrid();
- const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper();
-
- for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
- int cartDofIdx = cartMapper.cartesianIndex(globalDofIdx);
- const auto& center = eclGrid.getCellCenter(cartDofIdx);
- tracerConcentration_[tracerIdx][globalDofIdx] = tracer.tvdpf.evaluate("TRACER_CONCENTRATION", center[2]);
- }
- }
- ++tracerIdx;
- }
-
- // initial tracer concentration
- tracerConcentrationInitial_ = tracerConcentration_;
-
- // residual of tracers
- tracerResidual_.resize(numGridDof);
-
- // allocate matrix for storing the Jacobian of the tracer residual
- tracerMatrix_ = new TracerMatrix(numGridDof, numGridDof, TracerMatrix::random);
-
- // find the sparsity pattern of the tracer matrix
- typedef std::set NeighborSet;
- std::vector neighbors(numGridDof);
-
- Stencil stencil(simulator_.gridView(), simulator_.model().dofMapper() );
- ElementIterator elemIt = simulator_.gridView().template begin<0>();
- const ElementIterator elemEndIt = simulator_.gridView().template end<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- const Element& elem = *elemIt;
- stencil.update(elem);
-
- for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
- unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
-
- for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
- unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
- neighbors[myIdx].insert(neighborIdx);
- }
- }
- }
-
- // allocate space for the rows of the matrix
- for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx)
- tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
- tracerMatrix_->endrowsizes();
-
- // fill the rows with indices. each degree of freedom talks to
- // all of its neighbors. (it also talks to itself since
- // degrees of freedom are sometimes quite egocentric.)
- for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
- typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
- typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
- for (; nIt != nEndIt; ++nIt)
- tracerMatrix_->addindex(dofIdx, *nIt);
- }
- tracerMatrix_->endindices();
-
- const int sizeCartGrid = simulator_.vanguard().cartesianSize();
- cartToGlobal_.resize(sizeCartGrid);
- for (unsigned i = 0; i < numGridDof; ++i) {
- int cartIdx = simulator_.vanguard().cartesianIndex(i);
- cartToGlobal_[cartIdx] = i;
- }
-
- }
-
- /*!
- * \brief Return the number of tracers considered by the tracerModel.
- */
- int numTracers() const
- { return tracerNames_.size(); }
-
- /*!
- * \brief Return the tracer name
- */
- const std::string& tracerName(int tracerIdx) const
- {
- if (numTracers()==0)
- throw std::logic_error("This method should never be called when there are no tracers in the model");
-
- return tracerNames_[tracerIdx];
- }
-
- /*!
- * \brief Return the tracer concentration for tracer index and global DofIdx
- */
- Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const
- {
- if (numTracers()==0)
- return 0.0;
-
- return tracerConcentration_[tracerIdx][globalDofIdx];
+ bool enabled = EWOMS_GET_PARAM(TypeTag, bool, EnableTracerModel);
+ this->doInit(enabled, simulator_.model().numGridDof(),
+ gasPhaseIdx, oilPhaseIdx, waterPhaseIdx);
}
void beginTimeStep()
{
- if (numTracers()==0)
+ if (this->numTracers()==0)
return;
- tracerConcentrationInitial_ = tracerConcentration_;
+ this->tracerConcentrationInitial_ = this->tracerConcentration_;
// compute storageCache
ElementContext elemCtx(simulator_);
@@ -269,10 +117,10 @@ public:
for (; elemIt != elemEndIt; ++ elemIt) {
elemCtx.updateAll(*elemIt);
int globalDofIdx = elemCtx.globalSpaceIndex(0, 0);
- for (int tracerIdx = 0; tracerIdx < numTracers(); ++ tracerIdx){
+ for (int tracerIdx = 0; tracerIdx < this->numTracers(); ++ tracerIdx){
Scalar storageOfTimeIndex1;
computeStorage_(storageOfTimeIndex1, elemCtx, 0, /*timIdx=*/0, tracerIdx);
- storageOfTimeIndex1_[tracerIdx][globalDofIdx] = storageOfTimeIndex1;
+ this->storageOfTimeIndex1_[tracerIdx][globalDofIdx] = storageOfTimeIndex1;
}
}
}
@@ -282,17 +130,17 @@ public:
*/
void endTimeStep()
{
- if (numTracers()==0)
+ if (this->numTracers()==0)
return;
- for (int tracerIdx = 0; tracerIdx < numTracers(); ++ tracerIdx){
+ for (int tracerIdx = 0; tracerIdx < this->numTracers(); ++ tracerIdx){
- TracerVector dx(tracerResidual_.size());
+ typename BaseType::TracerVector dx(this->tracerResidual_.size());
// Newton step (currently the system is linear, converge in one iteration)
for (int iter = 0; iter < 5; ++ iter){
linearize_(tracerIdx);
- linearSolve_(*tracerMatrix_, dx, tracerResidual_);
- tracerConcentration_[tracerIdx] -= dx;
+ this->linearSolve_(*this->tracerMatrix_, dx, this->tracerResidual_);
+ this->tracerConcentration_[tracerIdx] -= dx;
if (dx.two_norm()<1e-2)
break;
@@ -305,7 +153,7 @@ public:
* to the hard disk.
*/
template
- void serialize(Restarter& res OPM_UNUSED)
+ void serialize(Restarter&)
{ /* not implemented */ }
/*!
@@ -315,7 +163,7 @@ public:
* It is the inverse of the serialize() method.
*/
template
- void deserialize(Restarter& res OPM_UNUSED)
+ void deserialize(Restarter&)
{ /* not implemented */ }
protected:
@@ -332,19 +180,19 @@ protected:
const auto& intQuants = elemCtx.intensiveQuantities(scvIdx, timeIdx);
const auto& fs = intQuants.fluidState();
Scalar phaseVolume =
- decay(fs.saturation(tracerPhaseIdx_[tracerIdx]))
- *decay(fs.invB(tracerPhaseIdx_[tracerIdx]))
+ decay(fs.saturation(this->tracerPhaseIdx_[tracerIdx]))
+ *decay(fs.invB(this->tracerPhaseIdx_[tracerIdx]))
*decay(intQuants.porosity());
// avoid singular matrix if no water is present.
phaseVolume = max(phaseVolume, 1e-10);
if (std::is_same::value)
- tracerStorage = phaseVolume * tracerConcentrationInitial_[tracerIdx][globalDofIdx];
+ tracerStorage = phaseVolume * this->tracerConcentrationInitial_[tracerIdx][globalDofIdx];
else
tracerStorage =
phaseVolume
- * variable(tracerConcentration_[tracerIdx][globalDofIdx][0], 0);
+ * variable(this->tracerConcentration_[tracerIdx][globalDofIdx][0], 0);
}
// evaluate the tracerflux over one face
@@ -361,7 +209,7 @@ protected:
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
unsigned inIdx = extQuants.interiorIndex();
- const int tracerPhaseIdx = tracerPhaseIdx_[tracerIdx];
+ const int tracerPhaseIdx = this->tracerPhaseIdx_[tracerIdx];
unsigned upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
int globalUpIdx = elemCtx.globalSpaceIndex(upIdx, timeIdx);
@@ -371,8 +219,8 @@ protected:
Scalar A = scvf.area();
Scalar v = decay(extQuants.volumeFlux(tracerPhaseIdx));
- Scalar b = decay(fs.invB(tracerPhaseIdx_[tracerIdx]));
- Scalar c = tracerConcentration_[tracerIdx][globalUpIdx];
+ Scalar b = decay(fs.invB(this->tracerPhaseIdx_[tracerIdx]));
+ Scalar c = this->tracerConcentration_[tracerIdx][globalUpIdx];
if (inIdx == upIdx)
tracerFlux = A*v*b*variable(c, 0);
@@ -381,41 +229,10 @@ protected:
}
- bool linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
- {
-#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
- Dune::FMatrixPrecision::set_singular_limit(1.e-30);
- Dune::FMatrixPrecision::set_absolute_limit(1.e-30);
-#endif
- x = 0.0;
- Scalar tolerance = 1e-2;
- int maxIter = 100;
-
- int verbosity = 0;
- typedef Dune::BiCGSTABSolver TracerSolver;
- typedef Dune::MatrixAdapter TracerOperator;
- typedef Dune::SeqScalarProduct< TracerVector > TracerScalarProduct ;
- typedef Dune::SeqILU< TracerMatrix, TracerVector, TracerVector > TracerPreconditioner;
-
- TracerOperator tracerOperator(M);
- TracerScalarProduct tracerScalarProduct;
- TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
-
- TracerSolver solver (tracerOperator, tracerScalarProduct,
- tracerPreconditioner, tolerance, maxIter,
- verbosity);
-
- Dune::InverseOperatorResult result;
- solver.apply(x, b, result);
-
- // return the result of the solver
- return result.converged;
- }
-
void linearize_(int tracerIdx)
{
- (*tracerMatrix_) = 0.0;
- tracerResidual_ = 0.0;
+ (*this->tracerMatrix_) = 0.0;
+ this->tracerResidual_ = 0.0;
size_t numGridDof = simulator_.model().numGridDof();
std::vector volumes(numGridDof, 0.0);
@@ -442,13 +259,13 @@ protected:
Scalar storageOfTimeIndex1;
computeStorage_(storageOfTimeIndex0, elemCtx, 0, /*timIdx=*/0, tracerIdx);
if (elemCtx.enableStorageCache())
- storageOfTimeIndex1 = storageOfTimeIndex1_[tracerIdx][I];
+ storageOfTimeIndex1 = this->storageOfTimeIndex1_[tracerIdx][I];
else
computeStorage_(storageOfTimeIndex1, elemCtx, 0, /*timIdx=*/1, tracerIdx);
localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1) * scvVolume/dt;
- tracerResidual_[I][0] += localStorage.value(); //residual + flux
- (*tracerMatrix_)[I][I][0][0] = localStorage.derivative(0);
+ this->tracerResidual_[I][0] += localStorage.value(); //residual + flux
+ (*this->tracerMatrix_)[I][I][0][0] = localStorage.derivative(0);
size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
TracerEvaluation flux;
@@ -456,9 +273,9 @@ protected:
unsigned j = face.exteriorIndex();
unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
computeFlux_(flux, elemCtx, scvfIdx, 0, tracerIdx);
- tracerResidual_[I][0] += flux.value(); //residual + flux
- (*tracerMatrix_)[J][I][0][0] = -flux.derivative(0);
- (*tracerMatrix_)[I][J][0][0] = flux.derivative(0);
+ this->tracerResidual_[I][0] += flux.value(); //residual + flux
+ (*this->tracerMatrix_)[J][I][0][0] = -flux.derivative(0);
+ (*this->tracerMatrix_)[I][J][0][0] = flux.derivative(0);
}
}
@@ -471,7 +288,7 @@ protected:
if (well.getStatus() == Well::Status::SHUT)
continue;
- const double wtracer = well.getTracerProperties().getConcentration(tracerNames_[tracerIdx]);
+ const double wtracer = well.getTracerProperties().getConcentration(this->tracerNames_[tracerIdx]);
std::array cartesianCoordinate;
for (auto& connection : well.getConnections()) {
@@ -482,28 +299,19 @@ protected:
cartesianCoordinate[1] = connection.getJ();
cartesianCoordinate[2] = connection.getK();
const size_t cartIdx = simulator_.vanguard().cartesianIndex(cartesianCoordinate);
- const int I = cartToGlobal_[cartIdx];
- Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, tracerPhaseIdx_[tracerIdx]);
+ const int I = this->cartToGlobal_[cartIdx];
+ Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, this->tracerPhaseIdx_[tracerIdx]);
if (rate > 0)
- tracerResidual_[I][0] -= rate*wtracer;
+ this->tracerResidual_[I][0] -= rate*wtracer;
else if (rate < 0)
- tracerResidual_[I][0] -= rate*tracerConcentration_[tracerIdx][I];
+ this->tracerResidual_[I][0] -= rate*this->tracerConcentration_[tracerIdx][I];
}
}
}
Simulator& simulator_;
-
- std::vector tracerNames_;
- std::vector tracerPhaseIdx_;
- std::vector>> tracerConcentration_;
- std::vector>> tracerConcentrationInitial_;
- TracerMatrix *tracerMatrix_;
- TracerVector tracerResidual_;
- std::vector cartToGlobal_;
- std::vector>> storageOfTimeIndex1_;
-
};
+
} // namespace Opm
#endif