28#ifndef OPM_GENERIC_TRACER_MODEL_IMPL_HPP
29#define OPM_GENERIC_TRACER_MODEL_IMPL_HPP
31#include <dune/istl/operators.hh>
32#include <dune/istl/solvers.hh>
33#include <dune/istl/schwarz.hh>
34#include <dune/istl/preconditioners.hh>
35#include <dune/istl/schwarz.hh>
37#include <opm/common/OpmLog/OpmLog.hpp>
39#include <opm/grid/CpGrid.hpp>
41#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
42#include <opm/input/eclipse/EclipseState/Phase.hpp>
43#include <opm/input/eclipse/EclipseState/Tables/TracerVdTable.hpp>
44#include <opm/input/eclipse/Schedule/Well/Well.hpp>
45#include <opm/input/eclipse/Schedule/Well/WellTracerProperties.hpp>
47#include <opm/models/discretization/ecfv/ecfvstencil.hh>
50#include <opm/simulators/linalg/ilufirstelement.hh>
54#include <fmt/format.h>
66template<
class M,
class V>
69 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
74template<
class Vector,
class Gr
id,
class Matrix>
75std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
76 Dune::OwnerOverlapCopyCommunication<int,int>>>,
77 std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
80 OPM_THROW(std::logic_error,
"Grid not supported for parallel Tracers.");
81 return {
nullptr,
nullptr};
84template<
class Vector,
class Matrix>
85std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
86 Dune::OwnerOverlapCopyCommunication<int,int>>>,
87 std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
90 using TracerOperator = Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
91 Dune::OwnerOverlapCopyCommunication<int,int>>;
93 const auto& cellComm = grid.cellCommunication();
94 auto op = std::make_unique<TracerOperator>(M, cellComm);
95 auto dummyWeights = [](){
return Vector();};
96 return {std::move(op), std::make_unique<TracerSolver>(*op, cellComm, prm, dummyWeights, 0)};
100template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Scalar>
103 const EclipseState& eclState,
105 const DofMapper& dofMapper,
106 const std::function<std::array<double,dimWorld>(
int)> centroids)
107 : gridView_(gridView)
108 , eclState_(eclState)
109 , cartMapper_(cartMapper)
110 , dofMapper_(dofMapper)
111 , centroids_(centroids)
115template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Scalar>
119 if (tracerConcentration_.empty())
122 return tracerConcentration_[tracerIdx][globalDofIdx];
125template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Scalar>
129 this->tracerConcentration_[tracerIdx][globalDofIdx] = value;
132template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Scalar>
136 return this->eclState_.tracer().size();
139template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Scalar>
141fname(
int tracerIdx)
const
143 return this->eclState_.tracer()[tracerIdx].fname();
146template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Scalar>
150 return eclWell.getTracerProperties().getConcentration(name);
153template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Scalar>
155name(
int tracerIdx)
const
157 return this->eclState_.tracer()[tracerIdx].name;
160template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Scalar>
162doInit(
bool rst, std::size_t numGridDof,
163 std::size_t gasPhaseIdx, std::size_t oilPhaseIdx, std::size_t waterPhaseIdx)
165 const auto& tracers = eclState_.tracer();
167 if (tracers.size() == 0)
171 const std::size_t numTracers = tracers.size();
172 tracerConcentration_.resize(numTracers);
173 storageOfTimeIndex1_.resize(numTracers);
176 tracerPhaseIdx_.resize(numTracers);
177 for (std::size_t tracerIdx = 0; tracerIdx < numTracers; tracerIdx++) {
178 const auto& tracer = tracers[tracerIdx];
180 if (tracer.phase == Phase::WATER)
181 tracerPhaseIdx_[tracerIdx] = waterPhaseIdx;
182 else if (tracer.phase == Phase::OIL)
183 tracerPhaseIdx_[tracerIdx] = oilPhaseIdx;
184 else if (tracer.phase == Phase::GAS)
185 tracerPhaseIdx_[tracerIdx] = gasPhaseIdx;
187 tracerConcentration_[tracerIdx].resize(numGridDof);
188 storageOfTimeIndex1_[tracerIdx].resize(numGridDof);
195 if (tracer.free_concentration.has_value()){
196 const auto& free_concentration = tracer.free_concentration.value();
197 int tblkDatasize = free_concentration.size();
198 if (tblkDatasize < cartMapper_.cartesianSize()){
199 throw std::runtime_error(
"Wrong size of TBLK for" + tracer.name);
201 for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
202 int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
203 tracerConcentration_[tracerIdx][globalDofIdx] = free_concentration[cartDofIdx];
207 else if (tracer.free_tvdp.has_value()) {
208 const auto& free_tvdp = tracer.free_tvdp.value();
209 for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
210 tracerConcentration_[tracerIdx][globalDofIdx] =
211 free_tvdp.evaluate(
"TRACER_CONCENTRATION",
212 centroids_(globalDofIdx)[2]);
215 throw std::logic_error(fmt::format(
"Can not initialize tracer: {}", tracer.name));
219 tracerMatrix_ = std::make_unique<TracerMatrix>(numGridDof, numGridDof, TracerMatrix::random);
222 using NeighborSet = std::set<unsigned>;
223 std::vector<NeighborSet> neighbors(numGridDof);
225 Stencil stencil(gridView_, dofMapper_);
226 for (
const auto& elem : elements(gridView_)) {
227 stencil.update(elem);
229 for (
unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
230 unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
232 for (
unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
233 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
234 neighbors[myIdx].insert(neighborIdx);
240 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx)
241 tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
242 tracerMatrix_->endrowsizes();
247 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
248 typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
249 typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
250 for (; nIt != nEndIt; ++nIt)
251 tracerMatrix_->addindex(dofIdx, *nIt);
253 tracerMatrix_->endindices();
256template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Scalar>
261 Scalar tolerance = 1e-2;
266 prm.
put(
"maxiter", maxIter);
267 prm.
put(
"tol", tolerance);
268 prm.
put(
"verbosity", verbosity);
269 prm.
put(
"solver", std::string(
"bicgstab"));
270 prm.
put(
"preconditioner.type", std::string(
"ParOverILU0"));
273 if(gridView_.grid().comm().size() > 1)
275 auto [tracerOperator, solver] =
276 createParallelFlexibleSolver<TracerVector>(gridView_.grid(), M, prm);
277 (void) tracerOperator;
280 solver->apply(x, b, result);
283 return result.converged;
288 using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
289 using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
290 using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
291 using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
293 TracerOperator tracerOperator(M);
294 TracerScalarProduct tracerScalarProduct;
295 TracerPreconditioner tracerPreconditioner(M, 0, 1);
297 TracerSolver solver (tracerOperator, tracerScalarProduct,
298 tracerPreconditioner, tolerance, maxIter,
302 solver.apply(x, b, result);
305 return result.converged;
311template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Scalar>
315 Scalar tolerance = 1e-2;
320 prm.
put(
"maxiter", maxIter);
321 prm.
put(
"tol", tolerance);
322 prm.
put(
"verbosity", verbosity);
323 prm.
put(
"solver", std::string(
"bicgstab"));
324 prm.
put(
"preconditioner.type", std::string(
"ParOverILU0"));
327 if(gridView_.grid().comm().size() > 1)
329 auto [tracerOperator, solver] =
330 createParallelFlexibleSolver<TracerVector>(gridView_.grid(), M, prm);
331 (void) tracerOperator;
332 bool converged =
true;
333 for (std::size_t nrhs = 0; nrhs < b.size(); ++nrhs) {
336 solver->apply(x[nrhs], b[nrhs], result);
337 converged = (converged && result.converged);
344 using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
345 using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
346 using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
347 using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
349 TracerOperator tracerOperator(M);
350 TracerScalarProduct tracerScalarProduct;
351 TracerPreconditioner tracerPreconditioner(M, 0, 1);
353 TracerSolver solver (tracerOperator, tracerScalarProduct,
354 tracerPreconditioner, tolerance, maxIter,
357 bool converged =
true;
358 for (std::size_t nrhs = 0; nrhs < b.size(); ++nrhs) {
361 solver.apply(x[nrhs], b[nrhs], result);
362 converged = (converged && result.converged);
Definition: CollectDataOnIORank.hpp:49
Definition: FlexibleSolver.hpp:45
Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const
Return the tracer concentration for tracer index and global DofIdx.
Definition: GenericTracerModel_impl.hpp:117
bool linearSolve_(const TracerMatrix &M, TracerVector &x, TracerVector &b)
Definition: GenericTracerModel_impl.hpp:258
std::string fname(int tracerIdx) const
Definition: GenericTracerModel_impl.hpp:141
void doInit(bool rst, std::size_t numGridDof, std::size_t gasPhaseIdx, std::size_t oilPhaseIdx, std::size_t waterPhaseIdx)
Initialize all internal data structures needed by the tracer module.
Definition: GenericTracerModel_impl.hpp:162
const std::string & name(int tracerIdx) const
Return the tracer name.
Definition: GenericTracerModel_impl.hpp:155
Dune::BlockVector< Dune::FieldVector< Scalar, 1 > > TracerVector
Definition: GenericTracerModel.hpp:56
void setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value)
Definition: GenericTracerModel_impl.hpp:127
bool linearSolveBatchwise_(const TracerMatrix &M, std::vector< TracerVector > &x, std::vector< TracerVector > &b)
Definition: GenericTracerModel_impl.hpp:313
GenericTracerModel(const GridView &gridView, const EclipseState &eclState, const CartesianIndexMapper &cartMapper, const DofMapper &dofMapper, const std::function< std::array< double, dimWorld >(int)> centroids)
Definition: GenericTracerModel_impl.hpp:102
int numTracers() const
Return the number of tracers considered by the tracerModel.
Definition: GenericTracerModel_impl.hpp:134
Dune::BCRSMatrix< Opm::MatrixBlock< Scalar, 1, 1 > > TracerMatrix
Definition: GenericTracerModel.hpp:55
double currentConcentration_(const Well &eclWell, const std::string &name) const
Definition: GenericTracerModel_impl.hpp:148
Definition: PropertyTree.hpp:37
void put(const std::string &key, const T &data)
Definition: BlackoilPhases.hpp:27
Dune::InverseOperatorResult InverseOperatorResult
Definition: BdaBridge.hpp:32
std::tuple< std::unique_ptr< Dune::OverlappingSchwarzOperator< Matrix, Vector, Vector, Dune::OwnerOverlapCopyCommunication< int, int > > >, std::unique_ptr< typename TracerSolverSelector< Matrix, Vector >::type > > createParallelFlexibleSolver(const Grid &, const Matrix &, const PropertyTree &)
Definition: GenericTracerModel_impl.hpp:78
Definition: GenericTracerModel_impl.hpp:68
Dune::OverlappingSchwarzOperator< M, V, V, Comm > TracerOperator
Definition: GenericTracerModel_impl.hpp:70
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: GenericTracerModel_impl.hpp:69