My Project
Loading...
Searching...
No Matches
AquiferConstantFlux.hpp
1/*
2 Copyright (C) 2023 Equinor
3
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#ifndef OPM_AQUIFERCONSTANTFLUX_HPP
21#define OPM_AQUIFERCONSTANTFLUX_HPP
22
23#include <opm/simulators/aquifers/AquiferInterface.hpp>
24
25#include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
26#include <opm/input/eclipse/EclipseState/Aquifer/AquiferFlux.hpp>
27#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
28
29#include <opm/common/ErrorMacros.hpp>
30
31#include <opm/material/common/MathToolbox.hpp>
32#include <opm/material/densead/Evaluation.hpp>
33
34#include <algorithm>
35#include <cassert>
36#include <numeric>
37#include <stdexcept>
38#include <vector>
39
40namespace Opm {
41
42template<typename TypeTag>
44{
45public:
51
52 static constexpr int numEq = BlackoilIndices::numEq;
53 using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
54
55 AquiferConstantFlux(const std::vector<Aquancon::AquancCell>& connections,
59 , connections_ (connections)
60 , aquifer_data_ (aquifer)
61 , connection_flux_ (connections_.size(), Eval{0})
62 {
63 this->total_face_area_ = this->initializeConnections();
64 }
65
66 static AquiferConstantFlux serializationTestObject(const Simulator& ebos_simulator)
67 {
69 result.cumulative_flux_ = 1.0;
70
71 return result;
72 }
73
74 virtual ~AquiferConstantFlux() = default;
75
76 void computeFaceAreaFraction(const std::vector<double>& total_face_area) override
77 {
78 assert (total_face_area.size() >= static_cast<std::vector<double>::size_type>(this->aquiferID()));
79
80 this->area_fraction_ = this->totalFaceArea()
81 / total_face_area[this->aquiferID() - 1];
82 }
83
84 double totalFaceArea() const override
85 {
86 return this->total_face_area_;
87 }
88
89 void updateAquifer(const SingleAquiferFlux& aquifer)
90 {
91 aquifer_data_ = aquifer;
92 }
93
94 void initFromRestart(const data::Aquifers& aquiferSoln) override
95 {
96 auto xaqPos = aquiferSoln.find(this->aquiferID());
97 if (xaqPos == aquiferSoln.end()) {
98 return;
99 }
100
101 this->cumulative_flux_ = this->area_fraction_ * xaqPos->second.volume;
102 }
103
104 void initialSolutionApplied() override
105 {}
106
107 void beginTimeStep() override
108 {}
109
110 void endTimeStep() override
111 {
112 this->flux_rate_ = this->totalFluxRate();
113 this->cumulative_flux_ +=
114 this->flux_rate_ * this->ebos_simulator_.timeStepSize();
115 }
116
117 data::AquiferData aquiferData() const override
118 {
119 data::AquiferData data;
120
121 data.aquiferID = this->aquifer_data_.id;
122
123 // Pressure for constant flux aquifer is 0
124 data.pressure = 0.0;
125 data.fluxRate = this->totalFluxRate();
126
127 data.volume = this->cumulative_flux_;
128
129 // not totally sure whether initPressure matters
130 data.initPressure = 0.0;
131
132 return data;
133 }
134
135 void addToSource(RateVector& rates,
136 const unsigned cellIdx,
137 [[maybe_unused]] const unsigned timeIdx) override
138 {
139 const int idx = this->cellToConnectionIdx_[cellIdx];
140 if (idx < 0) {
141 return;
142 }
143
144 const auto& model = this->ebos_simulator_.model();
145
146 const auto fw = this->aquifer_data_.flux;
147
148 this->connection_flux_[idx] = fw * this->connections_[idx].effective_facearea;
149
150 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
151 += this->connection_flux_[idx] / model.dofTotalVolume(cellIdx);
152 }
153
154 template<class Serializer>
155 void serializeOp(Serializer& serializer)
156 {
157 serializer(cumulative_flux_);
158 }
159
160 bool operator==(const AquiferConstantFlux& rhs) const
161 {
162 return this->cumulative_flux_ == rhs.cumulative_flux_;
163 }
164
165private:
166 const std::vector<Aquancon::AquancCell>& connections_;
167
168 SingleAquiferFlux aquifer_data_;
169 std::vector<Eval> connection_flux_{};
170 std::vector<int> cellToConnectionIdx_{};
171 double flux_rate_{};
172 double cumulative_flux_{};
173 double total_face_area_{0.0};
174 double area_fraction_{1.0};
175
176 double initializeConnections()
177 {
178 auto connected_face_area = 0.0;
179
180 this->cellToConnectionIdx_
181 .resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
182
183 for (std::size_t idx = 0; idx < this->connections_.size(); ++idx) {
184 const auto global_index = this->connections_[idx].global_index;
185 const int cell_index = this->ebos_simulator_.vanguard()
186 .compressedIndexForInterior(global_index);
187
188 if (cell_index < 0) {
189 continue;
190 }
191
192 this->cellToConnectionIdx_[cell_index] = idx;
193
194 connected_face_area += this->connections_[idx].effective_facearea;
195 }
196
197 // TODO: At the moment, we are using the effective_facearea from the
198 // parser. Should we update the facearea here if the grid changed
199 // during the preprocessing?
200
201 return connected_face_area;
202 }
203
204 double computeFaceAreaFraction(const double connected_face_area) const
205 {
206 const auto tot_face_area = this->ebos_simulator_.vanguard()
207 .grid().comm().sum(connected_face_area);
208
209 return (tot_face_area > 0.0)
211 : 0.0;
212 }
213
214 // TODO: this is a function from AquiferAnalytical
215 int compIdx_() const
216 {
217 if (this->co2store_or_h2store_())
218 return FluidSystem::oilCompIdx;
219
220 return FluidSystem::waterCompIdx;
221 }
222
223 double totalFluxRate() const
224 {
225 return std::accumulate(this->connection_flux_.begin(),
226 this->connection_flux_.end(), 0.0,
227 [](const double rate, const auto& q)
228 {
229 return rate + getValue(q);
230 });
231 }
232};
233
234} // namespace Opm
235
236#endif //OPM_AQUIFERCONSTANTFLUX_HPP
Definition AquiferConstantFlux.hpp:44
Definition AquiferInterface.hpp:35
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27