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/simulators/wells/RegionAttributeHelpers.hpp>
24
25#include <opm/simulators/utils/BlackoilPhases.hpp>
26#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
27
28#include <dune/grid/common/gridenums.hh>
29#include <dune/grid/common/rangegenerators.hh>
30
31#include <algorithm>
32#include <cassert>
33#include <unordered_map>
34
45namespace Opm {
46 namespace RegionAverageCalculator {
47
59 template <class FluidSystem, class Region>
61 public:
62 using Scalar = typename FluidSystem::Scalar;
71 const Region& region)
72 : phaseUsage_(phaseUsage)
73 , rmap_ (region)
74 , attr_ (rmap_, Attributes())
75 {
76 }
77
78
82 template <typename ElementContext, class Simulator>
83 void defineState(const Simulator& simulator)
84 {
85 int numRegions = 0;
86 const auto& gridView = simulator.gridView();
87 const auto& comm = gridView.comm();
88 for (const auto& reg : rmap_.activeRegions()) {
89 numRegions = std::max(numRegions, reg);
90 }
91 numRegions = comm.max(numRegions);
92 // reg = 0 is used for field
93 for (int reg = 0; reg <= numRegions ; ++ reg) {
94 if (!attr_.has(reg))
95 attr_.insert(reg, Attributes());
96 }
97
98 // quantities for pore volume average
99 std::unordered_map<RegionId, Attributes> attributes_pv;
100
101 // quantities for hydrocarbon volume average
102 std::unordered_map<RegionId, Attributes> attributes_hpv;
103
104 for (int reg = 1; reg <= numRegions ; ++ reg) {
105 attributes_pv.insert({reg, Attributes()});
106 attributes_hpv.insert({reg, Attributes()});
107 }
108
109 ElementContext elemCtx( simulator );
110 OPM_BEGIN_PARALLEL_TRY_CATCH();
111 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
112 elemCtx.updatePrimaryStencil(elem);
113 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
114 const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
115 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
116 const auto& fs = intQuants.fluidState();
117 // use pore volume weighted averages.
118 const Scalar pv_cell =
119 simulator.model().dofTotalVolume(cellIdx)
120 * intQuants.porosity().value();
121
122 // only count oil and gas filled parts of the domain
123 Scalar hydrocarbon = 1.0;
124 const auto& pu = phaseUsage_;
126 hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
127 }
128
129 const int reg = rmap_.region(cellIdx);
130 assert(reg >= 0);
131
132 // sum p, rs, rv, and T.
133 const Scalar hydrocarbonPV = pv_cell*hydrocarbon;
134 if (hydrocarbonPV > 0.) {
135 auto& attr = attributes_hpv[reg];
136 attr.pv += hydrocarbonPV;
138 attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
139 } else {
141 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
142 }
143 }
144
145 if (pv_cell > 0.) {
146 auto& attr = attributes_pv[reg];
147 attr.pv += pv_cell;
149 attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
151 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
152 } else {
154 attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
155 }
156 }
157 }
158 OPM_END_PARALLEL_TRY_CATCH("AverageRegionalPressure::defineState(): ", simulator.vanguard().grid().comm());
159
160 Scalar field_p_hpv = 0.0;
161 Scalar field_hpv = 0.0;
162 Scalar field_p_pv = 0.0;
163 Scalar field_pv = 0.0;
164 for (int reg = 1; reg <= numRegions ; ++ reg) {
165 auto& ra = attr_.attributes(reg);
166 const Scalar hpv_sum = comm.sum(attributes_hpv[reg].pv);
167 // TODO: should we have some epsilon here instead of zero?
168 if (hpv_sum > 0.) {
169 const auto& attri_hpv = attributes_hpv[reg];
170 const Scalar p_hpv_sum = comm.sum(attri_hpv.pressure);
171 ra.pressure = p_hpv_sum / hpv_sum;
172 ra.pv = hpv_sum;
175 } else {
176 // using the pore volume to do the averaging
177 const auto& attri_pv = attributes_pv[reg];
178 const Scalar pv_sum = comm.sum(attri_pv.pv);
179 // pore volums can be zero if a fipnum region is empty
180 if (pv_sum > 0) {
181 const Scalar p_pv_sum = comm.sum(attri_pv.pressure);
182 ra.pressure = p_pv_sum / pv_sum;
183 ra.pv = pv_sum;
185 field_pv += pv_sum;
186 }
187 }
188 }
189 // update field pressure
190 auto& ra_field = attr_.attributes(0);
191 if (field_hpv > 0)
192 ra_field.pressure = field_p_hpv / field_hpv;
193 else if (field_hpv > 0)
194 ra_field.pressure = field_p_pv / field_pv;
195 }
196
202 typedef typename RegionMapping<Region>::RegionId RegionId;
203
208 Scalar
209 pressure(const RegionId r) const
210 {
211 const auto& ra = attr_.attributes(r);
212 return ra.pressure;
213 }
214
215
216 private:
220 const PhaseUsage phaseUsage_;
221
225 const RegionMapping<Region> rmap_;
226
230 struct Attributes {
231 Attributes()
232 : pressure(0.0)
233 , pv(0.0)
234
235 {}
236
237 Scalar pressure;
238 Scalar pv;
239 };
240
241 RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
242
243 };
244 } // namespace RegionAverageCalculator
245} // namespace Opm
246
247#endif /* OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED */
Computes hydrocarbon weighed average pressures over regions.
Definition RegionAverageCalculator.hpp:60
AverageRegionalPressure(const PhaseUsage &phaseUsage, const Region &region)
Constructor.
Definition RegionAverageCalculator.hpp:70
void defineState(const Simulator &simulator)
Compute pore volume averaged hydrocarbon state pressure.
Definition RegionAverageCalculator.hpp:83
RegionMapping< Region >::RegionId RegionId
Region identifier.
Definition RegionAverageCalculator.hpp:202
Scalar pressure(const RegionId r) const
Average pressure.
Definition RegionAverageCalculator.hpp:209
Manages the initializing and running of time dependent problems.
Definition simulator.hh:97
Vanguard & vanguard()
Return a reference to the grid manager of simulation.
Definition simulator.hh:281
const GridView & gridView() const
Return the grid view for which the simulation is done.
Definition simulator.hh:293
Model & model()
Return the physical model used in the simulation.
Definition simulator.hh:299
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 blackoilboundaryratevector.hh:37
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition phaseUsageFromDeck.cpp:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Definition BlackoilPhases.hpp:46