My Project
Loading...
Searching...
No Matches
AquiferCarterTracy.hpp
1/*
2 Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface
3 Copyright 2017 Statoil ASA.
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_AQUIFERCT_HEADER_INCLUDED
22#define OPM_AQUIFERCT_HEADER_INCLUDED
23
24#include <opm/simulators/aquifers/AquiferAnalytical.hpp>
25
26#include <opm/input/eclipse/EclipseState/Aquifer/AquiferCT.hpp>
27
28#include <opm/output/data/Aquifer.hpp>
29
30#include <exception>
31#include <memory>
32#include <stdexcept>
33#include <utility>
34
35namespace Opm
36{
37
38template <typename TypeTag>
40{
41public:
43
46 using typename Base::Eval;
54
55 AquiferCarterTracy(const std::vector<Aquancon::AquancCell>& connections,
56 const Simulator& ebosSimulator,
57 const AquiferCT::AQUCT_data& aquct_data)
58 : Base(aquct_data.aquiferID, connections, ebosSimulator)
59 , aquct_data_(aquct_data)
60 {}
61
62 static AquiferCarterTracy serializationTestObject(const Simulator& ebosSimulator)
63 {
64 AquiferCarterTracy result({}, ebosSimulator, {});
65
66 result.pressure_previous_ = {1.0, 2.0, 3.0};
67 result.pressure_current_ = {4.0, 5.0};
68 result.Qai_ = {{6.0}};
69 result.rhow_ = 7.0;
70 result.W_flux_ = 8.0;
71 result.fluxValue_ = 9.0;
72 result.dimensionless_time_ = 10.0;
73 result.dimensionless_pressure_ = 11.0;
74
75 return result;
76 }
77
78 void endTimeStep() override
79 {
80 for (const auto& q : this->Qai_) {
81 this->W_flux_ += q * this->ebos_simulator_.timeStepSize();
82 }
83 this->fluxValue_ = this->W_flux_.value();
84 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
85 comm.sum(&this->fluxValue_, 1);
86 }
87
88 data::AquiferData aquiferData() const override
89 {
90 data::AquiferData data;
91 data.aquiferID = this->aquiferID();
92 // TODO: not sure how to get this pressure value yet
93 data.pressure = this->pa0_;
94 data.fluxRate = 0.;
95 for (const auto& q : this->Qai_) {
96 data.fluxRate += q.value();
97 }
98 data.volume = this->W_flux_.value();
99 data.initPressure = this->pa0_;
100
101 auto* aquCT = data.typeData.template create<data::AquiferType::CarterTracy>();
102
103 aquCT->dimensionless_time = this->dimensionless_time_;
104 aquCT->dimensionless_pressure = this->dimensionless_pressure_;
105 aquCT->influxConstant = this->aquct_data_.influxConstant();
106
107 if (!this->co2store_or_h2store_()) {
108 aquCT->timeConstant = this->aquct_data_.timeConstant();
109 aquCT->waterDensity = this->aquct_data_.waterDensity();
110 aquCT->waterViscosity = this->aquct_data_.waterViscosity();
111 } else {
112 aquCT->waterDensity = this->rhow_;
113 aquCT->timeConstant = this->Tc_;
114 const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
115 aquCT->waterViscosity = this->Tc_ * this->aquct_data_.permeability / x;
116 }
117
118 return data;
119 }
120
121 template<class Serializer>
122 void serializeOp(Serializer& serializer)
123 {
124 serializer(static_cast<Base&>(*this));
125 serializer(fluxValue_);
126 serializer(dimensionless_time_);
127 serializer(dimensionless_pressure_);
128 }
129
130 bool operator==(const AquiferCarterTracy& rhs) const
131 {
132 return static_cast<const AquiferAnalytical<TypeTag>&>(*this) == rhs &&
133 this->fluxValue_ == rhs.fluxValue_ &&
134 this->dimensionless_time_ == rhs.dimensionless_time_ &&
135 this->dimensionless_pressure_ == rhs.dimensionless_pressure_;
136 }
137
138protected:
139 // Variables constants
140 AquiferCT::AQUCT_data aquct_data_;
141
142 Scalar beta_; // Influx constant
143 // TODO: it is possible it should be a AD variable
144 Scalar fluxValue_{0}; // value of flux
145
146 Scalar dimensionless_time_{0};
147 Scalar dimensionless_pressure_{0};
148
149 void assignRestartData(const data::AquiferData& xaq) override
150 {
151 this->fluxValue_ = xaq.volume;
152 this->rhow_ = this->aquct_data_.waterDensity();
153 }
154
155 std::pair<Scalar, Scalar>
156 getInfluenceTableValues(const Scalar td_plus_dt)
157 {
158 // We use the opm-common numeric linear interpolator
159 this->dimensionless_pressure_ =
160 linearInterpolation(this->aquct_data_.dimensionless_time,
161 this->aquct_data_.dimensionless_pressure,
162 this->dimensionless_time_);
163
164 const auto PItd =
165 linearInterpolation(this->aquct_data_.dimensionless_time,
166 this->aquct_data_.dimensionless_pressure,
167 td_plus_dt);
168
169 const auto PItdprime =
170 linearInterpolationDerivative(this->aquct_data_.dimensionless_time,
171 this->aquct_data_.dimensionless_pressure,
172 td_plus_dt);
173
174 return std::make_pair(PItd, PItdprime);
175 }
176
177 Scalar dpai(const int idx) const
178 {
179 const auto gdz =
180 this->gravity_() * (this->cell_depth_.at(idx) - this->aquiferDepth());
181
182 const auto dp = this->pa0_ + this->rhow_*gdz
183 - this->pressure_previous_.at(idx);
184
185 return dp;
186 }
187
188 // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription
189 std::pair<Scalar, Scalar>
190 calculateEqnConstants(const int idx, const Simulator& simulator)
191 {
192 const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / this->Tc_;
193 this->dimensionless_time_ = simulator.time() / this->Tc_;
194
195 const auto [PItd, PItdprime] = this->getInfluenceTableValues(td_plus_dt);
196
197 const auto denom = this->Tc_ * (PItd - this->dimensionless_time_*PItdprime);
198 const auto a = (this->beta_*dpai(idx) - this->fluxValue_*PItdprime) / denom;
199 const auto b = this->beta_ / denom;
200
201 return std::make_pair(a, b);
202 }
203
204 std::size_t pvtRegionIdx() const
205 {
206 return this->aquct_data_.pvttableID - 1;
207 }
208
209 // This function implements Eq 5.7 of the EclipseTechnicalDescription
210 void calculateInflowRate(int idx, const Simulator& simulator) override
211 {
212 const auto [a, b] = this->calculateEqnConstants(idx, simulator);
213
214 this->Qai_.at(idx) = this->alphai_.at(idx) *
215 (a - b*(this->pressure_current_.at(idx) - this->pressure_previous_.at(idx)));
216 }
217
218 void calculateAquiferConstants() override
219 {
220 this->Tc_ = this->co2store_or_h2store_()
221 ? this->timeConstantCO2Store()
222 : this->aquct_data_.timeConstant();
223
224 this->beta_ = this->aquct_data_.influxConstant();
225 }
226
227 void calculateAquiferCondition() override
228 {
229 if (this->solution_set_from_restart_) {
230 return;
231 }
232
233 if (! this->aquct_data_.initial_pressure.has_value()) {
234 this->aquct_data_.initial_pressure =
235 this->calculateReservoirEquilibrium();
236
237 const auto& tables = this->ebos_simulator_.vanguard()
238 .eclState().getTableManager();
239
240 this->aquct_data_.finishInitialisation(tables);
241 }
242
243 this->pa0_ = this->aquct_data_.initial_pressure.value();
244 if (this->aquct_data_.initial_temperature.has_value()) {
245 this->Ta0_ = this->aquct_data_.initial_temperature.value();
246 }
247
248 this->rhow_ = this->co2store_or_h2store_()
249 ? this->waterDensityCO2Store()
250 : this->aquct_data_.waterDensity();
251 }
252
253 Scalar aquiferDepth() const override
254 {
255 return this->aquct_data_.datum_depth;
256 }
257
258private:
259 Scalar timeConstantCO2Store() const
260 {
261 const auto press = this->aquct_data_.initial_pressure.value();
262 const auto temp = this->reservoirTemperatureCO2Store();
263
264 auto waterViscosity = Scalar { 0 };
265 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
266 const auto rs = Scalar { 0 }; // no dissolved CO2
267 waterViscosity = FluidSystem::oilPvt()
268 .viscosity(pvtRegionIdx(), temp, press, rs);
269 }
270 else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
271 const auto salt = Scalar { 0 };
272 const auto rsw = Scalar { 0 };
273 waterViscosity = FluidSystem::waterPvt()
274 .viscosity(pvtRegionIdx(), temp, press, rsw, salt);
275 }
276 else {
277 OPM_THROW(std::runtime_error, "water or oil phase is needed to run CO2Store.");
278 }
279
280 const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr
281 * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
282
283 return waterViscosity * x / this->aquct_data_.permeability;
284 }
285
286 Scalar waterDensityCO2Store() const
287 {
288 const auto press = this->aquct_data_.initial_pressure.value();
289 const auto temp = this->reservoirTemperatureCO2Store();
290
291 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
292 const auto& pvt = FluidSystem::oilPvt();
293 const auto reg = this->pvtRegionIdx();
294
295 const auto rs = Scalar { 0 }; // no dissolved CO2
296 return pvt.inverseFormationVolumeFactor(reg, temp, press, rs)
297 * pvt.oilReferenceDensity(reg);
298 }
299 else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
300 const auto& pvt = FluidSystem::waterPvt();
301 const auto reg = this->pvtRegionIdx();
302
303 const auto salinity = Scalar { 0 };
304 const auto rsw = Scalar { 0 };
305
306 return pvt.inverseFormationVolumeFactor(reg, temp, press, rsw, salinity)
307 * pvt.waterReferenceDensity(reg);
308 }
309 else {
310 OPM_THROW(std::runtime_error, "water or oil phase is needed to run CO2Store.");
311 }
312 }
313
314 Scalar reservoirTemperatureCO2Store() const
315 {
316 return this->aquct_data_.initial_temperature.has_value()
317 ? this->aquct_data_.initial_temperature.value()
318 : FluidSystem::reservoirTemperature();
319 }
320
321}; // class AquiferCarterTracy
322
323} // namespace Opm
324
325#endif
Definition AquiferAnalytical.hpp:55
Definition AquiferCarterTracy.hpp:40
Definition AquiferInterface.hpp:35
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27