My Project
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/linalg/ParallelIstlInformation.hpp>
26 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
27 
28 #include <dune/grid/common/gridenums.hh>
29 #include <algorithm>
30 #include <cmath>
31 #include <memory>
32 #include <stdexcept>
33 #include <type_traits>
34 #include <unordered_map>
35 #include <utility>
36 #include <vector>
47 namespace Opm {
48  namespace RegionAverageCalculator {
49 
61  template <class FluidSystem, class Region>
63  public:
72  const Region& region)
73  : phaseUsage_(phaseUsage)
74  , rmap_ (region)
75  , attr_ (rmap_, Attributes())
76  {
77  }
78 
79 
83  template <typename ElementContext, class EbosSimulator>
84  void defineState(const EbosSimulator& simulator)
85  {
86 
87 
88  int numRegions = 0;
89  const auto& gridView = simulator.gridView();
90  const auto& comm = gridView.comm();
91  for (const auto& reg : rmap_.activeRegions()) {
92  numRegions = std::max(numRegions, reg);
93  }
94  numRegions = comm.max(numRegions);
95  for (int reg = 1; reg <= numRegions ; ++ reg) {
96  if(!attr_.has(reg))
97  attr_.insert(reg, Attributes());
98  }
99  // create map from cell to region
100  // and set all attributes to zero
101  for (int reg = 1; reg <= numRegions ; ++ reg) {
102  auto& ra = attr_.attributes(reg);
103  ra.pressure = 0.0;
104  ra.pv = 0.0;
105 
106  }
107 
108  // quantities for pore volume average
109  std::unordered_map<RegionId, Attributes> attributes_pv;
110 
111  // quantities for hydrocarbon volume average
112  std::unordered_map<RegionId, Attributes> attributes_hpv;
113 
114  for (int reg = 1; reg <= numRegions ; ++ reg) {
115  attributes_pv.insert({reg, Attributes()});
116  attributes_hpv.insert({reg, Attributes()});
117  }
118 
119  ElementContext elemCtx( simulator );
120  const auto& elemEndIt = gridView.template end</*codim=*/0>();
121  OPM_BEGIN_PARALLEL_TRY_CATCH();
122 
123  for (auto elemIt = gridView.template begin</*codim=*/0>();
124  elemIt != elemEndIt;
125  ++elemIt)
126  {
127 
128  const auto& elem = *elemIt;
129  if (elem.partitionType() != Dune::InteriorEntity)
130  continue;
131 
132  elemCtx.updatePrimaryStencil(elem);
133  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
134  const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
135  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
136  const auto& fs = intQuants.fluidState();
137  // use pore volume weighted averages.
138  const double pv_cell =
139  simulator.model().dofTotalVolume(cellIdx)
140  * intQuants.porosity().value();
141 
142  // only count oil and gas filled parts of the domain
143  double hydrocarbon = 1.0;
144  const auto& pu = phaseUsage_;
146  hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
147  }
148 
149  const int reg = rmap_.region(cellIdx);
150  assert(reg >= 0);
151 
152  // sum p, rs, rv, and T.
153  const double hydrocarbonPV = pv_cell*hydrocarbon;
154  if (hydrocarbonPV > 0.) {
155  auto& attr = attributes_hpv[reg];
156  attr.pv += hydrocarbonPV;
158  attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
159  } else {
161  attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
162  }
163  }
164 
165  if (pv_cell > 0.) {
166  auto& attr = attributes_pv[reg];
167  attr.pv += pv_cell;
169  attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
171  attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
172  } else {
174  attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
175  }
176  }
177  }
178  OPM_END_PARALLEL_TRY_CATCH("AverageRegionalPressure::defineState(): ", simulator.vanguard().grid().comm());
179 
180  for (int reg = 1; reg <= numRegions ; ++ reg) {
181  auto& ra = attr_.attributes(reg);
182  const double hpv_sum = comm.sum(attributes_hpv[reg].pv);
183  // TODO: should we have some epsilon here instead of zero?
184  if (hpv_sum > 0.) {
185  const auto& attri_hpv = attributes_hpv[reg];
186  const double p_hpv_sum = comm.sum(attri_hpv.pressure);
187  ra.pressure = p_hpv_sum / hpv_sum;
188  } else {
189  // using the pore volume to do the averaging
190  const auto& attri_pv = attributes_pv[reg];
191  const double pv_sum = comm.sum(attri_pv.pv);
192  // pore volums can be zero if a fipnum region is empty
193  if (pv_sum > 0) {
194  const double p_pv_sum = comm.sum(attri_pv.pressure);
195  ra.pressure = p_pv_sum / pv_sum;
196  }
197  }
198  }
199  }
200 
206  typedef typename RegionMapping<Region>::RegionId RegionId;
207 
212  double
213  pressure(const RegionId r) const
214  {
215  if (r == 0 ) // region 0 is the whole field
216  {
217  double pressure = 0.0;
218  int num_active_regions = 0;
219  for (const auto& attr : attr_.attributes()) {
220  const auto& value = *attr.second;
221  const auto& ra = value.attr_;
222  pressure += ra.pressure;
223  num_active_regions ++;
224  }
225  return pressure / num_active_regions;
226  }
227 
228  const auto& ra = attr_.attributes(r);
229  return ra.pressure;
230  }
231 
232 
233  private:
237  const PhaseUsage phaseUsage_;
238 
242  const RegionMapping<Region> rmap_;
243 
247  struct Attributes {
248  Attributes()
249  : pressure(0.0)
250  , pv(0.0)
251 
252  {}
253 
254  double pressure;
255  double pv;
256 
257  };
258 
259  RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
260 
261  };
262  } // namespace RegionAverageCalculator
263 } // namespace Opm
264 
265 #endif /* OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED */
Computes hydrocarbon weighed average pressures over regions.
Definition: RegionAverageCalculator.hpp:62
void defineState(const EbosSimulator &simulator)
Compute pore volume averaged hydrocarbon state pressure, *.
Definition: RegionAverageCalculator.hpp:84
AverageRegionalPressure(const PhaseUsage &phaseUsage, const Region &region)
Constructor.
Definition: RegionAverageCalculator.hpp:71
double pressure(const RegionId r) const
Average pressure.
Definition: RegionAverageCalculator.hpp:213
RegionMapping< Region >::RegionId RegionId
Region identifier.
Definition: RegionAverageCalculator.hpp:206
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