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/aquifers/AquiferInterface.hpp>
29#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
30
31#include <algorithm>
32#include <cassert>
33#include <cstddef>
34#include <unordered_map>
35#include <utility>
36#include <vector>
37
38namespace Opm
39{
40template <typename TypeTag>
41class AquiferNumerical : public AquiferInterface<TypeTag>
42{
43public:
44 using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
45 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
46 using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
47 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
48 using GridView = GetPropType<TypeTag, Properties::GridView>;
49 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
50 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
51 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
52
53 enum { dimWorld = GridView::dimensionworld };
54 enum { numPhases = FluidSystem::numPhases };
55 static constexpr int numEq = BlackoilIndices::numEq;
56
57 using Eval = DenseAd::Evaluation<double, numEq>;
58 using Toolbox = MathToolbox<Eval>;
59
60 using typename AquiferInterface<TypeTag>::RateVector;
61
62 // Constructor
63 AquiferNumerical(const SingleNumericalAquifer& aquifer,
64 const Simulator& ebos_simulator)
65 : AquiferInterface<TypeTag>(aquifer.id(), ebos_simulator)
66 , flux_rate_ (0.0)
67 , cumulative_flux_(0.0)
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 const int compressed_idx = ebos_simulator.vanguard().compressedIndexForInterior(cell->global_index);
78 if (compressed_idx >= 0) {
79 this->cell_to_aquifer_cell_idx_[compressed_idx] = 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) override
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 beginTimeStep() override {}
109 void addToSource(RateVector&, const unsigned, const unsigned) override {}
110
111 void endTimeStep() override
112 {
113 this->pressure_ = this->calculateAquiferPressure();
114 this->flux_rate_ = this->calculateAquiferFluxRate();
115 this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
116 }
117
118 data::AquiferData aquiferData() const override
119 {
120 data::AquiferData data;
121 data.aquiferID = this->aquiferID();
122 data.pressure = this->pressure_;
123 data.fluxRate = this->flux_rate_;
124 data.volume = this->cumulative_flux_;
125
126 auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
127 aquNum->initPressure = this->init_pressure_;
128
129 return data;
130 }
131
132 void initialSolutionApplied() override
133 {
134 if (this->solution_set_from_restart_) {
135 return;
136 }
137
138 this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
139 this->flux_rate_ = 0.;
140 this->cumulative_flux_ = 0.;
141 }
142
143private:
144 void checkConnectsToReservoir()
145 {
146 ElementContext elem_ctx(this->ebos_simulator_);
147 auto elemIt = std::find_if(this->ebos_simulator_.gridView().template begin</*codim=*/0>(),
148 this->ebos_simulator_.gridView().template end</*codim=*/0>(),
149 [&elem_ctx, this](const auto& elem) -> bool
150 {
151 elem_ctx.updateStencil(elem);
152
153 const auto cell_index = elem_ctx
154 .globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
155
156 return this->cell_to_aquifer_cell_idx_[cell_index] == 0;
157 });
158
159 assert ((elemIt != this->ebos_simulator_.gridView().template end</*codim=*/0>())
160 && "Internal error locating numerical aquifer's connecting cell");
161
162 this->connects_to_reservoir_ =
163 elemIt->partitionType() == Dune::InteriorEntity;
164 }
165
166 double calculateAquiferPressure() const
167 {
168 auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
169 return this->calculateAquiferPressure(capture);
170 }
171
172 double calculateAquiferPressure(std::vector<double>& cell_pressure) const
173 {
174 double sum_pressure_watervolume = 0.;
175 double sum_watervolume = 0.;
176
177 ElementContext elem_ctx(this->ebos_simulator_);
178 const auto& gridView = this->ebos_simulator_.gridView();
179 OPM_BEGIN_PARALLEL_TRY_CATCH();
180
181 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
182 elem_ctx.updatePrimaryStencil(elem);
183
184 const size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
185 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
186 if (idx < 0) {
187 continue;
188 }
189
190 elem_ctx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
191 const auto& iq0 = elem_ctx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
192 const auto& fs = iq0.fluidState();
193
194 // TODO: the porosity of the cells are still wrong for numerical aquifer cells
195 // Because the dofVolume still based on the grid information.
196 // The pore volume is correct. Extra efforts will be done to get sensible porosity value here later.
197 const double water_saturation = fs.saturation(this->phaseIdx_()).value();
198 const double porosity = iq0.porosity().value();
199 const double volume = elem_ctx.dofTotalVolume(0, 0);
200 // TODO: not sure we should use water pressure here
201 const double water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
202 const double water_volume = volume * porosity * water_saturation;
203 sum_pressure_watervolume += water_volume * water_pressure_reservoir;
204 sum_watervolume += water_volume;
205
206 cell_pressure[idx] = water_pressure_reservoir;
207 }
208 OPM_END_PARALLEL_TRY_CATCH("AquiferNumerical::calculateAquiferPressure() failed: ", this->ebos_simulator_.vanguard().grid().comm());
209 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
210 comm.sum(&sum_pressure_watervolume, 1);
211 comm.sum(&sum_watervolume, 1);
212
213 // Ensure all processes have same notion of the aquifer cells' pressure values.
214 comm.sum(cell_pressure.data(), cell_pressure.size());
215
216 return sum_pressure_watervolume / sum_watervolume;
217 }
218
219 template <class ElemCtx>
220 double getWaterFlux(const ElemCtx& elem_ctx, unsigned face_idx) const
221 {
222 const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
223 const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
224 return water_flux;
225 }
226
227 double calculateAquiferFluxRate() const
228 {
229 double aquifer_flux = 0.0;
230
231 if (! this->connects_to_reservoir_) {
232 return aquifer_flux;
233 }
234
235 ElementContext elem_ctx(this->ebos_simulator_);
236 const auto& gridView = this->ebos_simulator_.gridView();
237 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
238 // elem_ctx.updatePrimaryStencil(elem);
239 elem_ctx.updateStencil(elem);
240
241 const std::size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
242 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
243 // we only need the first aquifer cell
244 if (idx != 0) {
245 continue;
246 }
247
248 const std::size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0);
249 // const auto &problem = elem_ctx.problem();
250 const auto& stencil = elem_ctx.stencil(0);
251 // const auto& inQuants = elem_ctx.intensiveQuantities(0, /*timeIdx*/ 0);
252
253 for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
254 const auto& face = stencil.interiorFace(face_idx);
255 // dof index
256 const std::size_t i = face.interiorIndex();
257 const std::size_t j = face.exteriorIndex();
258 // compressed index
259 // const size_t I = stencil.globalSpaceIndex(i);
260 const std::size_t J = stencil.globalSpaceIndex(j);
261
262 assert(stencil.globalSpaceIndex(i) == cell_index);
263
264 // we do not consider the flux within aquifer cells
265 // we only need the flux to the connections
266 if (this->cell_to_aquifer_cell_idx_[J] > 0) {
267 continue;
268 }
269 elem_ctx.updateAllIntensiveQuantities();
270 elem_ctx.updateAllExtensiveQuantities();
271
272 const double water_flux = getWaterFlux(elem_ctx,face_idx);
273 const std::size_t up_id = water_flux >= 0.0 ? i : j;
274 const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
275 const double invB = Toolbox::value(intQuantsIn.fluidState().invB(this->phaseIdx_()));
276 const double face_area = face.area();
277 aquifer_flux += water_flux * invB * face_area;
278 }
279
280 // we only need to handle the first aquifer cell, we can exit loop here
281 break;
282 }
283
284 return aquifer_flux;
285 }
286
287 double flux_rate_; // aquifer influx rate
288 double cumulative_flux_; // cumulative aquifer influx
289 std::vector<double> init_pressure_{};
290 double pressure_; // aquifer pressure
291 bool solution_set_from_restart_ {false};
292 bool connects_to_reservoir_ {false};
293
294 // TODO: maybe unordered_map can also do the work to save memory?
295 std::vector<int> cell_to_aquifer_cell_idx_;
296};
297
298} // namespace Opm
299
300#endif
Definition: AquiferInterface.hpp:32
Definition: AquiferNumerical.hpp:42
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27