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/parser/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 auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
54 
55  static const int numEq = BlackoilIndices::numEq;
56 
57  using Eval = DenseAd::Evaluation<double, numEq>;
58  using Toolbox = MathToolbox<Eval>;
59 
60  // Constructor
61  AquiferNumerical(const SingleNumericalAquifer& aquifer,
62  const std::unordered_map<int, int>& cartesian_to_compressed,
63  const Simulator& ebos_simulator,
64  const int* global_cell)
65  : id_ (aquifer.id())
66  , ebos_simulator_ (ebos_simulator)
67  , flux_rate_ (0.0)
68  , cumulative_flux_(0.0)
69  , global_cell_ (global_cell)
70  , init_pressure_ (aquifer.numCells(), 0.0)
71  {
72  this->cell_to_aquifer_cell_idx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
73 
74  auto aquifer_on_process = false;
75  for (std::size_t idx = 0; idx < aquifer.numCells(); ++idx) {
76  const auto* cell = aquifer.getCellPrt(idx);
77 
78  // Due to parallelisation, the cell might not exist in the current process
79  auto search = cartesian_to_compressed.find(cell->global_index);
80  if (search != cartesian_to_compressed.end()) {
81  this->cell_to_aquifer_cell_idx_[search->second] = idx;
82  aquifer_on_process = true;
83  }
84  }
85 
86  if (aquifer_on_process) {
87  this->checkConnectsToReservoir();
88  }
89  }
90 
91  void initFromRestart(const data::Aquifers& aquiferSoln)
92  {
93  auto xaqPos = aquiferSoln.find(this->aquiferID());
94  if (xaqPos == aquiferSoln.end())
95  return;
96 
97  if (this->connects_to_reservoir_) {
98  this->cumulative_flux_ = xaqPos->second.volume;
99  }
100 
101  if (const auto* aqData = xaqPos->second.typeData.template get<data::AquiferType::Numerical>();
102  aqData != nullptr)
103  {
104  this->init_pressure_ = aqData->initPressure;
105  }
106 
107  this->solution_set_from_restart_ = true;
108  }
109 
110  void endTimeStep()
111  {
112  this->pressure_ = this->calculateAquiferPressure();
113  this->flux_rate_ = this->calculateAquiferFluxRate();
114  this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
115  }
116 
117  data::AquiferData aquiferData() const
118  {
119  data::AquiferData data;
120  data.aquiferID = this->aquiferID();
121  data.pressure = this->pressure_;
122  data.fluxRate = this->flux_rate_;
123  data.volume = this->cumulative_flux_;
124 
125  auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
126  aquNum->initPressure = this->init_pressure_;
127 
128  return data;
129  }
130 
131  void initialSolutionApplied()
132  {
133  if (this->solution_set_from_restart_) {
134  return;
135  }
136 
137  this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
138  this->flux_rate_ = 0.;
139  this->cumulative_flux_ = 0.;
140  }
141 
142  int aquiferID() const
143  {
144  return static_cast<int>(this->id_);
145  }
146 
147 private:
148  const std::size_t id_;
149  const Simulator& ebos_simulator_;
150  double flux_rate_; // aquifer influx rate
151  double cumulative_flux_; // cumulative aquifer influx
152  const int* global_cell_; // mapping to global index
153  std::vector<double> init_pressure_{};
154  double pressure_; // aquifer pressure
155  bool solution_set_from_restart_ {false};
156  bool connects_to_reservoir_ {false};
157 
158  // TODO: maybe unordered_map can also do the work to save memory?
159  std::vector<int> cell_to_aquifer_cell_idx_;
160 
161  void checkConnectsToReservoir()
162  {
163  ElementContext elem_ctx(this->ebos_simulator_);
164  auto elemIt = std::find_if(this->ebos_simulator_.gridView().template begin</*codim=*/0>(),
165  this->ebos_simulator_.gridView().template end</*codim=*/0>(),
166  [&elem_ctx, this](const auto& elem) -> bool
167  {
168  elem_ctx.updateStencil(elem);
169 
170  const auto cell_index = elem_ctx
171  .globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
172 
173  return this->cell_to_aquifer_cell_idx_[cell_index] == 0;
174  });
175 
176  assert ((elemIt != this->ebos_simulator_.gridView().template end</*codim=*/0>())
177  && "Internal error locating numerical aquifer's connecting cell");
178 
179  this->connects_to_reservoir_ =
180  elemIt->partitionType() == Dune::InteriorEntity;
181  }
182 
183  double calculateAquiferPressure() const
184  {
185  auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
186  return this->calculateAquiferPressure(capture);
187  }
188 
189  double calculateAquiferPressure(std::vector<double>& cell_pressure) const
190  {
191  double sum_pressure_watervolume = 0.;
192  double sum_watervolume = 0.;
193 
194  ElementContext elem_ctx(this->ebos_simulator_);
195  const auto& gridView = this->ebos_simulator_.gridView();
196  auto elemIt = gridView.template begin</*codim=*/0>();
197  const auto& elemEndIt = gridView.template end</*codim=*/0>();
198  OPM_BEGIN_PARALLEL_TRY_CATCH();
199 
200  for (; elemIt != elemEndIt; ++elemIt) {
201  const auto& elem = *elemIt;
202  if (elem.partitionType() != Dune::InteriorEntity) {
203  continue;
204  }
205  elem_ctx.updatePrimaryStencil(elem);
206 
207  const size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
208  const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
209  if (idx < 0) {
210  continue;
211  }
212 
213  elem_ctx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
214  const auto& iq0 = elem_ctx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
215  const auto& fs = iq0.fluidState();
216 
217  // TODO: the porosity of the cells are still wrong for numerical aquifer cells
218  // Because the dofVolume still based on the grid information.
219  // The pore volume is correct. Extra efforts will be done to get sensible porosity value here later.
220  const double water_saturation = fs.saturation(waterPhaseIdx).value();
221  const double porosity = iq0.porosity().value();
222  const double volume = elem_ctx.dofTotalVolume(0, 0);
223  // TODO: not sure we should use water pressure here
224  const double water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
225  const double water_volume = volume * porosity * water_saturation;
226  sum_pressure_watervolume += water_volume * water_pressure_reservoir;
227  sum_watervolume += water_volume;
228 
229  cell_pressure[idx] = water_pressure_reservoir;
230  }
231  OPM_END_PARALLEL_TRY_CATCH("AquiferNumerical::calculateAquiferPressure() failed: ", this->ebos_simulator_.vanguard().grid().comm());
232  const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
233  comm.sum(&sum_pressure_watervolume, 1);
234  comm.sum(&sum_watervolume, 1);
235 
236  // Ensure all processes have same notion of the aquifer cells' pressure values.
237  comm.sum(cell_pressure.data(), cell_pressure.size());
238 
239  return sum_pressure_watervolume / sum_watervolume;
240  }
241 
242  double calculateAquiferFluxRate() const
243  {
244  double aquifer_flux = 0.0;
245 
246  if (! this->connects_to_reservoir_) {
247  return aquifer_flux;
248  }
249 
250  ElementContext elem_ctx(this->ebos_simulator_);
251  const auto& gridView = this->ebos_simulator_.gridView();
252  auto elemIt = gridView.template begin</*codim=*/0>();
253  const auto& elemEndIt = gridView.template end</*codim=*/0>();
254  for (; elemIt != elemEndIt; ++elemIt) {
255  const auto& elem = *elemIt;
256  if (elem.partitionType() != Dune::InteriorEntity) {
257  continue;
258  }
259  // elem_ctx.updatePrimaryStencil(elem);
260  elem_ctx.updateStencil(elem);
261 
262  const std::size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
263  const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
264  // we only need the first aquifer cell
265  if (idx != 0) {
266  continue;
267  }
268  elem_ctx.updateAllIntensiveQuantities();
269  elem_ctx.updateAllExtensiveQuantities();
270 
271  const std::size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0);
272  // const auto &problem = elem_ctx.problem();
273  const auto& stencil = elem_ctx.stencil(0);
274  // const auto& inQuants = elem_ctx.intensiveQuantities(0, /*timeIdx*/ 0);
275 
276  for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
277  const auto& face = stencil.interiorFace(face_idx);
278  // dof index
279  const std::size_t i = face.interiorIndex();
280  const std::size_t j = face.exteriorIndex();
281  // compressed index
282  // const size_t I = stencil.globalSpaceIndex(i);
283  const std::size_t J = stencil.globalSpaceIndex(j);
284 
285  assert(stencil.globalSpaceIndex(i) == cell_index);
286 
287  // we do not consider the flux within aquifer cells
288  // we only need the flux to the connections
289  if (this->cell_to_aquifer_cell_idx_[J] > 0) {
290  continue;
291  }
292  const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
293  const double water_flux = Toolbox::value(exQuants.volumeFlux(waterPhaseIdx));
294 
295  const std::size_t up_id = water_flux >= 0.0 ? i : j;
296  const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
297  const double invB = Toolbox::value(intQuantsIn.fluidState().invB(waterPhaseIdx));
298  const double face_area = face.area();
299  aquifer_flux += water_flux * invB * face_area;
300  }
301 
302  // we only need to handle the first aquifer cell, we can exit loop here
303  break;
304  }
305 
306  return aquifer_flux;
307  }
308 };
309 } // namespace Opm
310 #endif
Definition: AquiferNumerical.hpp:41
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:26