My Project
RateConverter.hpp
Go to the documentation of this file.
1 /*
2  Copyright 2014, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2014, 2015 Statoil ASA.
4  Copyright 2017, IRIS
5 
6  This file is part of the Open Porous Media Project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 #ifndef OPM_RATECONVERTER_HPP_HEADER_INCLUDED
23 #define OPM_RATECONVERTER_HPP_HEADER_INCLUDED
24 
25 #include <opm/core/props/BlackoilPhases.hpp>
26 #include <opm/grid/utility/RegionMapping.hpp>
27 #include <opm/simulators/wells/RegionAttributeHelpers.hpp>
28 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
29 #include <dune/grid/common/gridenums.hh>
30 #include <dune/grid/common/rangegenerators.hh>
31 #include <algorithm>
32 #include <cmath>
33 #include <memory>
34 #include <stdexcept>
35 #include <type_traits>
36 #include <unordered_map>
37 #include <utility>
38 #include <vector>
49 namespace Opm {
50  namespace RateConverter {
51 
67  template <class FluidSystem, class Region>
69  public:
78  const Region& region)
79  : phaseUsage_(phaseUsage)
80  , rmap_ (region)
81  , attr_ (rmap_, Attributes())
82  {
83  }
84 
85 
94  template <typename ElementContext, class EbosSimulator>
95  void defineState(const EbosSimulator& simulator)
96  {
97 
98  // create map from cell to region
99  // and set all attributes to zero
100  for (const auto& reg : rmap_.activeRegions()) {
101  auto& ra = attr_.attributes(reg);
102  ra.pressure = 0.0;
103  ra.temperature = 0.0;
104  ra.rs = 0.0;
105  ra.rv = 0.0;
106  ra.pv = 0.0;
107  ra.saltConcentration = 0.0;
108 
109  }
110 
111  // quantities for pore volume average
112  std::unordered_map<RegionId, Attributes> attributes_pv;
113 
114  // quantities for hydrocarbon volume average
115  std::unordered_map<RegionId, Attributes> attributes_hpv;
116 
117  for (const auto& reg : rmap_.activeRegions()) {
118  attributes_pv.insert({reg, Attributes()});
119  attributes_hpv.insert({reg, Attributes()});
120  }
121 
122  ElementContext elemCtx( simulator );
123  const auto& gridView = simulator.gridView();
124  const auto& comm = gridView.comm();
125 
126  OPM_BEGIN_PARALLEL_TRY_CATCH();
127  for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
128  elemCtx.updatePrimaryStencil(elem);
129  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
130  const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
131  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
132  const auto& fs = intQuants.fluidState();
133  // use pore volume weighted averages.
134  const double pv_cell =
135  simulator.model().dofTotalVolume(cellIdx)
136  * intQuants.porosity().value();
137 
138  // only count oil and gas filled parts of the domain
139  double hydrocarbon = 1.0;
140  const auto& pu = phaseUsage_;
142  hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
143  }
144 
145  const int reg = rmap_.region(cellIdx);
146  assert(reg >= 0);
147 
148  // sum p, rs, rv, and T.
149  const double hydrocarbonPV = pv_cell*hydrocarbon;
150  if (hydrocarbonPV > 0.) {
151  auto& attr = attributes_hpv[reg];
152  attr.pv += hydrocarbonPV;
154  attr.rs += fs.Rs().value() * hydrocarbonPV;
155  attr.rv += fs.Rv().value() * hydrocarbonPV;
156  }
158  attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
159  attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
160  } else {
162  attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
163  attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
164  }
165  attr.saltConcentration += fs.saltConcentration().value() * hydrocarbonPV;
166  }
167 
168  if (pv_cell > 0.) {
169  auto& attr = attributes_pv[reg];
170  attr.pv += pv_cell;
172  attr.rs += fs.Rs().value() * pv_cell;
173  attr.rv += fs.Rv().value() * pv_cell;
174  }
176  attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
177  attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * pv_cell;
179  attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
180  attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * pv_cell;
181  } else {
183  attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
184  attr.temperature += fs.temperature(FluidSystem::waterPhaseIdx).value() * pv_cell;
185  }
186  attr.saltConcentration += fs.saltConcentration().value() * pv_cell;
187  }
188  }
189 
190  OPM_END_PARALLEL_TRY_CATCH("SurfaceToReservoirVoidage::defineState() failed: ", simulator.vanguard().grid().comm());
191 
192  for (const auto& reg : rmap_.activeRegions()) {
193  auto& ra = attr_.attributes(reg);
194  const double hpv_sum = comm.sum(attributes_hpv[reg].pv);
195  // TODO: should we have some epsilon here instead of zero?
196  if (hpv_sum > 0.) {
197  const auto& attri_hpv = attributes_hpv[reg];
198  const double p_hpv_sum = comm.sum(attri_hpv.pressure);
199  const double T_hpv_sum = comm.sum(attri_hpv.temperature);
200  const double rs_hpv_sum = comm.sum(attri_hpv.rs);
201  const double rv_hpv_sum = comm.sum(attri_hpv.rv);
202  const double sc_hpv_sum = comm.sum(attri_hpv.saltConcentration);
203 
204  ra.pressure = p_hpv_sum / hpv_sum;
205  ra.temperature = T_hpv_sum / hpv_sum;
206  ra.rs = rs_hpv_sum / hpv_sum;
207  ra.rv = rv_hpv_sum / hpv_sum;
208  ra.pv = hpv_sum;
209  ra.saltConcentration = sc_hpv_sum / hpv_sum;
210  } else {
211  // using the pore volume to do the averaging
212  const auto& attri_pv = attributes_pv[reg];
213  const double pv_sum = comm.sum(attri_pv.pv);
214  assert(pv_sum > 0.);
215  const double p_pv_sum = comm.sum(attri_pv.pressure);
216  const double T_pv_sum = comm.sum(attri_pv.temperature);
217  const double rs_pv_sum = comm.sum(attri_pv.rs);
218  const double rv_pv_sum = comm.sum(attri_pv.rv);
219  const double sc_pv_sum = comm.sum(attri_pv.saltConcentration);
220 
221  ra.pressure = p_pv_sum / pv_sum;
222  ra.temperature = T_pv_sum / pv_sum;
223  ra.rs = rs_pv_sum / pv_sum;
224  ra.rv = rv_pv_sum / pv_sum;
225  ra.pv = pv_sum;
226  ra.saltConcentration = sc_pv_sum / pv_sum;
227  }
228  }
229  }
230 
236  typedef typename RegionMapping<Region>::RegionId RegionId;
237 
265  template <class Coeff>
266  void
267  calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
268  {
269  const auto& pu = phaseUsage_;
270  const auto& ra = attr_.attributes(r);
271  const double p = ra.pressure;
272  const double T = ra.temperature;
273  const double saltConcentration = ra.saltConcentration;
274 
275  const int iw = RegionAttributeHelpers::PhasePos::water(pu);
276  const int io = RegionAttributeHelpers::PhasePos::oil (pu);
277  const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
278 
279  std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
280 
282  // q[w]_r = q[w]_s / bw
283 
284  const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
285 
286  coeff[iw] = 1.0 / bw;
287  }
288 
289  // Actual Rs and Rv:
290  double Rs = ra.rs;
291  double Rv = ra.rv;
292 
293  // Determinant of 'R' matrix
294  const double detR = 1.0 - (Rs * Rv);
295 
297  // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
298 
299  const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
300  const double den = bo * detR;
301 
302  coeff[io] += 1.0 / den;
303 
305  coeff[ig] -= ra.rv / den;
306  }
307  }
308 
310  // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
311 
312  const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv, 0.0 /*=Rvw*/);
313  const double den = bg * detR;
314 
315  coeff[ig] += 1.0 / den;
316 
318  coeff[io] -= ra.rs / den;
319  }
320  }
321  }
322 
323  template <class Coeff>
324  void
325  calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
326  {
327  const auto& pu = phaseUsage_;
328  const auto& ra = attr_.attributes(r);
329  const double p = ra.pressure;
330  const double T = ra.temperature;
331  const double saltConcentration = ra.saltConcentration;
332 
333  const int iw = RegionAttributeHelpers::PhasePos::water(pu);
334  const int io = RegionAttributeHelpers::PhasePos::oil (pu);
335  const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
336 
337  std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
338 
340  // q[w]_r = q[w]_s / bw
341 
342  const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
343 
344  coeff[iw] = 1.0 / bw;
345  }
346 
348  const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0);
349  coeff[io] += 1.0 / bo;
350  }
351 
353  const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, 0.0);
354  coeff[ig] += 1.0 / bg;
355  }
356  }
357 
358 
376  template <class Rates >
377  void
378  calcReservoirVoidageRates(const RegionId r, const int pvtRegionIdx, const Rates& surface_rates,
379  Rates& voidage_rates) const
380  {
381  assert(voidage_rates.size() == surface_rates.size());
382 
383  std::fill(voidage_rates.begin(), voidage_rates.end(), 0.0);
384 
385  const auto& pu = phaseUsage_;
386  const auto& ra = attr_.attributes(r);
387  const double p = ra.pressure;
388  const double T = ra.temperature;
389  const double saltConcentration = ra.saltConcentration;
390 
391  const int iw = RegionAttributeHelpers::PhasePos::water(pu);
392  const int io = RegionAttributeHelpers::PhasePos::oil (pu);
393  const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
394 
396  // q[w]_r = q[w]_s / bw
397 
398  const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
399 
400  voidage_rates[iw] = surface_rates[iw] / bw;
401  }
402 
403  // Use average Rs and Rv:
404  auto a = ra.rs;
405  auto b = a;
406  if (io >= 0 && ig >= 0) {
407  b = surface_rates[ig]/(surface_rates[io]+1.0e-15);
408  }
409 
410  double Rs = std::min(a, b);
411 
412  a = ra.rv;
413  b = a;
414  if (io >= 0 && ig >= 0) {
415  b = surface_rates[io]/(surface_rates[ig]+1.0e-15);
416  }
417 
418  double Rv = std::min(a, b);
419 
420  // Determinant of 'R' matrix
421  const double detR = 1.0 - (Rs * Rv);
422 
424  // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
425 
426  const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
427  const double den = bo * detR;
428 
429  voidage_rates[io] = surface_rates[io];
430 
432  voidage_rates[io] -= Rv * surface_rates[ig];
433  }
434 
435  voidage_rates[io] /= den;
436  }
437 
439  // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
440 
441  const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv, 0.0 /*=Rvw*/);
442  const double den = bg * detR;
443 
444  voidage_rates[ig] = surface_rates[ig];
445 
447  voidage_rates[ig] -= Rs * surface_rates[io];
448  }
449 
450  voidage_rates[ig] /= den;
451  }
452  }
453 
454 
467  template <class SolventModule>
468  void
469  calcCoeffSolvent(const RegionId r, const int pvtRegionIdx, double& coeff) const
470  {
471  const auto& ra = attr_.attributes(r);
472  const double p = ra.pressure;
473  const double T = ra.temperature;
474  const auto& solventPvt = SolventModule::solventPvt();
475  const double bs = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
476  coeff = 1.0 / bs;
477  }
478 
479  private:
483  const PhaseUsage phaseUsage_;
484 
488  const RegionMapping<Region> rmap_;
489 
493  struct Attributes {
494  Attributes()
495  : pressure (0.0)
496  , temperature(0.0)
497  , rs(0.0)
498  , rv(0.0)
499  , pv(0.0)
500  , saltConcentration(0.0)
501  {}
502 
503  double pressure;
504  double temperature;
505  double rs;
506  double rv;
507  double pv;
508  double saltConcentration;
509  };
510 
511  RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
512 
513  };
514  } // namespace RateConverter
515 } // namespace Opm
516 
517 #endif /* OPM_RATECONVERTER_HPP_HEADER_INCLUDED */
Convert component rates at surface conditions to phase (voidage) rates at reservoir conditions.
Definition: RateConverter.hpp:68
RegionMapping< Region >::RegionId RegionId
Region identifier.
Definition: RateConverter.hpp:236
SurfaceToReservoirVoidage(const PhaseUsage &phaseUsage, const Region &region)
Constructor.
Definition: RateConverter.hpp:77
void calcCoeffSolvent(const RegionId r, const int pvtRegionIdx, double &coeff) const
Compute coefficients for surface-to-reservoir voidage conversion for solvent.
Definition: RateConverter.hpp:469
void calcReservoirVoidageRates(const RegionId r, const int pvtRegionIdx, const Rates &surface_rates, Rates &voidage_rates) const
Converting surface volume rates to reservoir voidage rates.
Definition: RateConverter.hpp:378
void calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff &coeff) const
Compute coefficients for surface-to-reservoir voidage conversion.
Definition: RateConverter.hpp:267
void defineState(const EbosSimulator &simulator)
Compute pore volume averaged hydrocarbon state pressure, rs and rv.
Definition: RateConverter.hpp:95
int gas(const PhaseUsage &pu)
Numerical ID of active gas phase.
Definition: RegionAttributeHelpers.hpp:394
int oil(const PhaseUsage &pu)
Numerical ID of active oil phase.
Definition: RegionAttributeHelpers.hpp:374
int water(const PhaseUsage &pu)
Numerical ID of active water phase.
Definition: RegionAttributeHelpers.hpp:354
bool water(const PhaseUsage &pu)
Active water predicate.
Definition: RegionAttributeHelpers.hpp:308
bool oil(const PhaseUsage &pu)
Active oil predicate.
Definition: RegionAttributeHelpers.hpp:321
bool gas(const PhaseUsage &pu)
Active gas predicate.
Definition: RegionAttributeHelpers.hpp:334
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition: phaseUsageFromDeck.cpp:37
Definition: BlackoilPhases.hpp:46