My Project
AquiferAnalytical.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_AQUIFERANALYTICAL_HEADER_INCLUDED
23 #define OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
24 
25 #include <opm/simulators/aquifers/AquiferInterface.hpp>
26 
27 #include <opm/common/utility/numeric/linearInterpolation.hpp>
28 #include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
29 
30 #include <opm/output/data/Aquifer.hpp>
31 
32 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
33 
34 #include <opm/material/common/MathToolbox.hpp>
35 #include <opm/material/densead/Evaluation.hpp>
36 #include <opm/material/densead/Math.hpp>
37 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
38 
39 #include <algorithm>
40 #include <cmath>
41 #include <cstddef>
42 #include <limits>
43 #include <numeric>
44 #include <unordered_map>
45 #include <vector>
46 
47 namespace Opm
48 {
49 template <typename TypeTag>
50 class AquiferAnalytical : public AquiferInterface<TypeTag>
51 {
52 public:
53  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
54  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
55  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
56  using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
57  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
58  using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
59  using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
60 
61  enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
62  enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
63  enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
64  enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
65  enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
66 
67  static constexpr int numEq = BlackoilIndices::numEq;
68  using Scalar = double;
69 
70  using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
71 
72  using FluidState = BlackOilFluidState<Eval,
73  FluidSystem,
74  enableTemperature,
75  enableEnergy,
76  BlackoilIndices::gasEnabled,
77  enableEvaporation,
78  enableBrine,
79  enableSaltPrecipitation,
80  BlackoilIndices::numPhases>;
81 
82  // Constructor
83  AquiferAnalytical(int aqID,
84  const std::vector<Aquancon::AquancCell>& connections,
85  const Simulator& ebosSimulator)
86  : AquiferInterface<TypeTag>(aqID, ebosSimulator)
87  , connections_(connections)
88  {
89  }
90 
91  // Destructor
92  virtual ~AquiferAnalytical()
93  {
94  }
95 
96  void initFromRestart(const data::Aquifers& aquiferSoln) override
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() override
110  {
111  initQuantities();
112  }
113 
114  void beginTimeStep() override
115  {
116  ElementContext elemCtx(this->ebos_simulator_);
117  OPM_BEGIN_PARALLEL_TRY_CATCH();
118 
119  for (const auto& elem : elements(this->ebos_simulator_.gridView())) {
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(this->phaseIdx_()));
130  }
131 
132  OPM_END_PARALLEL_TRY_CATCH("AquiferAnalytical::beginTimeStep() failed: ",
133  this->ebos_simulator_.vanguard().grid().comm());
134  }
135 
136  void addToSource(RateVector& rates,
137  const unsigned cellIdx,
138  const unsigned timeIdx) override
139  {
140  const auto& model = this->ebos_simulator_.model();
141 
142  const int idx = this->cellToConnectionIdx_[cellIdx];
143  if (idx < 0)
144  return;
145 
146  const auto* intQuantsPtr = model.cachedIntensiveQuantities(cellIdx, timeIdx);
147  if (intQuantsPtr == nullptr) {
148  throw std::logic_error("Invalid intensive quantities cache detected in AquiferAnalytical::addToSource()");
149  }
150 
151  // This is the pressure at td + dt
152  this->updateCellPressure(this->pressure_current_, idx, *intQuantsPtr);
153  this->calculateInflowRate(idx, this->ebos_simulator_);
154 
155  rates[BlackoilIndices::conti0EqIdx + compIdx_()]
156  += this->Qai_[idx] / model.dofTotalVolume(cellIdx);
157 
158  if constexpr (enableEnergy) {
159  auto fs = intQuantsPtr->fluidState();
160  if (this->Ta0_.has_value() && this->Qai_[idx] > 0)
161  {
162  fs.setTemperature(this->Ta0_.value());
163  typedef typename std::decay<decltype(fs)>::type::Scalar FsScalar;
164  typename FluidSystem::template ParameterCache<FsScalar> paramCache;
165  const unsigned pvtRegionIdx = intQuantsPtr->pvtRegionIndex();
166  paramCache.setRegionIndex(pvtRegionIdx);
167  paramCache.setMaxOilSat(this->ebos_simulator_.problem().maxOilSaturation(cellIdx));
168  paramCache.updatePhase(fs, this->phaseIdx_());
169  const auto& h = FluidSystem::enthalpy(fs, paramCache, this->phaseIdx_());
170  fs.setEnthalpy(this->phaseIdx_(), h);
171  }
172  rates[BlackoilIndices::contiEnergyEqIdx]
173  += this->Qai_[idx] *fs.enthalpy(this->phaseIdx_()) * FluidSystem::referenceDensity( this->phaseIdx_(), intQuantsPtr->pvtRegionIndex()) / model.dofTotalVolume(cellIdx);
174 
175  }
176  }
177 
178  std::size_t size() const
179  {
180  return this->connections_.size();
181  }
182 
183 protected:
184  virtual void assignRestartData(const data::AquiferData& xaq) = 0;
185  virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0;
186  virtual void calculateAquiferCondition() = 0;
187  virtual void calculateAquiferConstants() = 0;
188  virtual Scalar aquiferDepth() const = 0;
189 
190  Scalar gravity_() const
191  {
192  return this->ebos_simulator_.problem().gravity()[2];
193  }
194 
195  int compIdx_() const
196  {
197  if (this->co2store_())
198  return FluidSystem::oilCompIdx;
199 
200  return FluidSystem::waterCompIdx;
201  }
202 
203 
204  void initQuantities()
205  {
206  // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
207  if (!this->solution_set_from_restart_) {
208  W_flux_ = Scalar{0};
209  }
210 
211  // We next get our connections to the aquifer and initialize these quantities using the initialize_connections
212  // function
213  initializeConnections();
214  calculateAquiferCondition();
215  calculateAquiferConstants();
216 
217  pressure_previous_.resize(this->connections_.size(), Scalar{0});
218  pressure_current_.resize(this->connections_.size(), Scalar{0});
219  Qai_.resize(this->connections_.size(), Scalar{0});
220  }
221 
222  void updateCellPressure(std::vector<Eval>& pressure_water,
223  const int idx,
224  const IntensiveQuantities& intQuants)
225  {
226  const auto& fs = intQuants.fluidState();
227  pressure_water.at(idx) = fs.pressure(this->phaseIdx_());
228  }
229 
230  void updateCellPressure(std::vector<Scalar>& pressure_water,
231  const int idx,
232  const IntensiveQuantities& intQuants)
233  {
234  const auto& fs = intQuants.fluidState();
235  pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value();
236  }
237 
238  void initializeConnections()
239  {
240  this->cell_depth_.resize(this->size(), this->aquiferDepth());
241  this->alphai_.resize(this->size(), 1.0);
242  this->faceArea_connected_.resize(this->size(), Scalar{0});
243 
244  // Translate the C face tag into the enum used by opm-parser's TransMult class
245  FaceDir::DirEnum faceDirection;
246 
247  bool has_active_connection_on_proc = false;
248 
249  // denom_face_areas is the sum of the areas connected to an aquifer
250  Scalar denom_face_areas{0};
251  this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
252  const auto& gridView = this->ebos_simulator_.vanguard().gridView();
253  for (std::size_t idx = 0; idx < this->size(); ++idx) {
254  const auto global_index = this->connections_[idx].global_index;
255  const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
256  auto elemIt = gridView.template begin</*codim=*/ 0>();
257  if (cell_index > 0)
258  std::advance(elemIt, cell_index);
259 
260  //the global_index is not part of this grid
261  if ( cell_index < 0 || elemIt->partitionType() != Dune::InteriorEntity)
262  continue;
263 
264  has_active_connection_on_proc = true;
265 
266  this->cellToConnectionIdx_[cell_index] = idx;
267  this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
268  }
269  // get areas for all connections
270  ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
271  for (const auto& elem : elements(gridView)) {
272  unsigned cell_index = elemMapper.index(elem);
273  int idx = this->cellToConnectionIdx_[cell_index];
274 
275  // only deal with connections given by the aquifer
276  if( idx < 0)
277  continue;
278 
279  for (const auto& intersection : intersections(gridView, elem)) {
280  // only deal with grid boundaries
281  if (!intersection.boundary())
282  continue;
283 
284  int insideFaceIdx = intersection.indexInInside();
285  switch (insideFaceIdx) {
286  case 0:
287  faceDirection = FaceDir::XMinus;
288  break;
289  case 1:
290  faceDirection = FaceDir::XPlus;
291  break;
292  case 2:
293  faceDirection = FaceDir::YMinus;
294  break;
295  case 3:
296  faceDirection = FaceDir::YPlus;
297  break;
298  case 4:
299  faceDirection = FaceDir::ZMinus;
300  break;
301  case 5:
302  faceDirection = FaceDir::ZPlus;
303  break;
304  default:
305  OPM_THROW(std::logic_error,
306  "Internal error in initialization of aquifer.");
307  }
308 
309 
310  if (faceDirection == this->connections_[idx].face_dir) {
311  this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
312  break;
313  }
314  }
315  denom_face_areas += this->faceArea_connected_.at(idx);
316  }
317 
318  const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
319  comm.sum(&denom_face_areas, 1);
320  const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
321  for (std::size_t idx = 0; idx < this->size(); ++idx) {
322  // Protect against division by zero NaNs.
323  this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
324  ? Scalar{0}
325  : this->faceArea_connected_.at(idx) / denom_face_areas;
326  }
327 
328  if (this->solution_set_from_restart_) {
329  this->rescaleProducedVolume(has_active_connection_on_proc);
330  }
331  }
332 
333  void rescaleProducedVolume(const bool has_active_connection_on_proc)
334  {
335  // Needed in parallel restart to approximate influence of aquifer
336  // being "owned" by a subset of the parallel processes. If the
337  // aquifer is fully owned by a single process--i.e., if all cells
338  // connecting to the aquifer are on a single process--then this_area
339  // is tot_area on that process and zero elsewhere.
340 
341  const auto this_area = has_active_connection_on_proc
342  ? std::accumulate(this->alphai_.begin(),
343  this->alphai_.end(),
344  Scalar{0})
345  : Scalar{0};
346 
347  const auto tot_area = this->ebos_simulator_.vanguard()
348  .grid().comm().sum(this_area);
349 
350  this->W_flux_ *= this_area / tot_area;
351  }
352 
353  // This function is for calculating the aquifer properties from equilibrium state with the reservoir
354  Scalar calculateReservoirEquilibrium()
355  {
356  // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
357  std::vector<Scalar> pw_aquifer;
358  Scalar water_pressure_reservoir;
359 
360  ElementContext elemCtx(this->ebos_simulator_);
361  const auto& gridView = this->ebos_simulator_.gridView();
362  for (const auto& elem : elements(gridView)) {
363  elemCtx.updatePrimaryStencil(elem);
364 
365  const auto cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
366  const auto idx = this->cellToConnectionIdx_[cellIdx];
367  if (idx < 0)
368  continue;
369 
370  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
371  const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
372  const auto& fs = iq0.fluidState();
373 
374  water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
375  const auto water_density = fs.density(this->phaseIdx_());
376 
377  const auto gdz =
378  this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
379 
380  pw_aquifer.push_back(this->alphai_[idx] *
381  (water_pressure_reservoir - water_density.value()*gdz));
382  }
383 
384  // We take the average of the calculated equilibrium pressures.
385  const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
386 
387  Scalar vals[2];
388  vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
389  vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), Scalar{0});
390 
391  comm.sum(vals, 2);
392 
393  return vals[1] / vals[0];
394  }
395 
396  const std::vector<Aquancon::AquancCell> connections_;
397 
398  // Grid variables
399  std::vector<Scalar> faceArea_connected_;
400  std::vector<int> cellToConnectionIdx_;
401 
402  // Quantities at each grid id
403  std::vector<Scalar> cell_depth_;
404  std::vector<Scalar> pressure_previous_;
405  std::vector<Eval> pressure_current_;
406  std::vector<Eval> Qai_;
407  std::vector<Scalar> alphai_;
408 
409  Scalar Tc_{}; // Time constant
410  Scalar pa0_{}; // initial aquifer pressure
411  std::optional<Scalar> Ta0_{}; // initial aquifer temperature
412  Scalar rhow_{};
413 
414  Eval W_flux_;
415 
416  bool solution_set_from_restart_ {false};
417 };
418 
419 } // namespace Opm
420 
421 #endif
Definition: AquiferAnalytical.hpp:51
Definition: AquiferInterface.hpp:32
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27