My Project
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/output/data/Aquifer.hpp>
27
28#include <exception>
29#include <memory>
30#include <stdexcept>
31#include <utility>
32
33namespace Opm
34{
35
36template <typename TypeTag>
38{
39public:
41
42 using typename Base::BlackoilIndices;
43 using typename Base::ElementContext;
44 using typename Base::Eval;
45 using typename Base::FluidState;
46 using typename Base::FluidSystem;
47 using typename Base::IntensiveQuantities;
48 using typename Base::RateVector;
49 using typename Base::Scalar;
50 using typename Base::Simulator;
51 using typename Base::ElementMapper;
52
53 AquiferCarterTracy(const std::vector<Aquancon::AquancCell>& connections,
54 const Simulator& ebosSimulator,
55 const AquiferCT::AQUCT_data& aquct_data)
56 : Base(aquct_data.aquiferID, connections, ebosSimulator)
57 , aquct_data_(aquct_data)
58 {}
59
60 void endTimeStep() override
61 {
62 for (const auto& q : this->Qai_) {
63 this->W_flux_ += q * this->ebos_simulator_.timeStepSize();
64 }
65 this->fluxValue_ = this->W_flux_.value();
66 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
67 comm.sum(&this->fluxValue_, 1);
68 }
69
70 data::AquiferData aquiferData() const override
71 {
72 data::AquiferData data;
73 data.aquiferID = this->aquiferID();
74 // TODO: not sure how to get this pressure value yet
75 data.pressure = this->pa0_;
76 data.fluxRate = 0.;
77 for (const auto& q : this->Qai_) {
78 data.fluxRate += q.value();
79 }
80 data.volume = this->W_flux_.value();
81 data.initPressure = this->pa0_;
82
83 auto* aquCT = data.typeData.template create<data::AquiferType::CarterTracy>();
84
85 aquCT->dimensionless_time = this->dimensionless_time_;
86 aquCT->dimensionless_pressure = this->dimensionless_pressure_;
87 aquCT->influxConstant = this->aquct_data_.influxConstant();
88
89 if (!this->co2store_()) {
90 aquCT->timeConstant = this->aquct_data_.timeConstant();
91 aquCT->waterDensity = this->aquct_data_.waterDensity();
92 aquCT->waterViscosity = this->aquct_data_.waterViscosity();
93 } else {
94 aquCT->waterDensity = this->rhow_;
95 aquCT->timeConstant = this->Tc_;
96 const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
97 aquCT->waterViscosity = this->Tc_ * this->aquct_data_.permeability / x;
98 }
99
100 return data;
101 }
102
103protected:
104 // Variables constants
105 AquiferCT::AQUCT_data aquct_data_;
106
107 Scalar beta_; // Influx constant
108 // TODO: it is possible it should be a AD variable
109 Scalar fluxValue_{0}; // value of flux
110
111 Scalar dimensionless_time_{0};
112 Scalar dimensionless_pressure_{0};
113
114 void assignRestartData(const data::AquiferData& xaq) override
115 {
116 this->fluxValue_ = xaq.volume;
117 this->rhow_ = this->aquct_data_.waterDensity();
118 }
119
120 std::pair<Scalar, Scalar>
121 getInfluenceTableValues(const Scalar td_plus_dt)
122 {
123 // We use the opm-common numeric linear interpolator
124 this->dimensionless_pressure_ =
125 linearInterpolation(this->aquct_data_.dimensionless_time,
126 this->aquct_data_.dimensionless_pressure,
127 this->dimensionless_time_);
128
129 const auto PItd =
130 linearInterpolation(this->aquct_data_.dimensionless_time,
131 this->aquct_data_.dimensionless_pressure,
132 td_plus_dt);
133
134 const auto PItdprime =
135 linearInterpolationDerivative(this->aquct_data_.dimensionless_time,
136 this->aquct_data_.dimensionless_pressure,
137 td_plus_dt);
138
139 return std::make_pair(PItd, PItdprime);
140 }
141
142 Scalar dpai(const int idx) const
143 {
144 const auto gdz =
145 this->gravity_() * (this->cell_depth_.at(idx) - this->aquiferDepth());
146
147 const auto dp = this->pa0_ + this->rhow_*gdz
148 - this->pressure_previous_.at(idx);
149
150 return dp;
151 }
152
153 // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription
154 std::pair<Scalar, Scalar>
155 calculateEqnConstants(const int idx, const Simulator& simulator)
156 {
157 const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / this->Tc_;
158 this->dimensionless_time_ = simulator.time() / this->Tc_;
159
160 const auto [PItd, PItdprime] = this->getInfluenceTableValues(td_plus_dt);
161
162 const auto denom = this->Tc_ * (PItd - this->dimensionless_time_*PItdprime);
163 const auto a = (this->beta_*dpai(idx) - this->fluxValue_*PItdprime) / denom;
164 const auto b = this->beta_ / denom;
165
166 return std::make_pair(a, b);
167 }
168
169 std::size_t pvtRegionIdx() const
170 {
171 return this->aquct_data_.pvttableID - 1;
172 }
173
174 // This function implements Eq 5.7 of the EclipseTechnicalDescription
175 inline void calculateInflowRate(int idx, const Simulator& simulator) override
176 {
177 const auto [a, b] = this->calculateEqnConstants(idx, simulator);
178
179 this->Qai_.at(idx) = this->alphai_.at(idx) *
180 (a - b*(this->pressure_current_.at(idx) - this->pressure_previous_.at(idx)));
181 }
182
183 inline void calculateAquiferConstants() override
184 {
185 if(this->co2store_()) {
186 const auto press = this->aquct_data_.initial_pressure.value();
187 Scalar temp = FluidSystem::reservoirTemperature();
188 if (this->aquct_data_.initial_temperature.has_value())
189 temp = this->aquct_data_.initial_temperature.value();
190
191 Scalar rs = 0.0; // no dissolved CO2
192 Scalar waterViscosity = FluidSystem::oilPvt().viscosity(pvtRegionIdx(), temp, press, rs);
193 const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
194 this->Tc_ = waterViscosity * x / this->aquct_data_.permeability;
195 } else {
196 this->Tc_ = this->aquct_data_.timeConstant();
197 }
198 this->beta_ = this->aquct_data_.influxConstant();
199 }
200
201 inline void calculateAquiferCondition() override
202 {
203 if (this->solution_set_from_restart_) {
204 return;
205 }
206
207 if (! this->aquct_data_.initial_pressure.has_value()) {
208 this->aquct_data_.initial_pressure =
209 this->calculateReservoirEquilibrium();
210
211 const auto& tables = this->ebos_simulator_.vanguard()
212 .eclState().getTableManager();
213
214 this->aquct_data_.finishInitialisation(tables);
215 }
216
217 this->pa0_ = this->aquct_data_.initial_pressure.value();
218 if (this->aquct_data_.initial_temperature.has_value())
219 this->Ta0_ = this->aquct_data_.initial_temperature.value();
220
221 if(this->co2store_()) {
222 const auto press = this->aquct_data_.initial_pressure.value();
223
224 Scalar temp = FluidSystem::reservoirTemperature();
225 if (this->aquct_data_.initial_temperature.has_value())
226 temp = this->aquct_data_.initial_temperature.value();
227
228 Scalar rs = 0.0; // no dissolved CO2
229 Scalar waterDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx(), temp, press, rs)
230 * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIdx());
231 this->rhow_ = waterDensity;
232 } else {
233 this->rhow_ = this->aquct_data_.waterDensity();
234 }
235 }
236
237 virtual Scalar aquiferDepth() const override
238 {
239 return this->aquct_data_.datum_depth;
240 }
241}; // class AquiferCarterTracy
242
243} // namespace Opm
244
245#endif
Definition: AquiferAnalytical.hpp:51
Definition: AquiferCarterTracy.hpp:38
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27