My Project
EclThermalLawManager.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
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 2 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  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_ECL_THERMAL_LAW_MANAGER_HPP
28 #define OPM_ECL_THERMAL_LAW_MANAGER_HPP
29 
30 #if ! HAVE_ECL_INPUT
31 #error "Eclipse input support in opm-common is required to use the ECL thermal law manager!"
32 #endif
33 
36 
39 
40 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
41 #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
42 
43 #include <cassert>
44 #include <stdexcept>
45 #include <vector>
46 
47 namespace Opm {
48 
55 template <class Scalar, class FluidSystem>
57 {
58 public:
60  using SolidEnergyLawParams = typename SolidEnergyLaw::Params;
61  using HeatcrLawParams = typename SolidEnergyLawParams::HeatcrLawParams;
62  using SpecrockLawParams = typename SolidEnergyLawParams::SpecrockLawParams;
63 
65  using ThermalConductionLawParams = typename ThermalConductionLaw::Params;
66 
68  {
69  solidEnergyApproach_ = SolidEnergyLawParams::undefinedApproach;
70  thermalConductivityApproach_ = ThermalConductionLawParams::undefinedApproach;
71  }
72 
73  void initParamsForElements(const EclipseState& eclState, size_t numElems)
74  {
75  const auto& fp = eclState.fieldProps();
76  const auto& tableManager = eclState.getTableManager();
77  bool has_heatcr = fp.has_double("HEATCR");
78  bool has_thconr = fp.has_double("THCONR");
79  bool has_thc = fp.has_double("THCROCK") || fp.has_double("THCOIL") || fp.has_double("THCGAS") || fp.has_double("THCWATER");
80 
81  if (has_heatcr)
82  initHeatcr_(eclState, numElems);
83  else if (tableManager.hasTables("SPECROCK"))
84  initSpecrock_(eclState, numElems);
85  else
86  initNullRockEnergy_();
87 
88  if (has_thconr)
89  initThconr_(eclState, numElems);
90  else if (has_thc)
91  initThc_(eclState, numElems);
92  else
93  initNullCond_();
94  }
95 
96  const SolidEnergyLawParams& solidEnergyLawParams(unsigned elemIdx) const
97  {
98  switch (solidEnergyApproach_) {
99  case SolidEnergyLawParams::heatcrApproach:
100  assert(elemIdx < solidEnergyLawParams_.size());
101  return solidEnergyLawParams_[elemIdx];
102 
103  case SolidEnergyLawParams::specrockApproach:
104  {
105  assert(elemIdx < elemToSatnumIdx_.size());
106  unsigned satnumIdx = elemToSatnumIdx_[elemIdx];
107  assert(satnumIdx < solidEnergyLawParams_.size());
108  return solidEnergyLawParams_[satnumIdx];
109  }
110 
111  case SolidEnergyLawParams::nullApproach:
112  return solidEnergyLawParams_[0];
113 
114  default:
115  throw std::runtime_error("Attempting to retrieve solid energy storage parameters "
116  "without a known approach being defined by the deck.");
117  }
118  }
119 
120  const ThermalConductionLawParams& thermalConductionLawParams(unsigned elemIdx) const
121  {
122  switch (thermalConductivityApproach_) {
123  case ThermalConductionLawParams::thconrApproach:
124  case ThermalConductionLawParams::thcApproach:
125  assert(elemIdx < thermalConductionLawParams_.size());
126  return thermalConductionLawParams_[elemIdx];
127 
128  case ThermalConductionLawParams::nullApproach:
129  return thermalConductionLawParams_[0];
130 
131  default:
132  throw std::runtime_error("Attempting to retrieve thermal conduction parameters without "
133  "a known approach being defined by the deck.");
134  }
135  }
136 
137 private:
141  void initHeatcr_(const EclipseState& eclState,
142  size_t numElems)
143  {
144  solidEnergyApproach_ = SolidEnergyLawParams::heatcrApproach;
145  // actually the value of the reference temperature does not matter for energy
146  // conservation. We set it anyway to faciliate comparisons with ECL
147  HeatcrLawParams::setReferenceTemperature(FluidSystem::surfaceTemperature);
148 
149  const auto& fp = eclState.fieldProps();
150  const std::vector<double>& heatcrData = fp.get_double("HEATCR");
151  const std::vector<double>& heatcrtData = fp.get_double("HEATCRT");
152  solidEnergyLawParams_.resize(numElems);
153  for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
154  auto& elemParam = solidEnergyLawParams_[elemIdx];
155  elemParam.setSolidEnergyApproach(SolidEnergyLawParams::heatcrApproach);
156  auto& heatcrElemParams = elemParam.template getRealParams<SolidEnergyLawParams::heatcrApproach>();
157 
158  heatcrElemParams.setReferenceRockHeatCapacity(heatcrData[elemIdx]);
159  heatcrElemParams.setDRockHeatCapacity_dT(heatcrtData[elemIdx]);
160  heatcrElemParams.finalize();
161  elemParam.finalize();
162  }
163  }
164 
168  void initSpecrock_(const EclipseState& eclState,
169  size_t numElems)
170  {
171  solidEnergyApproach_ = SolidEnergyLawParams::specrockApproach;
172 
173  // initialize the element index -> SATNUM index mapping
174  const auto& fp = eclState.fieldProps();
175  const std::vector<int>& satnumData = fp.get_int("SATNUM");
176  elemToSatnumIdx_.resize(numElems);
177  for (unsigned elemIdx = 0; elemIdx < numElems; ++ elemIdx) {
178  // satnumData contains Fortran-style indices, i.e., they start with 1 instead
179  // of 0!
180  elemToSatnumIdx_[elemIdx] = satnumData[elemIdx] - 1;
181  }
182  // internalize the SPECROCK table
183  unsigned numSatRegions = eclState.runspec().tabdims().getNumSatTables();
184  const auto& tableManager = eclState.getTableManager();
185  solidEnergyLawParams_.resize(numSatRegions);
186  for (unsigned satnumIdx = 0; satnumIdx < numSatRegions; ++satnumIdx) {
187  const auto& specrockTable = tableManager.getSpecrockTables()[satnumIdx];
188 
189  auto& multiplexerParams = solidEnergyLawParams_[satnumIdx];
190 
191  multiplexerParams.setSolidEnergyApproach(SolidEnergyLawParams::specrockApproach);
192 
193  auto& specrockParams = multiplexerParams.template getRealParams<SolidEnergyLawParams::specrockApproach>();
194  const auto& temperatureColumn = specrockTable.getColumn("TEMPERATURE");
195  const auto& cvRockColumn = specrockTable.getColumn("CV_ROCK");
196  specrockParams.setHeatCapacities(temperatureColumn, cvRockColumn);
197  specrockParams.finalize();
198 
199  multiplexerParams.finalize();
200  }
201  }
202 
206  void initNullRockEnergy_()
207  {
208  solidEnergyApproach_ = SolidEnergyLawParams::nullApproach;
209 
210  solidEnergyLawParams_.resize(1);
211  solidEnergyLawParams_[0].finalize();
212  }
213 
217  void initThconr_(const EclipseState& eclState,
218  size_t numElems)
219  {
220  thermalConductivityApproach_ = ThermalConductionLawParams::thconrApproach;
221 
222  const auto& fp = eclState.fieldProps();
223  std::vector<double> thconrData;
224  std::vector<double> thconsfData;
225  if (fp.has_double("THCONR"))
226  thconrData = fp.get_double("THCONR");
227 
228  if (fp.has_double("THCONSF"))
229  thconsfData = fp.get_double("THCONSF");
230 
231  thermalConductionLawParams_.resize(numElems);
232  for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
233  auto& elemParams = thermalConductionLawParams_[elemIdx];
234  elemParams.setThermalConductionApproach(ThermalConductionLawParams::thconrApproach);
235  auto& thconrElemParams = elemParams.template getRealParams<ThermalConductionLawParams::thconrApproach>();
236 
237  double thconr = thconrData.empty() ? 0.0 : thconrData[elemIdx];
238  double thconsf = thconsfData.empty() ? 0.0 : thconsfData[elemIdx];
239  thconrElemParams.setReferenceTotalThermalConductivity(thconr);
240  thconrElemParams.setDTotalThermalConductivity_dSg(thconsf);
241 
242  thconrElemParams.finalize();
243  elemParams.finalize();
244  }
245  }
246 
250  void initThc_(const EclipseState& eclState,
251  size_t numElems)
252  {
253  thermalConductivityApproach_ = ThermalConductionLawParams::thcApproach;
254 
255  const auto& fp = eclState.fieldProps();
256  std::vector<double> thcrockData;
257  std::vector<double> thcoilData;
258  std::vector<double> thcgasData;
259  std::vector<double> thcwaterData = fp.get_double("THCWATER");
260 
261  if (fp.has_double("THCROCK"))
262  thcrockData = fp.get_double("THCROCK");
263 
264  if (fp.has_double("THCOIL"))
265  thcoilData = fp.get_double("THCOIL");
266 
267  if (fp.has_double("THCGAS"))
268  thcgasData = fp.get_double("THCGAS");
269 
270  if (fp.has_double("THCWATER"))
271  thcwaterData = fp.get_double("THCWATER");
272 
273  const std::vector<double>& poroData = fp.get_double("PORO");
274 
275  thermalConductionLawParams_.resize(numElems);
276  for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
277  auto& elemParams = thermalConductionLawParams_[elemIdx];
278  elemParams.setThermalConductionApproach(ThermalConductionLawParams::thcApproach);
279  auto& thcElemParams = elemParams.template getRealParams<ThermalConductionLawParams::thcApproach>();
280 
281  thcElemParams.setPorosity(poroData[elemIdx]);
282  double thcrock = thcrockData.empty() ? 0.0 : thcrockData[elemIdx];
283  double thcoil = thcoilData.empty() ? 0.0 : thcoilData[elemIdx];
284  double thcgas = thcgasData.empty() ? 0.0 : thcgasData[elemIdx];
285  double thcwater = thcwaterData.empty() ? 0.0 : thcwaterData[elemIdx];
286  thcElemParams.setThcrock(thcrock);
287  thcElemParams.setThcoil(thcoil);
288  thcElemParams.setThcgas(thcgas);
289  thcElemParams.setThcwater(thcwater);
290 
291  thcElemParams.finalize();
292  elemParams.finalize();
293  }
294  }
295 
299  void initNullCond_()
300  {
301  thermalConductivityApproach_ = ThermalConductionLawParams::nullApproach;
302 
303  thermalConductionLawParams_.resize(1);
304  thermalConductionLawParams_[0].finalize();
305  }
306 
307 private:
308  typename ThermalConductionLawParams::ThermalConductionApproach thermalConductivityApproach_;
309  typename SolidEnergyLawParams::SolidEnergyApproach solidEnergyApproach_;
310 
311  std::vector<unsigned> elemToSatnumIdx_;
312 
313  std::vector<SolidEnergyLawParams> solidEnergyLawParams_;
314  std::vector<ThermalConductionLawParams> thermalConductionLawParams_;
315 };
316 } // namespace Opm
317 
318 #endif
The default implementation of a parameter object for the ECL thermal law.
Provides the energy storage relation of rock.
The default implementation of a parameter object for the ECL thermal law.
Implements the total thermal conductivity and rock enthalpy relations used by ECL.
Provides the energy storage relation of rock.
Definition: EclSolidEnergyLawMultiplexer.hpp:50
Implements the total thermal conductivity and rock enthalpy relations used by ECL.
Definition: EclThermalConductionLawMultiplexer.hpp:50
Provides an simple way to create and manage the thermal law objects for a complete ECL deck.
Definition: EclThermalLawManager.hpp:57