My Project
Loading...
Searching...
No Matches
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/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
25
26#include <opm/material/common/MathToolbox.hpp>
27#include <opm/material/densead/Evaluation.hpp>
28
29#include <opm/output/data/Aquifer.hpp>
30
31#include <opm/simulators/aquifers/AquiferInterface.hpp>
32#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
33
34#include <algorithm>
35#include <cassert>
36#include <cstddef>
37#include <unordered_map>
38#include <utility>
39#include <vector>
40
41namespace Opm
42{
43template <typename TypeTag>
44class AquiferNumerical : public AquiferInterface<TypeTag>
45{
46public:
55
56 enum { dimWorld = GridView::dimensionworld };
57 enum { numPhases = FluidSystem::numPhases };
58 static constexpr int numEq = BlackoilIndices::numEq;
59
60 using Eval = DenseAd::Evaluation<double, numEq>;
62
64
65 // Constructor
69 , flux_rate_ (0.0)
70 , cumulative_flux_(0.0)
71 , init_pressure_ (aquifer.numCells(), 0.0)
72 {
73 this->cell_to_aquifer_cell_idx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
74
75 auto aquifer_on_process = false;
76 for (std::size_t idx = 0; idx < aquifer.numCells(); ++idx) {
77 const auto* cell = aquifer.getCellPrt(idx);
78
79 // Due to parallelisation, the cell might not exist in the current process
80 const int compressed_idx = ebos_simulator.vanguard().compressedIndexForInterior(cell->global_index);
81 if (compressed_idx >= 0) {
82 this->cell_to_aquifer_cell_idx_[compressed_idx] = idx;
83 aquifer_on_process = true;
84 }
85 }
86
88 this->checkConnectsToReservoir();
89 }
90 }
91
92 static AquiferNumerical serializationTestObject(const Simulator& ebos_simulator)
93 {
95 result.flux_rate_ = 1.0;
96 result.cumulative_flux_ = 2.0;
97 result.init_pressure_ = {3.0, 4.0};
98 result.pressure_ = 5.0;
99
100 return result;
101 }
102
103 void initFromRestart(const data::Aquifers& aquiferSoln) override
104 {
105 auto xaqPos = aquiferSoln.find(this->aquiferID());
106 if (xaqPos == aquiferSoln.end())
107 return;
108
109 if (this->connects_to_reservoir_) {
110 this->cumulative_flux_ = xaqPos->second.volume;
111 }
112
113 if (const auto* aqData = xaqPos->second.typeData.template get<data::AquiferType::Numerical>();
114 aqData != nullptr)
115 {
116 this->init_pressure_ = aqData->initPressure;
117 }
118
119 this->solution_set_from_restart_ = true;
120 }
121
122 void beginTimeStep() override {}
123 void addToSource(RateVector&, const unsigned, const unsigned) override {}
124
125 void endTimeStep() override
126 {
127 this->pressure_ = this->calculateAquiferPressure();
128 this->flux_rate_ = this->calculateAquiferFluxRate();
129 this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
130 }
131
132 data::AquiferData aquiferData() const override
133 {
134 data::AquiferData data;
135 data.aquiferID = this->aquiferID();
136 data.pressure = this->pressure_;
137 data.fluxRate = this->flux_rate_;
138 data.volume = this->cumulative_flux_;
139
140 auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
141 aquNum->initPressure = this->init_pressure_;
142
143 return data;
144 }
145
146 void initialSolutionApplied() override
147 {
148 if (this->solution_set_from_restart_) {
149 return;
150 }
151
152 this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
153 this->flux_rate_ = 0.;
154 this->cumulative_flux_ = 0.;
155 }
156
157 void computeFaceAreaFraction(const std::vector<double>& /*total_face_area*/) override
158 {}
159
160 double totalFaceArea() const override
161 {
162 return 1.0;
163 }
164
165 template<class Serializer>
166 void serializeOp(Serializer& serializer)
167 {
168 serializer(flux_rate_);
169 serializer(cumulative_flux_);
170 serializer(init_pressure_);
171 serializer(pressure_);
172 }
173
174 bool operator==(const AquiferNumerical& rhs) const
175 {
176 return this->flux_rate_ == rhs.flux_rate_ &&
177 this->cumulative_flux_ == rhs.cumulative_flux_ &&
178 this->init_pressure_ == rhs.init_pressure_ &&
179 this->pressure_ == rhs.pressure_;
180 }
181
182 double cumulativeFlux() const
183 {
184 return this->cumulative_flux_;
185 }
186
187private:
188 void checkConnectsToReservoir()
189 {
190 ElementContext elem_ctx(this->ebos_simulator_);
191 auto elemIt = std::find_if(this->ebos_simulator_.gridView().template begin</*codim=*/0>(),
192 this->ebos_simulator_.gridView().template end</*codim=*/0>(),
193 [&elem_ctx, this](const auto& elem) -> bool
194 {
195 elem_ctx.updateStencil(elem);
196
197 const auto cell_index = elem_ctx
198 .globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
199
200 return this->cell_to_aquifer_cell_idx_[cell_index] == 0;
201 });
202
203 assert ((elemIt != this->ebos_simulator_.gridView().template end</*codim=*/0>())
204 && "Internal error locating numerical aquifer's connecting cell");
205
206 this->connects_to_reservoir_ =
207 elemIt->partitionType() == Dune::InteriorEntity;
208 }
209
210 double calculateAquiferPressure() const
211 {
212 auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
213 return this->calculateAquiferPressure(capture);
214 }
215
216 double calculateAquiferPressure(std::vector<double>& cell_pressure) const
217 {
218 double sum_pressure_watervolume = 0.;
219 double sum_watervolume = 0.;
220
221 ElementContext elem_ctx(this->ebos_simulator_);
222 const auto& gridView = this->ebos_simulator_.gridView();
223 OPM_BEGIN_PARALLEL_TRY_CATCH();
224
225 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
226 elem_ctx.updatePrimaryStencil(elem);
227
228 const std::size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
229 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
230 if (idx < 0) {
231 continue;
232 }
233
234 elem_ctx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
235 const auto& iq0 = elem_ctx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
236 const auto& fs = iq0.fluidState();
237
238 // TODO: the porosity of the cells are still wrong for numerical aquifer cells
239 // Because the dofVolume still based on the grid information.
240 // The pore volume is correct. Extra efforts will be done to get sensible porosity value here later.
241 const double water_saturation = fs.saturation(this->phaseIdx_()).value();
242 const double porosity = iq0.porosity().value();
243 const double volume = elem_ctx.dofTotalVolume(0, 0);
244 // TODO: not sure we should use water pressure here
245 const double water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
246 const double water_volume = volume * porosity * water_saturation;
249
251 }
252 OPM_END_PARALLEL_TRY_CATCH("AquiferNumerical::calculateAquiferPressure() failed: ", this->ebos_simulator_.vanguard().grid().comm());
253 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
254 comm.sum(&sum_pressure_watervolume, 1);
255 comm.sum(&sum_watervolume, 1);
256
257 // Ensure all processes have same notion of the aquifer cells' pressure values.
258 comm.sum(cell_pressure.data(), cell_pressure.size());
259
261 }
262
263 template <class ElemCtx>
264 double getWaterFlux(const ElemCtx& elem_ctx, unsigned face_idx) const
265 {
266 const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
267 const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
268 return water_flux;
269 }
270
271 double calculateAquiferFluxRate() const
272 {
273 double aquifer_flux = 0.0;
274
275 if (! this->connects_to_reservoir_) {
276 return aquifer_flux;
277 }
278
279 ElementContext elem_ctx(this->ebos_simulator_);
280 const auto& gridView = this->ebos_simulator_.gridView();
281 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
282 // elem_ctx.updatePrimaryStencil(elem);
283 elem_ctx.updateStencil(elem);
284
285 const std::size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
286 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
287 // we only need the first aquifer cell
288 if (idx != 0) {
289 continue;
290 }
291
292 const std::size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0);
293 // const auto &problem = elem_ctx.problem();
294 const auto& stencil = elem_ctx.stencil(0);
295 // const auto& inQuants = elem_ctx.intensiveQuantities(0, /*timeIdx*/ 0);
296
297 for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
298 const auto& face = stencil.interiorFace(face_idx);
299 // dof index
300 const std::size_t i = face.interiorIndex();
301 const std::size_t j = face.exteriorIndex();
302 // compressed index
303 // const std::size_t I = stencil.globalSpaceIndex(i);
304 const std::size_t J = stencil.globalSpaceIndex(j);
305
306 assert(stencil.globalSpaceIndex(i) == cell_index);
307
308 // we do not consider the flux within aquifer cells
309 // we only need the flux to the connections
310 if (this->cell_to_aquifer_cell_idx_[J] > 0) {
311 continue;
312 }
313 elem_ctx.updateAllIntensiveQuantities();
314 elem_ctx.updateAllExtensiveQuantities();
315
316 const double water_flux = getWaterFlux(elem_ctx,face_idx);
317 const std::size_t up_id = water_flux >= 0.0 ? i : j;
318 const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
319 const double invB = Toolbox::value(intQuantsIn.fluidState().invB(this->phaseIdx_()));
320 const double face_area = face.area();
322 }
323
324 // we only need to handle the first aquifer cell, we can exit loop here
325 break;
326 }
327
328 return aquifer_flux;
329 }
330
331 double flux_rate_; // aquifer influx rate
332 double cumulative_flux_; // cumulative aquifer influx
333 std::vector<double> init_pressure_{};
334 double pressure_; // aquifer pressure
335 bool solution_set_from_restart_ {false};
336 bool connects_to_reservoir_ {false};
337
338 // TODO: maybe unordered_map can also do the work to save memory?
339 std::vector<int> cell_to_aquifer_cell_idx_;
340};
341
342} // namespace Opm
343
344#endif
Definition AquiferInterface.hpp:35
Definition AquiferNumerical.hpp:45
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27