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>
49namespace 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
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
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
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