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/linalg/ParallelIstlInformation.hpp>
28 #include <opm/simulators/wells/RegionAttributeHelpers.hpp>
29 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
30 #include <dune/grid/common/gridenums.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  OPM_BEGIN_PARALLEL_TRY_CATCH();
126 
127  const auto& elemEndIt = gridView.template end</*codim=*/0>();
128  for (auto elemIt = gridView.template begin</*codim=*/0>();
129  elemIt != elemEndIt;
130  ++elemIt)
131  {
132 
133  const auto& elem = *elemIt;
134  if (elem.partitionType() != Dune::InteriorEntity)
135  continue;
136 
137  elemCtx.updatePrimaryStencil(elem);
138  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
139  const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
140  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
141  const auto& fs = intQuants.fluidState();
142  // use pore volume weighted averages.
143  const double pv_cell =
144  simulator.model().dofTotalVolume(cellIdx)
145  * intQuants.porosity().value();
146 
147  // only count oil and gas filled parts of the domain
148  double hydrocarbon = 1.0;
149  const auto& pu = phaseUsage_;
151  hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
152  }
153 
154  const int reg = rmap_.region(cellIdx);
155  assert(reg >= 0);
156 
157  // sum p, rs, rv, and T.
158  const double hydrocarbonPV = pv_cell*hydrocarbon;
159  if (hydrocarbonPV > 0.) {
160  auto& attr = attributes_hpv[reg];
161  attr.pv += hydrocarbonPV;
163  attr.rs += fs.Rs().value() * hydrocarbonPV;
164  attr.rv += fs.Rv().value() * hydrocarbonPV;
165  }
167  attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
168  attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
169  } else {
171  attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
172  attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
173  }
174  attr.saltConcentration += fs.saltConcentration().value() * hydrocarbonPV;
175  }
176 
177  if (pv_cell > 0.) {
178  auto& attr = attributes_pv[reg];
179  attr.pv += pv_cell;
181  attr.rs += fs.Rs().value() * pv_cell;
182  attr.rv += fs.Rv().value() * pv_cell;
183  }
185  attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
186  attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * pv_cell;
188  attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
189  attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * pv_cell;
190  } else {
192  attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
193  attr.temperature += fs.temperature(FluidSystem::waterPhaseIdx).value() * pv_cell;
194  }
195  attr.saltConcentration += fs.saltConcentration().value() * pv_cell;
196  }
197  }
198 
199  OPM_END_PARALLEL_TRY_CATCH("SurfaceToReservoirVoidage::defineState() failed: ", simulator.vanguard().grid().comm());
200 
201  for (const auto& reg : rmap_.activeRegions()) {
202  auto& ra = attr_.attributes(reg);
203  const double hpv_sum = comm.sum(attributes_hpv[reg].pv);
204  // TODO: should we have some epsilon here instead of zero?
205  if (hpv_sum > 0.) {
206  const auto& attri_hpv = attributes_hpv[reg];
207  const double p_hpv_sum = comm.sum(attri_hpv.pressure);
208  const double T_hpv_sum = comm.sum(attri_hpv.temperature);
209  const double rs_hpv_sum = comm.sum(attri_hpv.rs);
210  const double rv_hpv_sum = comm.sum(attri_hpv.rv);
211  const double sc_hpv_sum = comm.sum(attri_hpv.saltConcentration);
212 
213  ra.pressure = p_hpv_sum / hpv_sum;
214  ra.temperature = T_hpv_sum / hpv_sum;
215  ra.rs = rs_hpv_sum / hpv_sum;
216  ra.rv = rv_hpv_sum / hpv_sum;
217  ra.pv = hpv_sum;
218  ra.saltConcentration = sc_hpv_sum / hpv_sum;
219  } else {
220  // using the pore volume to do the averaging
221  const auto& attri_pv = attributes_pv[reg];
222  const double pv_sum = comm.sum(attri_pv.pv);
223  assert(pv_sum > 0.);
224  const double p_pv_sum = comm.sum(attri_pv.pressure);
225  const double T_pv_sum = comm.sum(attri_pv.temperature);
226  const double rs_pv_sum = comm.sum(attri_pv.rs);
227  const double rv_pv_sum = comm.sum(attri_pv.rv);
228  const double sc_pv_sum = comm.sum(attri_pv.saltConcentration);
229 
230  ra.pressure = p_pv_sum / pv_sum;
231  ra.temperature = T_pv_sum / pv_sum;
232  ra.rs = rs_pv_sum / pv_sum;
233  ra.rv = rv_pv_sum / pv_sum;
234  ra.pv = pv_sum;
235  ra.saltConcentration = sc_pv_sum / pv_sum;
236  }
237  }
238  }
239 
245  typedef typename RegionMapping<Region>::RegionId RegionId;
246 
274  template <class Coeff>
275  void
276  calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
277  {
278  const auto& pu = phaseUsage_;
279  const auto& ra = attr_.attributes(r);
280  const double p = ra.pressure;
281  const double T = ra.temperature;
282  const double saltConcentration = ra.saltConcentration;
283 
284  const int iw = RegionAttributeHelpers::PhasePos::water(pu);
285  const int io = RegionAttributeHelpers::PhasePos::oil (pu);
286  const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
287 
288  std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
289 
291  // q[w]_r = q[w]_s / bw
292 
293  const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
294 
295  coeff[iw] = 1.0 / bw;
296  }
297 
298  // Actual Rs and Rv:
299  double Rs = ra.rs;
300  double Rv = ra.rv;
301 
302  // Determinant of 'R' matrix
303  const double detR = 1.0 - (Rs * Rv);
304 
306  // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
307 
308  const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
309  const double den = bo * detR;
310 
311  coeff[io] += 1.0 / den;
312 
314  coeff[ig] -= ra.rv / den;
315  }
316  }
317 
319  // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
320 
321  const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
322  const double den = bg * detR;
323 
324  coeff[ig] += 1.0 / den;
325 
327  coeff[io] -= ra.rs / den;
328  }
329  }
330  }
331 
332  template <class Coeff>
333  void
334  calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
335  {
336  const auto& pu = phaseUsage_;
337  const auto& ra = attr_.attributes(r);
338  const double p = ra.pressure;
339  const double T = ra.temperature;
340  const double saltConcentration = ra.saltConcentration;
341 
342  const int iw = RegionAttributeHelpers::PhasePos::water(pu);
343  const int io = RegionAttributeHelpers::PhasePos::oil (pu);
344  const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
345 
346  std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
347 
349  // q[w]_r = q[w]_s / bw
350 
351  const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
352 
353  coeff[iw] = 1.0 / bw;
354  }
355 
357  const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0);
358  coeff[io] += 1.0 / bo;
359  }
360 
362  const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0);
363  coeff[ig] += 1.0 / bg;
364  }
365  }
366 
367 
385  template <class Rates >
386  void
387  calcReservoirVoidageRates(const RegionId r, const int pvtRegionIdx, const Rates& surface_rates,
388  Rates& voidage_rates) const
389  {
390  assert(voidage_rates.size() == surface_rates.size());
391 
392  std::fill(voidage_rates.begin(), voidage_rates.end(), 0.0);
393 
394  const auto& pu = phaseUsage_;
395  const auto& ra = attr_.attributes(r);
396  const double p = ra.pressure;
397  const double T = ra.temperature;
398  const double saltConcentration = ra.saltConcentration;
399 
400  const int iw = RegionAttributeHelpers::PhasePos::water(pu);
401  const int io = RegionAttributeHelpers::PhasePos::oil (pu);
402  const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
403 
405  // q[w]_r = q[w]_s / bw
406 
407  const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
408 
409  voidage_rates[iw] = surface_rates[iw] / bw;
410  }
411 
412  // Use average Rs and Rv:
413  auto a = ra.rs;
414  auto b = a;
415  if (io >= 0 && ig >= 0) {
416  b = surface_rates[ig]/(surface_rates[io]+1.0e-15);
417  }
418 
419  double Rs = std::min(a, b);
420 
421  a = ra.rv;
422  b = a;
423  if (io >= 0 && ig >= 0) {
424  b = surface_rates[io]/(surface_rates[ig]+1.0e-15);
425  }
426 
427  double Rv = std::min(a, b);
428 
429  // Determinant of 'R' matrix
430  const double detR = 1.0 - (Rs * Rv);
431 
433  // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
434 
435  const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
436  const double den = bo * detR;
437 
438  voidage_rates[io] = surface_rates[io];
439 
441  voidage_rates[io] -= Rv * surface_rates[ig];
442  }
443 
444  voidage_rates[io] /= den;
445  }
446 
448  // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
449 
450  const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
451  const double den = bg * detR;
452 
453  voidage_rates[ig] = surface_rates[ig];
454 
456  voidage_rates[ig] -= Rs * surface_rates[io];
457  }
458 
459  voidage_rates[ig] /= den;
460  }
461  }
462 
463 
476  template <class SolventModule>
477  void
478  calcCoeffSolvent(const RegionId r, const int pvtRegionIdx, double& coeff) const
479  {
480  const auto& ra = attr_.attributes(r);
481  const double p = ra.pressure;
482  const double T = ra.temperature;
483  const auto& solventPvt = SolventModule::solventPvt();
484  const double bs = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
485  coeff = 1.0 / bs;
486  }
487 
488  private:
492  const PhaseUsage phaseUsage_;
493 
497  const RegionMapping<Region> rmap_;
498 
502  struct Attributes {
503  Attributes()
504  : pressure (0.0)
505  , temperature(0.0)
506  , rs(0.0)
507  , rv(0.0)
508  , pv(0.0)
509  , saltConcentration(0.0)
510  {}
511 
512  double pressure;
513  double temperature;
514  double rs;
515  double rv;
516  double pv;
517  double saltConcentration;
518  };
519 
520  RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
521 
522  };
523  } // namespace RateConverter
524 } // namespace Opm
525 
526 #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:245
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:478
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:387
void calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff &coeff) const
Compute coefficients for surface-to-reservoir voidage conversion.
Definition: RateConverter.hpp:276
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:395
int oil(const PhaseUsage &pu)
Numerical ID of active oil phase.
Definition: RegionAttributeHelpers.hpp:375
int water(const PhaseUsage &pu)
Numerical ID of active water phase.
Definition: RegionAttributeHelpers.hpp:355
bool water(const PhaseUsage &pu)
Active water predicate.
Definition: RegionAttributeHelpers.hpp:309
bool oil(const PhaseUsage &pu)
Active oil predicate.
Definition: RegionAttributeHelpers.hpp:322
bool gas(const PhaseUsage &pu)
Active gas predicate.
Definition: RegionAttributeHelpers.hpp:335
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:33
Definition: BlackoilPhases.hpp:46