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