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 
33 namespace Opm
34 {
35 
36 template <typename TypeTag>
37 class AquiferCarterTracy : public AquiferAnalytical<TypeTag>
38 {
39 public:
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 
103 protected:
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