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
28#include <opm/material/common/MathToolbox.hpp>
29#include <opm/material/densead/Evaluation.hpp>
30
31#include <cassert>
32#include <numeric>
33#include <vector>
34
35namespace Opm {
36
37template<typename TypeTag>
39{
40public:
41 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
42 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
43 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
44 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
45 using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
46
47 static constexpr int numEq = BlackoilIndices::numEq;
48 using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
49
50 AquiferConstantFlux(const std::vector<Aquancon::AquancCell>& connections,
51 const Simulator& simulator,
52 const SingleAquiferFlux& aquifer)
53 : AquiferInterface<TypeTag>(aquifer.id, simulator)
54 , connections_ (connections)
55 , aquifer_data_ (aquifer)
56 , connection_flux_ (connections_.size(), Eval{0})
57 {
58 this->total_face_area_ = this->initializeConnections();
59 }
60
61 static AquiferConstantFlux serializationTestObject(const Simulator& simulator)
62 {
63 AquiferConstantFlux<TypeTag> result({}, simulator, {});
64 result.cumulative_flux_ = 1.0;
65
66 return result;
67 }
68
69 virtual ~AquiferConstantFlux() = default;
70
71 void computeFaceAreaFraction(const std::vector<double>& total_face_area) override
72 {
73 assert (total_face_area.size() >= static_cast<std::vector<double>::size_type>(this->aquiferID()));
74
75 this->area_fraction_ = this->totalFaceArea()
76 / total_face_area[this->aquiferID() - 1];
77 }
78
79 double totalFaceArea() const override
80 {
81 return this->total_face_area_;
82 }
83
84 void updateAquifer(const SingleAquiferFlux& aquifer)
85 {
86 aquifer_data_ = aquifer;
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
96 this->cumulative_flux_ = this->area_fraction_ * xaqPos->second.volume;
97 }
98
99 void initialSolutionApplied() override
100 {}
101
102 void beginTimeStep() override
103 {}
104
105 void endTimeStep() override
106 {
107 this->flux_rate_ = this->totalFluxRate();
108 this->cumulative_flux_ +=
109 this->flux_rate_ * this->simulator_.timeStepSize();
110 }
111
112 data::AquiferData aquiferData() const override
113 {
114 data::AquiferData data;
115
116 data.aquiferID = this->aquifer_data_.id;
117
118 // Pressure for constant flux aquifer is 0
119 data.pressure = 0.0;
120 data.fluxRate = this->totalFluxRate();
121
122 data.volume = this->cumulative_flux_;
123
124 // not totally sure whether initPressure matters
125 data.initPressure = 0.0;
126
127 return data;
128 }
129
130 void addToSource(RateVector& rates,
131 const unsigned cellIdx,
132 [[maybe_unused]] const unsigned timeIdx) override
133 {
134 const int idx = this->cellToConnectionIdx_[cellIdx];
135 if (idx < 0) {
136 return;
137 }
138
139 const auto& model = this->simulator_.model();
140
141 const auto fw = this->aquifer_data_.flux;
142
143 this->connection_flux_[idx] = fw * this->connections_[idx].effective_facearea;
144
145 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
146 += this->connection_flux_[idx] / model.dofTotalVolume(cellIdx);
147 }
148
149 template<class Serializer>
150 void serializeOp(Serializer& serializer)
151 {
152 serializer(cumulative_flux_);
153 }
154
155 bool operator==(const AquiferConstantFlux& rhs) const
156 {
157 return this->cumulative_flux_ == rhs.cumulative_flux_;
158 }
159
160private:
161 const std::vector<Aquancon::AquancCell>& connections_;
162
163 SingleAquiferFlux aquifer_data_;
164 std::vector<Eval> connection_flux_{};
165 std::vector<int> cellToConnectionIdx_{};
166 double flux_rate_{};
167 double cumulative_flux_{};
168 double total_face_area_{0.0};
169 double area_fraction_{1.0};
170
171 double initializeConnections()
172 {
173 auto connected_face_area = 0.0;
174
175 this->cellToConnectionIdx_
176 .resize(this->simulator_.gridView().size(/*codim=*/0), -1);
177
178 for (std::size_t idx = 0; idx < this->connections_.size(); ++idx) {
179 const auto global_index = this->connections_[idx].global_index;
180 const int cell_index = this->simulator_.vanguard()
181 .compressedIndexForInterior(global_index);
182
183 if (cell_index < 0) {
184 continue;
185 }
186
187 this->cellToConnectionIdx_[cell_index] = idx;
188
189 connected_face_area += this->connections_[idx].effective_facearea;
190 }
191
192 // TODO: At the moment, we are using the effective_facearea from the
193 // parser. Should we update the facearea here if the grid changed
194 // during the preprocessing?
195
196 return connected_face_area;
197 }
198
199 double computeFaceAreaFraction(const double connected_face_area) const
200 {
201 const auto tot_face_area = this->simulator_.vanguard()
202 .grid().comm().sum(connected_face_area);
203
204 return (tot_face_area > 0.0)
205 ? connected_face_area / tot_face_area
206 : 0.0;
207 }
208
209 // TODO: this is a function from AquiferAnalytical
210 int compIdx_() const
211 {
212 if (this->co2store_or_h2store_())
213 return FluidSystem::oilCompIdx;
214
215 return FluidSystem::waterCompIdx;
216 }
217
218 double totalFluxRate() const
219 {
220 return std::accumulate(this->connection_flux_.begin(),
221 this->connection_flux_.end(), 0.0,
222 [](const double rate, const auto& q)
223 {
224 return rate + getValue(q);
225 });
226 }
227};
228
229} // namespace Opm
230
231#endif //OPM_AQUIFERCONSTANTFLUX_HPP
Definition AquiferConstantFlux.hpp:39
Definition AquiferInterface.hpp:35
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27