My Project
AquiferInterface.hpp
1 /*
2  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3  Copyright 2017 Statoil ASA.
4  Copyright 2017 IRIS
5 
6  This file is part of the Open Porous Media project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 #ifndef OPM_AQUIFERINTERFACE_HEADER_INCLUDED
23 #define OPM_AQUIFERINTERFACE_HEADER_INCLUDED
24 
25 #include <opm/common/utility/numeric/linearInterpolation.hpp>
26 #include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
27 
28 #include <opm/output/data/Aquifer.hpp>
29 
30 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
31 
32 #include <opm/material/common/MathToolbox.hpp>
33 #include <opm/material/densead/Evaluation.hpp>
34 #include <opm/material/densead/Math.hpp>
35 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
36 
37 #include <algorithm>
38 #include <cmath>
39 #include <cstddef>
40 #include <limits>
41 #include <numeric>
42 #include <unordered_map>
43 #include <vector>
44 
45 namespace Opm
46 {
47 template <typename TypeTag>
49 {
50 public:
51  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
52  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
53  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
54  using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
55  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
56  using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
57  using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
58 
59  enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
60  enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
61  enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
62  enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
63  enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
64 
65  static const int numEq = BlackoilIndices::numEq;
66  typedef double Scalar;
67 
68  typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
69 
70  typedef BlackOilFluidState<Eval,
71  FluidSystem,
72  enableTemperature,
73  enableEnergy,
74  BlackoilIndices::gasEnabled,
75  enableEvaporation,
76  enableBrine,
77  enableSaltPrecipitation,
78  BlackoilIndices::numPhases>
79  FluidState;
80 
81  // Constructor
82  AquiferInterface(int aqID,
83  const std::vector<Aquancon::AquancCell>& connections,
84  const Simulator& ebosSimulator)
85  : aquiferID_(aqID)
86  , connections_(connections)
87  , ebos_simulator_(ebosSimulator)
88  {
89  }
90 
91  // Destructor
92  virtual ~AquiferInterface()
93  {
94  }
95 
96  void initFromRestart(const data::Aquifers& aquiferSoln)
97  {
98  auto xaqPos = aquiferSoln.find(this->aquiferID());
99  if (xaqPos == aquiferSoln.end())
100  return;
101 
102  this->assignRestartData(xaqPos->second);
103 
104  this->W_flux_ = xaqPos->second.volume;
105  this->pa0_ = xaqPos->second.initPressure;
106  this->solution_set_from_restart_ = true;
107  }
108 
109  void initialSolutionApplied()
110  {
111  initQuantities();
112  }
113 
114  void beginTimeStep()
115  {
116  ElementContext elemCtx(ebos_simulator_);
117  auto elemIt = ebos_simulator_.gridView().template begin<0>();
118  const auto& elemEndIt = ebos_simulator_.gridView().template end<0>();
119  OPM_BEGIN_PARALLEL_TRY_CATCH();
120 
121  for (; elemIt != elemEndIt; ++elemIt) {
122  const auto& elem = *elemIt;
123 
124  elemCtx.updatePrimaryStencil(elem);
125 
126  const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
127  const int idx = cellToConnectionIdx_[cellIdx];
128  if (idx < 0)
129  continue;
130 
131  elemCtx.updateIntensiveQuantities(0);
132  const auto& iq = elemCtx.intensiveQuantities(0, 0);
133  pressure_previous_[idx] = getValue(iq.fluidState().pressure(phaseIdx_()));
134  }
135 
136  OPM_END_PARALLEL_TRY_CATCH("AquiferInterface::beginTimeStep() failed: ", ebos_simulator_.vanguard().grid().comm());
137  }
138 
139  template <class Context>
140  void addToSource(RateVector& rates,
141  const Context& context,
142  const unsigned spaceIdx,
143  const unsigned timeIdx)
144  {
145  const unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
146 
147  const int idx = this->cellToConnectionIdx_[cellIdx];
148  if (idx < 0)
149  return;
150 
151  // We are dereferencing the value of IntensiveQuantities because
152  // cachedIntensiveQuantities return a const pointer to
153  // IntensiveQuantities of that particular cell_id
154  const auto& intQuants = context.intensiveQuantities(spaceIdx, timeIdx);
155 
156  // This is the pressure at td + dt
157  this->updateCellPressure(this->pressure_current_, idx, intQuants);
158  this->calculateInflowRate(idx, context.simulator());
159 
160  rates[BlackoilIndices::conti0EqIdx + compIdx_()]
161  += this->Qai_[idx] / context.dofVolume(spaceIdx, timeIdx);
162  }
163 
164  std::size_t size() const {
165  return this->connections_.size();
166  }
167 
168  int aquiferID() const { return this->aquiferID_; }
169 
170 protected:
171  inline Scalar gravity_() const
172  {
173  return ebos_simulator_.problem().gravity()[2];
174  }
175 
176  inline bool co2store_() const
177  {
178  return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
179  }
180 
181  inline int phaseIdx_() const
182  {
183  if(co2store_())
184  return FluidSystem::oilPhaseIdx;
185 
186  return FluidSystem::waterPhaseIdx;
187  }
188 
189  inline int compIdx_() const
190  {
191  if(co2store_())
192  return FluidSystem::oilCompIdx;
193 
194  return FluidSystem::waterCompIdx;
195  }
196 
197 
198  inline void initQuantities()
199  {
200  // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
201  if (!this->solution_set_from_restart_) {
202  W_flux_ = Scalar{0};
203  }
204 
205  // We next get our connections to the aquifer and initialize these quantities using the initialize_connections
206  // function
207  initializeConnections();
208  calculateAquiferCondition();
209  calculateAquiferConstants();
210 
211  pressure_previous_.resize(this->connections_.size(), Scalar{0});
212  pressure_current_.resize(this->connections_.size(), Scalar{0});
213  Qai_.resize(this->connections_.size(), Scalar{0});
214  }
215 
216  inline void
217  updateCellPressure(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
218  {
219  const auto& fs = intQuants.fluidState();
220  pressure_water.at(idx) = fs.pressure(phaseIdx_());
221  }
222 
223  inline void
224  updateCellPressure(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
225  {
226  const auto& fs = intQuants.fluidState();
227  pressure_water.at(idx) = fs.pressure(phaseIdx_()).value();
228  }
229 
230  virtual void endTimeStep() = 0;
231 
232  const int aquiferID_{};
233  const std::vector<Aquancon::AquancCell> connections_;
234  const Simulator& ebos_simulator_;
235 
236  // Grid variables
237  std::vector<Scalar> faceArea_connected_;
238  std::vector<int> cellToConnectionIdx_;
239 
240  // Quantities at each grid id
241  std::vector<Scalar> cell_depth_;
242  std::vector<Scalar> pressure_previous_;
243  std::vector<Eval> pressure_current_;
244  std::vector<Eval> Qai_;
245  std::vector<Scalar> alphai_;
246 
247  Scalar Tc_{}; // Time constant
248  Scalar pa0_{}; // initial aquifer pressure
249  Scalar rhow_{};
250 
251  Eval W_flux_;
252 
253  bool solution_set_from_restart_ {false};
254 
255  void initializeConnections()
256  {
257  this->cell_depth_.resize(this->size(), this->aquiferDepth());
258  this->alphai_.resize(this->size(), 1.0);
259  this->faceArea_connected_.resize(this->size(), Scalar{0});
260 
261  // Translate the C face tag into the enum used by opm-parser's TransMult class
262  FaceDir::DirEnum faceDirection;
263 
264  bool has_active_connection_on_proc = false;
265 
266  // denom_face_areas is the sum of the areas connected to an aquifer
267  Scalar denom_face_areas{0};
268  this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
269  const auto& gridView = this->ebos_simulator_.vanguard().gridView();
270  for (std::size_t idx = 0; idx < this->size(); ++idx) {
271  const auto global_index = this->connections_[idx].global_index;
272  const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
273  auto elemIt = gridView.template begin</*codim=*/ 0>();
274  if (cell_index > 0)
275  std::advance(elemIt, cell_index);
276 
277  //the global_index is not part of this grid
278  if ( cell_index < 0 || elemIt->partitionType() != Dune::InteriorEntity)
279  continue;
280 
281  has_active_connection_on_proc = true;
282 
283  this->cellToConnectionIdx_[cell_index] = idx;
284  this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
285  }
286  // get areas for all connections
287  ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
288  auto elemIt = gridView.template begin</*codim=*/ 0>();
289  const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
290  for (; elemIt != elemEndIt; ++elemIt) {
291  const auto& elem = *elemIt;
292  unsigned cell_index = elemMapper.index(elem);
293  int idx = this->cellToConnectionIdx_[cell_index];
294 
295  // only deal with connections given by the aquifer
296  if( idx < 0)
297  continue;
298 
299  auto isIt = gridView.ibegin(elem);
300  const auto& isEndIt = gridView.iend(elem);
301  for (; isIt != isEndIt; ++ isIt) {
302  // store intersection, this might be costly
303  const auto& intersection = *isIt;
304 
305  // only deal with grid boundaries
306  if (!intersection.boundary())
307  continue;
308 
309  int insideFaceIdx = intersection.indexInInside();
310  switch (insideFaceIdx) {
311  case 0:
312  faceDirection = FaceDir::XMinus;
313  break;
314  case 1:
315  faceDirection = FaceDir::XPlus;
316  break;
317  case 2:
318  faceDirection = FaceDir::YMinus;
319  break;
320  case 3:
321  faceDirection = FaceDir::YPlus;
322  break;
323  case 4:
324  faceDirection = FaceDir::ZMinus;
325  break;
326  case 5:
327  faceDirection = FaceDir::ZPlus;
328  break;
329  default:
330  OPM_THROW(std::logic_error,
331  "Internal error in initialization of aquifer.");
332  }
333 
334 
335  if (faceDirection == this->connections_[idx].face_dir) {
336  this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
337  break;
338  }
339  }
340  denom_face_areas += this->faceArea_connected_.at(idx);
341  }
342 
343  const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
344  comm.sum(&denom_face_areas, 1);
345  const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
346  for (std::size_t idx = 0; idx < this->size(); ++idx) {
347  // Protect against division by zero NaNs.
348  this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
349  ? Scalar{0}
350  : this->faceArea_connected_.at(idx) / denom_face_areas;
351  }
352 
353  if (this->solution_set_from_restart_) {
354  this->rescaleProducedVolume(has_active_connection_on_proc);
355  }
356  }
357 
358  void rescaleProducedVolume(const bool has_active_connection_on_proc)
359  {
360  // Needed in parallel restart to approximate influence of aquifer
361  // being "owned" by a subset of the parallel processes. If the
362  // aquifer is fully owned by a single process--i.e., if all cells
363  // connecting to the aquifer are on a single process--then this_area
364  // is tot_area on that process and zero elsewhere.
365 
366  const auto this_area = has_active_connection_on_proc
367  ? std::accumulate(this->alphai_.begin(),
368  this->alphai_.end(),
369  Scalar{0})
370  : Scalar{0};
371 
372  const auto tot_area = this->ebos_simulator_.vanguard()
373  .grid().comm().sum(this_area);
374 
375  this->W_flux_ *= this_area / tot_area;
376  }
377 
378  virtual void assignRestartData(const data::AquiferData& xaq) = 0;
379 
380  virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0;
381 
382  virtual void calculateAquiferCondition() = 0;
383 
384  virtual void calculateAquiferConstants() = 0;
385 
386  virtual Scalar aquiferDepth() const = 0;
387 
388  // This function is for calculating the aquifer properties from equilibrium state with the reservoir
389  virtual Scalar calculateReservoirEquilibrium()
390  {
391  // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
392  std::vector<Scalar> pw_aquifer;
393  Scalar water_pressure_reservoir;
394 
395  ElementContext elemCtx(this->ebos_simulator_);
396  const auto& gridView = this->ebos_simulator_.gridView();
397  auto elemIt = gridView.template begin</*codim=*/0>();
398  const auto& elemEndIt = gridView.template end</*codim=*/0>();
399  for (; elemIt != elemEndIt; ++elemIt) {
400  const auto& elem = *elemIt;
401  elemCtx.updatePrimaryStencil(elem);
402 
403  const auto cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
404  const auto idx = this->cellToConnectionIdx_[cellIdx];
405  if (idx < 0)
406  continue;
407 
408  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
409  const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
410  const auto& fs = iq0.fluidState();
411 
412  water_pressure_reservoir = fs.pressure(phaseIdx_()).value();
413  const auto water_density = fs.density(phaseIdx_());
414 
415  const auto gdz =
416  this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
417 
418  pw_aquifer.push_back(this->alphai_[idx] *
419  (water_pressure_reservoir - water_density.value()*gdz));
420  }
421 
422  // We take the average of the calculated equilibrium pressures.
423  const auto& comm = ebos_simulator_.vanguard().grid().comm();
424 
425  Scalar vals[2];
426  vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
427  vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), Scalar{0});
428 
429  comm.sum(vals, 2);
430 
431  return vals[1] / vals[0];
432  }
433 
434  // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
435 };
436 } // namespace Opm
437 #endif
Definition: AquiferInterface.hpp:49
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27