73template<
class Vector,
class Gr
id,
class Matrix>
74std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
75 Dune::OwnerOverlapCopyCommunication<int,int>>>,
76 std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
79 OPM_THROW(std::logic_error,
"Grid not supported for parallel Tracers.");
80 return {
nullptr,
nullptr};
83template<
class Vector,
class Matrix>
84std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
85 Dune::OwnerOverlapCopyCommunication<int,int>>>,
86 std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
89 using TracerOperator = Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
90 Dune::OwnerOverlapCopyCommunication<int,int>>;
92 const auto&
cellComm = grid.cellCommunication();
93 auto op = std::make_unique<TracerOperator>(M,
cellComm);
99template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
102 const EclipseState& eclState,
104 const DofMapper& dofMapper,
105 const std::function<std::array<double,dimWorld>(
int)>
centroids)
106 : gridView_(gridView)
107 , eclState_(eclState)
109 , dofMapper_(dofMapper)
114template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
118 if (freeTracerConcentration_.empty())
121 return freeTracerConcentration_[tracerIdx][globalDofIdx];
124template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
128 if (solTracerConcentration_.empty())
131 return solTracerConcentration_[tracerIdx][globalDofIdx];
134template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
135void GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
136setFreeTracerConcentration(
int tracerIdx,
int globalDofIdx, Scalar value)
138 this->freeTracerConcentration_[tracerIdx][globalDofIdx] = value;
139 this->tracerConcentration_[tracerIdx][globalDofIdx][0] = value;
142template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
143void GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
144setSolTracerConcentration(
int tracerIdx,
int globalDofIdx, Scalar value)
146 this->solTracerConcentration_[tracerIdx][globalDofIdx] = value;
147 this->tracerConcentration_[tracerIdx][globalDofIdx][1] = value;
150template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
151void GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
152setEnableSolTracers(
int tracerIdx,
bool enableSolTracer)
154 this->enableSolTracers_[tracerIdx] = enableSolTracer;
157template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
161 return this->eclState_.tracer().size();
164template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
166fname(
int tracerIdx)
const
168 return this->eclState_.tracer()[tracerIdx].fname();
171template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
172std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
173sname(
int tracerIdx)
const
175 return this->eclState_.tracer()[tracerIdx].sname();
178template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
179std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
180wellfname(
int tracerIdx)
const
182 return this->eclState_.tracer()[tracerIdx].wellfname();
185template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
186std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
187wellsname(
int tracerIdx)
const
189 return this->eclState_.tracer()[tracerIdx].wellsname();
192template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
193Phase GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
194phase(
int tracerIdx)
const
196 return this->eclState_.tracer()[tracerIdx].phase;
199template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
200const std::vector<bool>& GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
201enableSolTracers()
const
203 return this->enableSolTracers_;
206template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
207Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>::
208currentConcentration_(
const Well& eclWell,
const std::string& name)
const
210 return eclWell.getTracerProperties().getConcentration(name);
213template<
class Gr
id,
class Gr
idView,
class DofMapper,
class Stencil,
class Flu
idSystem,
class Scalar>
215name(
int tracerIdx)
const
217 return this->eclState_.tracer()[tracerIdx].
name;
222doInit(
bool rst, std::size_t numGridDof,
223 std::size_t gasPhaseIdx, std::size_t oilPhaseIdx, std::size_t waterPhaseIdx)
225 const auto& tracers = eclState_.tracer();
227 if (tracers.size() == 0)
231 const std::size_t numTracers = tracers.size();
232 enableSolTracers_.resize(numTracers);
233 tracerConcentration_.resize(numTracers);
234 freeTracerConcentration_.resize(numTracers);
235 solTracerConcentration_.resize(numTracers);
238 tracerPhaseIdx_.resize(numTracers);
239 for (std::size_t tracerIdx = 0; tracerIdx < numTracers; tracerIdx++) {
240 const auto& tracer = tracers[tracerIdx];
242 if (tracer.phase == Phase::WATER)
243 tracerPhaseIdx_[tracerIdx] = waterPhaseIdx;
244 else if (tracer.phase == Phase::OIL)
245 tracerPhaseIdx_[tracerIdx] = oilPhaseIdx;
246 else if (tracer.phase == Phase::GAS)
247 tracerPhaseIdx_[tracerIdx] = gasPhaseIdx;
249 tracerConcentration_[tracerIdx].resize(numGridDof);
250 freeTracerConcentration_[tracerIdx].resize(numGridDof);
251 solTracerConcentration_[tracerIdx].resize(numGridDof);
258 if (tracer.free_concentration.has_value()){
259 const auto& free_concentration = tracer.free_concentration.value();
260 int tblkDatasize = free_concentration.size();
261 if (tblkDatasize < cartMapper_.cartesianSize()){
262 throw std::runtime_error(
"Wrong size of TBLKF for" + tracer.name);
264 for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
265 int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
266 tracerConcentration_[tracerIdx][globalDofIdx][0] = free_concentration[cartDofIdx];
267 freeTracerConcentration_[tracerIdx][globalDofIdx] = free_concentration[cartDofIdx];
271 else if (tracer.free_tvdp.has_value()) {
272 const auto& free_tvdp = tracer.free_tvdp.value();
273 for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
274 tracerConcentration_[tracerIdx][globalDofIdx][0] =
275 free_tvdp.evaluate(
"TRACER_CONCENTRATION",
276 centroids_(globalDofIdx)[2]);
277 freeTracerConcentration_[tracerIdx][globalDofIdx] =
278 free_tvdp.evaluate(
"TRACER_CONCENTRATION",
279 centroids_(globalDofIdx)[2]);
283 OpmLog::warning(fmt::format(
"No TBLKF or TVDPF given for free tracer {}. "
284 "Initial values set to zero. ", tracer.name));
285 for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
286 tracerConcentration_[tracerIdx][globalDofIdx][0] = 0.0;
287 freeTracerConcentration_[tracerIdx][globalDofIdx] = 0.0;
292 if (tracer.phase != Phase::WATER &&
293 ((tracer.phase == Phase::GAS && FluidSystem::enableDissolvedGas()) ||
294 (tracer.phase == Phase::OIL && FluidSystem::enableVaporizedOil()))) {
296 if (tracer.solution_concentration.has_value()){
297 enableSolTracers_[tracerIdx] =
true;
298 const auto& solution_concentration = tracer.solution_concentration.value();
299 int tblkDatasize = solution_concentration.size();
300 if (tblkDatasize < cartMapper_.cartesianSize()){
301 throw std::runtime_error(
"Wrong size of TBLKS for" + tracer.name);
303 for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
304 int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
305 tracerConcentration_[tracerIdx][globalDofIdx][1] = solution_concentration[cartDofIdx];
306 solTracerConcentration_[tracerIdx][globalDofIdx] = solution_concentration[cartDofIdx];
310 else if (tracer.solution_tvdp.has_value()) {
311 enableSolTracers_[tracerIdx] =
true;
312 const auto& solution_tvdp = tracer.solution_tvdp.value();
313 for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
314 tracerConcentration_[tracerIdx][globalDofIdx][1] =
315 solution_tvdp.evaluate(
"TRACER_CONCENTRATION",
316 centroids_(globalDofIdx)[2]);
317 solTracerConcentration_[tracerIdx][globalDofIdx] =
318 solution_tvdp.evaluate(
"TRACER_CONCENTRATION",
319 centroids_(globalDofIdx)[2]);
324 enableSolTracers_[tracerIdx] =
false;
325 OpmLog::warning(fmt::format(
"No TBLKS or TVDPS given for solution tracer {}. "
326 "Initial values set to zero. ", tracer.name));
327 for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
328 tracerConcentration_[tracerIdx][globalDofIdx][1] = 0.0;
329 solTracerConcentration_[tracerIdx][globalDofIdx] = 0.0;
335 enableSolTracers_[tracerIdx] =
false;
336 for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
337 tracerConcentration_[tracerIdx][globalDofIdx][1] = 0.0;
338 solTracerConcentration_[tracerIdx][globalDofIdx] = 0.0;
344 tracerMatrix_ = std::make_unique<TracerMatrix>(numGridDof, numGridDof, TracerMatrix::random);
347 using NeighborSet = std::set<unsigned>;
348 std::vector<NeighborSet> neighbors(numGridDof);
350 Stencil stencil(gridView_, dofMapper_);
351 for (
const auto& elem : elements(gridView_)) {
352 stencil.update(elem);
354 for (
unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
355 unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
357 for (
unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
358 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
359 neighbors[myIdx].insert(neighborIdx);
365 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
366 tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
368 tracerMatrix_->endrowsizes();
373 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
374 typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
375 typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
376 for (; nIt != nEndIt; ++nIt) {
377 tracerMatrix_->addindex(dofIdx, *nIt);
380 tracerMatrix_->endindices();