My Project
AquiferNumerical.hpp
1 /*
2  Copyright (C) 2020 Equinor ASA
3  Copyright (C) 2020 SINTEF Digital
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
22 #define OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
23 
24 #include <opm/output/data/Aquifer.hpp>
25 
26 #include <opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
27 
28 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
29 
30 #include <algorithm>
31 #include <cassert>
32 #include <cstddef>
33 #include <unordered_map>
34 #include <utility>
35 #include <vector>
36 
37 namespace Opm
38 {
39 template <typename TypeTag>
41 {
42 public:
43  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
44  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
45  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
46  using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
47 
48  using GridView = GetPropType<TypeTag, Properties::GridView>;
49  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
50 
51  enum { dimWorld = GridView::dimensionworld };
52 
53  static const int numEq = BlackoilIndices::numEq;
54 
55  using Eval = DenseAd::Evaluation<double, numEq>;
56  using Toolbox = MathToolbox<Eval>;
57 
58  // Constructor
59  AquiferNumerical(const SingleNumericalAquifer& aquifer,
60  const std::unordered_map<int, int>& cartesian_to_compressed,
61  const Simulator& ebos_simulator,
62  const int* global_cell)
63  : id_ (aquifer.id())
64  , ebos_simulator_ (ebos_simulator)
65  , flux_rate_ (0.0)
66  , cumulative_flux_(0.0)
67  , global_cell_ (global_cell)
68  , init_pressure_ (aquifer.numCells(), 0.0)
69  {
70  this->cell_to_aquifer_cell_idx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
71 
72  auto aquifer_on_process = false;
73  for (std::size_t idx = 0; idx < aquifer.numCells(); ++idx) {
74  const auto* cell = aquifer.getCellPrt(idx);
75 
76  // Due to parallelisation, the cell might not exist in the current process
77  auto search = cartesian_to_compressed.find(cell->global_index);
78  if (search != cartesian_to_compressed.end()) {
79  this->cell_to_aquifer_cell_idx_[search->second] = idx;
80  aquifer_on_process = true;
81  }
82  }
83 
84  if (aquifer_on_process) {
85  this->checkConnectsToReservoir();
86  }
87  }
88 
89  void initFromRestart(const data::Aquifers& aquiferSoln)
90  {
91  auto xaqPos = aquiferSoln.find(this->aquiferID());
92  if (xaqPos == aquiferSoln.end())
93  return;
94 
95  if (this->connects_to_reservoir_) {
96  this->cumulative_flux_ = xaqPos->second.volume;
97  }
98 
99  if (const auto* aqData = xaqPos->second.typeData.template get<data::AquiferType::Numerical>();
100  aqData != nullptr)
101  {
102  this->init_pressure_ = aqData->initPressure;
103  }
104 
105  this->solution_set_from_restart_ = true;
106  }
107 
108  void endTimeStep()
109  {
110  this->pressure_ = this->calculateAquiferPressure();
111  this->flux_rate_ = this->calculateAquiferFluxRate();
112  this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
113  }
114 
115  data::AquiferData aquiferData() const
116  {
117  data::AquiferData data;
118  data.aquiferID = this->aquiferID();
119  data.pressure = this->pressure_;
120  data.fluxRate = this->flux_rate_;
121  data.volume = this->cumulative_flux_;
122 
123  auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
124  aquNum->initPressure = this->init_pressure_;
125 
126  return data;
127  }
128 
129  void initialSolutionApplied()
130  {
131  if (this->solution_set_from_restart_) {
132  return;
133  }
134 
135  this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
136  this->flux_rate_ = 0.;
137  this->cumulative_flux_ = 0.;
138  }
139 
140  int aquiferID() const
141  {
142  return static_cast<int>(this->id_);
143  }
144 
145 private:
146  const std::size_t id_;
147  const Simulator& ebos_simulator_;
148  double flux_rate_; // aquifer influx rate
149  double cumulative_flux_; // cumulative aquifer influx
150  const int* global_cell_; // mapping to global index
151  std::vector<double> init_pressure_{};
152  double pressure_; // aquifer pressure
153  bool solution_set_from_restart_ {false};
154  bool connects_to_reservoir_ {false};
155 
156  // TODO: maybe unordered_map can also do the work to save memory?
157  std::vector<int> cell_to_aquifer_cell_idx_;
158 
159  inline bool co2store_() const
160  {
161  return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
162  }
163 
164  inline int phaseIdx_() const
165  {
166  if(co2store_())
167  return FluidSystem::oilPhaseIdx;
168 
169  return FluidSystem::waterPhaseIdx;
170  }
171 
172  void checkConnectsToReservoir()
173  {
174  ElementContext elem_ctx(this->ebos_simulator_);
175  auto elemIt = std::find_if(this->ebos_simulator_.gridView().template begin</*codim=*/0>(),
176  this->ebos_simulator_.gridView().template end</*codim=*/0>(),
177  [&elem_ctx, this](const auto& elem) -> bool
178  {
179  elem_ctx.updateStencil(elem);
180 
181  const auto cell_index = elem_ctx
182  .globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
183 
184  return this->cell_to_aquifer_cell_idx_[cell_index] == 0;
185  });
186 
187  assert ((elemIt != this->ebos_simulator_.gridView().template end</*codim=*/0>())
188  && "Internal error locating numerical aquifer's connecting cell");
189 
190  this->connects_to_reservoir_ =
191  elemIt->partitionType() == Dune::InteriorEntity;
192  }
193 
194  double calculateAquiferPressure() const
195  {
196  auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
197  return this->calculateAquiferPressure(capture);
198  }
199 
200  double calculateAquiferPressure(std::vector<double>& cell_pressure) const
201  {
202  double sum_pressure_watervolume = 0.;
203  double sum_watervolume = 0.;
204 
205  ElementContext elem_ctx(this->ebos_simulator_);
206  const auto& gridView = this->ebos_simulator_.gridView();
207  auto elemIt = gridView.template begin</*codim=*/0>();
208  const auto& elemEndIt = gridView.template end</*codim=*/0>();
209  OPM_BEGIN_PARALLEL_TRY_CATCH();
210 
211  for (; elemIt != elemEndIt; ++elemIt) {
212  const auto& elem = *elemIt;
213  if (elem.partitionType() != Dune::InteriorEntity) {
214  continue;
215  }
216  elem_ctx.updatePrimaryStencil(elem);
217 
218  const size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
219  const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
220  if (idx < 0) {
221  continue;
222  }
223 
224  elem_ctx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
225  const auto& iq0 = elem_ctx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
226  const auto& fs = iq0.fluidState();
227 
228  // TODO: the porosity of the cells are still wrong for numerical aquifer cells
229  // Because the dofVolume still based on the grid information.
230  // The pore volume is correct. Extra efforts will be done to get sensible porosity value here later.
231  const double water_saturation = fs.saturation(this->phaseIdx_()).value();
232  const double porosity = iq0.porosity().value();
233  const double volume = elem_ctx.dofTotalVolume(0, 0);
234  // TODO: not sure we should use water pressure here
235  const double water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
236  const double water_volume = volume * porosity * water_saturation;
237  sum_pressure_watervolume += water_volume * water_pressure_reservoir;
238  sum_watervolume += water_volume;
239 
240  cell_pressure[idx] = water_pressure_reservoir;
241  }
242  OPM_END_PARALLEL_TRY_CATCH("AquiferNumerical::calculateAquiferPressure() failed: ", this->ebos_simulator_.vanguard().grid().comm());
243  const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
244  comm.sum(&sum_pressure_watervolume, 1);
245  comm.sum(&sum_watervolume, 1);
246 
247  // Ensure all processes have same notion of the aquifer cells' pressure values.
248  comm.sum(cell_pressure.data(), cell_pressure.size());
249 
250  return sum_pressure_watervolume / sum_watervolume;
251  }
252 
253  double calculateAquiferFluxRate() const
254  {
255  double aquifer_flux = 0.0;
256 
257  if (! this->connects_to_reservoir_) {
258  return aquifer_flux;
259  }
260 
261  ElementContext elem_ctx(this->ebos_simulator_);
262  const auto& gridView = this->ebos_simulator_.gridView();
263  auto elemIt = gridView.template begin</*codim=*/0>();
264  const auto& elemEndIt = gridView.template end</*codim=*/0>();
265  for (; elemIt != elemEndIt; ++elemIt) {
266  const auto& elem = *elemIt;
267  if (elem.partitionType() != Dune::InteriorEntity) {
268  continue;
269  }
270  // elem_ctx.updatePrimaryStencil(elem);
271  elem_ctx.updateStencil(elem);
272 
273  const std::size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
274  const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
275  // we only need the first aquifer cell
276  if (idx != 0) {
277  continue;
278  }
279  elem_ctx.updateAllIntensiveQuantities();
280  elem_ctx.updateAllExtensiveQuantities();
281 
282  const std::size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0);
283  // const auto &problem = elem_ctx.problem();
284  const auto& stencil = elem_ctx.stencil(0);
285  // const auto& inQuants = elem_ctx.intensiveQuantities(0, /*timeIdx*/ 0);
286 
287  for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
288  const auto& face = stencil.interiorFace(face_idx);
289  // dof index
290  const std::size_t i = face.interiorIndex();
291  const std::size_t j = face.exteriorIndex();
292  // compressed index
293  // const size_t I = stencil.globalSpaceIndex(i);
294  const std::size_t J = stencil.globalSpaceIndex(j);
295 
296  assert(stencil.globalSpaceIndex(i) == cell_index);
297 
298  // we do not consider the flux within aquifer cells
299  // we only need the flux to the connections
300  if (this->cell_to_aquifer_cell_idx_[J] > 0) {
301  continue;
302  }
303  const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
304  const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
305 
306  const std::size_t up_id = water_flux >= 0.0 ? i : j;
307  const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
308  const double invB = Toolbox::value(intQuantsIn.fluidState().invB(this->phaseIdx_()));
309  const double face_area = face.area();
310  aquifer_flux += water_flux * invB * face_area;
311  }
312 
313  // we only need to handle the first aquifer cell, we can exit loop here
314  break;
315  }
316 
317  return aquifer_flux;
318  }
319 };
320 } // namespace Opm
321 #endif
Definition: AquiferNumerical.hpp:41
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27