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_()) {
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 inline 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 inline void calculateAquiferConstants() override
219 {
220 if(this->co2store_()) {
221 const auto press = this->aquct_data_.initial_pressure.value();
222 Scalar temp = FluidSystem::reservoirTemperature();
223 if (this->aquct_data_.initial_temperature.has_value())
224 temp = this->aquct_data_.initial_temperature.value();
225
227 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
228 Scalar rs = 0.0; // no dissolved CO2
229 waterViscosity = FluidSystem::oilPvt().viscosity(pvtRegionIdx(), temp, press, rs);
230 } else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
231 Scalar salt = 0.;
232 Scalar rsw = 0.;
233 waterViscosity = FluidSystem::waterPvt().viscosity(pvtRegionIdx(), temp, press, rsw, salt);
234 } else {
235 OPM_THROW(std::runtime_error, "water or oil phase is needed to run CO2Store.");
236 }
237 const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
238 this->Tc_ = waterViscosity * x / this->aquct_data_.permeability;
239 } else {
240 this->Tc_ = this->aquct_data_.timeConstant();
241 }
242 this->beta_ = this->aquct_data_.influxConstant();
243 }
244
245 inline void calculateAquiferCondition() override
246 {
247 if (this->solution_set_from_restart_) {
248 return;
249 }
250
251 if (! this->aquct_data_.initial_pressure.has_value()) {
252 this->aquct_data_.initial_pressure =
253 this->calculateReservoirEquilibrium();
254
255 const auto& tables = this->ebos_simulator_.vanguard()
256 .eclState().getTableManager();
257
258 this->aquct_data_.finishInitialisation(tables);
259 }
260
261 this->pa0_ = this->aquct_data_.initial_pressure.value();
262 if (this->aquct_data_.initial_temperature.has_value())
263 this->Ta0_ = this->aquct_data_.initial_temperature.value();
264
265 if(this->co2store_()) {
266 const auto press = this->aquct_data_.initial_pressure.value();
267
268 Scalar temp = FluidSystem::reservoirTemperature();
269 if (this->aquct_data_.initial_temperature.has_value())
270 temp = this->aquct_data_.initial_temperature.value();
271
272 Scalar waterDensity = 0.;
273 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
274 Scalar rs = 0.0; // no dissolved CO2
275 waterDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx(), temp, press, rs)
276 * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIdx());
277 } else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
278 Scalar salinity = 0.;
279 Scalar rsw = 0.;
280 waterDensity = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx(), temp, press, rsw, salinity)
281 * FluidSystem::waterPvt().waterReferenceDensity(pvtRegionIdx());
282 } else {
283 OPM_THROW(std::runtime_error, "water or oil phase is needed to run CO2Store.");
284 }
285 this->rhow_ = waterDensity;
286 } else {
287 this->rhow_ = this->aquct_data_.waterDensity();
288 }
289 }
290
291 virtual Scalar aquiferDepth() const override
292 {
293 return this->aquct_data_.datum_depth;
294 }
295}; // class AquiferCarterTracy
296
297} // namespace Opm
298
299#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