My Project
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 <numeric>
36#include <stdexcept>
37
38namespace Opm {
39
40template<typename TypeTag>
42{
43public:
44 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
45 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
46 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
47 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
48 using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
49
50 static constexpr int numEq = BlackoilIndices::numEq;
51 using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
52
53 AquiferConstantFlux(const std::vector<Aquancon::AquancCell>& connections,
54 const Simulator& ebos_simulator,
55 const SingleAquiferFlux& aquifer)
56 : AquiferInterface<TypeTag>(aquifer.id, ebos_simulator)
57 , connections_ (connections)
58 , aquifer_data_ (aquifer)
59 , connection_flux_ (connections_.size(), Eval{0})
60 {
61 this->initializeConnections();
62 }
63
64 static AquiferConstantFlux serializationTestObject(const Simulator& ebos_simulator)
65 {
66 AquiferConstantFlux<TypeTag> result({}, ebos_simulator, {});
67 result.cumulative_flux_ = 1.0;
68
69 return result;
70 }
71
72 virtual ~AquiferConstantFlux() = default;
73
74 void updateAquifer(const SingleAquiferFlux& aquifer)
75 {
76 aquifer_data_ = aquifer;
77 }
78
79 void initFromRestart(const data::Aquifers& /* aquiferSoln */) override {
80 }
81
82 void initialSolutionApplied() override {
83 }
84
85
86 void beginTimeStep() override
87 {}
88
89 void endTimeStep() override
90 {
91 this->flux_rate_ = this->totalFluxRate();
92 this->cumulative_flux_ +=
93 this->flux_rate_ * this->ebos_simulator_.timeStepSize();
94 }
95
96 data::AquiferData aquiferData() const override
97 {
98 data::AquiferData data;
99
100 data.aquiferID = this->aquifer_data_.id;
101
102 // Pressure for constant flux aquifer is 0
103 data.pressure = 0.0;
104 data.fluxRate = this->totalFluxRate();
105
106 data.volume = this->cumulative_flux_;
107
108 // not totally sure whether initPressure matters
109 data.initPressure = 0.0;
110
111 return data;
112 }
113
114 void addToSource(RateVector& rates,
115 const unsigned cellIdx,
116 [[maybe_unused]] const unsigned timeIdx) override
117 {
118 const int idx = this->cellToConnectionIdx_[cellIdx];
119 if (idx < 0) {
120 return;
121 }
122
123 const auto& model = this->ebos_simulator_.model();
124
125 const auto fw = this->aquifer_data_.flux;
126
127 this->connection_flux_[idx] = fw * this->connections_[idx].effective_facearea;
128
129 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
130 += this->connection_flux_[idx] / model.dofTotalVolume(cellIdx);
131 }
132
133 template<class Serializer>
134 void serializeOp(Serializer& serializer)
135 {
136 serializer(cumulative_flux_);
137 }
138
139 bool operator==(const AquiferConstantFlux& rhs) const
140 {
141 return this->cumulative_flux_ == rhs.cumulative_flux_;
142 }
143
144private:
145 const std::vector<Aquancon::AquancCell>& connections_;
146
147 SingleAquiferFlux aquifer_data_;
148 std::vector<Eval> connection_flux_{};
149 std::vector<int> cellToConnectionIdx_{};
150 double flux_rate_{};
151 double cumulative_flux_{};
152
153 void initializeConnections() {
154 this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
155 for (std::size_t idx = 0; idx < this->connections_.size(); ++idx) {
156 const auto global_index = this->connections_[idx].global_index;
157 const int cell_index = this->ebos_simulator_.vanguard()
158 .compressedIndexForInterior(global_index);
159
160 if (cell_index < 0) {
161 continue;
162 }
163
164 this->cellToConnectionIdx_[cell_index] = idx;
165 }
166
167 // TODO: At the moment, we are using the effective_facearea from the
168 // parser. Should we update the facearea here if the grid changed
169 // during the preprocessing?
170 }
171
172 // TODO: this is a function from AquiferAnalytical
173 int compIdx_() const
174 {
175 if (this->co2store_())
176 return FluidSystem::oilCompIdx;
177
178 return FluidSystem::waterCompIdx;
179 }
180
181 double totalFluxRate() const
182 {
183 return std::accumulate(this->connection_flux_.begin(),
184 this->connection_flux_.end(), 0.0,
185 [](const double rate, const auto& q)
186 {
187 return rate + getValue(q);
188 });
189 }
190};
191
192} // namespace Opm
193
194#endif //OPM_AQUIFERCONSTANTFLUX_HPP
Definition: AquiferConstantFlux.hpp:42
Definition: AquiferInterface.hpp:35
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27