My Project
Loading...
Searching...
No Matches
RegionAverageCalculator.hpp
Go to the documentation of this file.
1/*
2 Copyright 2021, Equinor
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 3 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
20#ifndef OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED
21#define OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED
22
23#include <opm/core/props/BlackoilPhases.hpp>
24#include <opm/simulators/wells/RegionAttributeHelpers.hpp>
25#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
26
27#include <dune/grid/common/gridenums.hh>
28#include <dune/grid/common/rangegenerators.hh>
29
30#include <algorithm>
31#include <cassert>
32#include <cmath>
33#include <memory>
34#include <stdexcept>
35#include <type_traits>
36#include <unordered_map>
37#include <utility>
38#include <vector>
39
50namespace Opm {
51 namespace RegionAverageCalculator {
52
64 template <class FluidSystem, class Region>
66 public:
75 const Region& region)
76 : phaseUsage_(phaseUsage)
77 , rmap_ (region)
78 , attr_ (rmap_, Attributes())
79 {
80 }
81
82
86 template <typename ElementContext, class EbosSimulator>
87 void defineState(const EbosSimulator& simulator)
88 {
89
90
91 int numRegions = 0;
92 const auto& gridView = simulator.gridView();
93 const auto& comm = gridView.comm();
94 for (const auto& reg : rmap_.activeRegions()) {
95 numRegions = std::max(numRegions, reg);
96 }
97 numRegions = comm.max(numRegions);
98 for (int reg = 1; reg <= numRegions ; ++ reg) {
99 if (!attr_.has(reg))
100 attr_.insert(reg, Attributes());
101 }
102 // create map from cell to region
103 // and set all attributes to zero
104 for (int reg = 1; reg <= numRegions ; ++ reg) {
105 auto& ra = attr_.attributes(reg);
106 ra.pressure = 0.0;
107 ra.pv = 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 (int reg = 1; reg <= numRegions ; ++ reg) {
118 attributes_pv.insert({reg, Attributes()});
119 attributes_hpv.insert({reg, Attributes()});
120 }
121
122 ElementContext elemCtx( simulator );
123
124 OPM_BEGIN_PARALLEL_TRY_CATCH();
125 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
126 elemCtx.updatePrimaryStencil(elem);
127 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
128 const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
129 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
130 const auto& fs = intQuants.fluidState();
131 // use pore volume weighted averages.
132 const double pv_cell =
133 simulator.model().dofTotalVolume(cellIdx)
134 * intQuants.porosity().value();
135
136 // only count oil and gas filled parts of the domain
137 double hydrocarbon = 1.0;
138 const auto& pu = phaseUsage_;
140 hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
141 }
142
143 const int reg = rmap_.region(cellIdx);
144 assert(reg >= 0);
145
146 // sum p, rs, rv, and T.
147 const double hydrocarbonPV = pv_cell*hydrocarbon;
148 if (hydrocarbonPV > 0.) {
149 auto& attr = attributes_hpv[reg];
150 attr.pv += hydrocarbonPV;
152 attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
153 } else {
155 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
156 }
157 }
158
159 if (pv_cell > 0.) {
160 auto& attr = attributes_pv[reg];
161 attr.pv += pv_cell;
163 attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
165 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
166 } else {
168 attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
169 }
170 }
171 }
172 OPM_END_PARALLEL_TRY_CATCH("AverageRegionalPressure::defineState(): ", simulator.vanguard().grid().comm());
173
174 for (int reg = 1; reg <= numRegions ; ++ reg) {
175 auto& ra = attr_.attributes(reg);
176 const double hpv_sum = comm.sum(attributes_hpv[reg].pv);
177 // TODO: should we have some epsilon here instead of zero?
178 if (hpv_sum > 0.) {
179 const auto& attri_hpv = attributes_hpv[reg];
180 const double p_hpv_sum = comm.sum(attri_hpv.pressure);
181 ra.pressure = p_hpv_sum / hpv_sum;
182 } else {
183 // using the pore volume to do the averaging
184 const auto& attri_pv = attributes_pv[reg];
185 const double pv_sum = comm.sum(attri_pv.pv);
186 // pore volums can be zero if a fipnum region is empty
187 if (pv_sum > 0) {
188 const double p_pv_sum = comm.sum(attri_pv.pressure);
189 ra.pressure = p_pv_sum / pv_sum;
190 }
191 }
192 }
193 }
194
201
206 double
207 pressure(const RegionId r) const
208 {
209 if (r == 0 ) // region 0 is the whole field
210 {
211 double pressure = 0.0;
212 int num_active_regions = 0;
213 for (const auto& attr : attr_.attributes()) {
214 const auto& value = *attr.second;
215 const auto& ra = value.attr_;
216 pressure += ra.pressure;
218 }
220 }
221
222 const auto& ra = attr_.attributes(r);
223 return ra.pressure;
224 }
225
226
227 private:
231 const PhaseUsage phaseUsage_;
232
236 const RegionMapping<Region> rmap_;
237
241 struct Attributes {
242 Attributes()
243 : pressure(0.0)
244 , pv(0.0)
245
246 {}
247
248 double pressure;
249 double pv;
250
251 };
252
253 RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
254
255 };
256 } // namespace RegionAverageCalculator
257} // namespace Opm
258
259#endif /* OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED */
Definition AquiferInterface.hpp:35
Computes hydrocarbon weighed average pressures over regions.
Definition RegionAverageCalculator.hpp:65
void defineState(const EbosSimulator &simulator)
Compute pore volume averaged hydrocarbon state pressure, *.
Definition RegionAverageCalculator.hpp:87
AverageRegionalPressure(const PhaseUsage &phaseUsage, const Region &region)
Constructor.
Definition RegionAverageCalculator.hpp:74
double pressure(const RegionId r) const
Average pressure.
Definition RegionAverageCalculator.hpp:207
RegionMapping< Region >::RegionId RegionId
Region identifier.
Definition RegionAverageCalculator.hpp:200
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