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
24
27
28#include <dune/grid/common/gridenums.hh>
29#include <dune/grid/common/rangegenerators.hh>
30
31#include <algorithm>
32#include <cassert>
33#include <numeric>
34#include <unordered_map>
35
46namespace Opm {
47 namespace RegionAverageCalculator {
48
60 template <class FluidSystem, class Region>
62 public:
63 using Scalar = typename FluidSystem::Scalar;
72 const Region& region)
73 : phaseUsage_(phaseUsage)
74 , rmap_ (region)
75 , attr_ (rmap_, Attributes())
76 {
77 }
78
79
83 template <typename ElementContext, class Simulator>
84 void defineState(const Simulator& simulator)
85 {
86 const auto& gridView = simulator.gridView();
87 const auto& comm = gridView.comm();
88 int numRegions =
89 std::accumulate(rmap_.activeRegions().begin(),
90 rmap_.activeRegions().end(), 0,
91 [](const auto acc, const auto& reg)
92 { return std::max(acc, reg); });
93 numRegions = comm.max(numRegions);
94 // reg = 0 is used for field
95 for (int reg = 0; reg <= numRegions ; ++ reg) {
96 if (!attr_.has(reg))
97 attr_.insert(reg, Attributes());
98 }
99
100 // quantities for pore volume average
101 std::unordered_map<RegionId, Attributes> attributes_pv;
102
103 // quantities for hydrocarbon volume average
104 std::unordered_map<RegionId, Attributes> attributes_hpv;
105
106 for (int reg = 1; reg <= numRegions ; ++ reg) {
107 attributes_pv.insert({reg, Attributes()});
108 attributes_hpv.insert({reg, Attributes()});
109 }
110
111 ElementContext elemCtx( simulator );
113 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
114 elemCtx.updatePrimaryStencil(elem);
115 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
116 const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
117 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
118 const auto& fs = intQuants.fluidState();
119 // use pore volume weighted averages.
120 const Scalar pv_cell =
121 simulator.model().dofTotalVolume(cellIdx)
122 * intQuants.porosity().value();
123
124 // only count oil and gas filled parts of the domain
125 Scalar hydrocarbon = 1.0;
126 const auto& pu = phaseUsage_;
128 hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
129 }
130
131 const int reg = rmap_.region(cellIdx);
132 assert(reg >= 0);
133
134 // sum p, rs, rv, and T.
135 const Scalar hydrocarbonPV = pv_cell*hydrocarbon;
136 if (hydrocarbonPV > 0.) {
137 auto& attr = attributes_hpv[reg];
138 attr.pv += hydrocarbonPV;
140 attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
141 } else {
143 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
144 }
145 }
146
147 if (pv_cell > 0.) {
148 auto& attr = attributes_pv[reg];
149 attr.pv += pv_cell;
151 attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
153 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
154 } else {
156 attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
157 }
158 }
159 }
160 OPM_END_PARALLEL_TRY_CATCH("AverageRegionalPressure::defineState(): ", simulator.vanguard().grid().comm());
161
162 Scalar field_p_hpv = 0.0;
163 Scalar field_hpv = 0.0;
164 Scalar field_p_pv = 0.0;
165 Scalar field_pv = 0.0;
166 for (int reg = 1; reg <= numRegions ; ++ reg) {
167 auto& ra = attr_.attributes(reg);
168 const Scalar hpv_sum = comm.sum(attributes_hpv[reg].pv);
169 // TODO: should we have some epsilon here instead of zero?
170 if (hpv_sum > 0.) {
171 const auto& attri_hpv = attributes_hpv[reg];
172 const Scalar p_hpv_sum = comm.sum(attri_hpv.pressure);
173 ra.pressure = p_hpv_sum / hpv_sum;
174 ra.pv = hpv_sum;
175 field_p_hpv += p_hpv_sum;
176 field_hpv += hpv_sum;
177 } else {
178 // using the pore volume to do the averaging
179 const auto& attri_pv = attributes_pv[reg];
180 const Scalar pv_sum = comm.sum(attri_pv.pv);
181 // pore volums can be zero if a fipnum region is empty
182 if (pv_sum > 0) {
183 const Scalar p_pv_sum = comm.sum(attri_pv.pressure);
184 ra.pressure = p_pv_sum / pv_sum;
185 ra.pv = pv_sum;
186 field_p_pv += p_pv_sum;
187 field_pv += pv_sum;
188 }
189 }
190 }
191 // update field pressure
192 auto& ra_field = attr_.attributes(0);
193 if (field_hpv > 0)
194 ra_field.pressure = field_p_hpv / field_hpv;
195 else if (field_hpv > 0)
196 ra_field.pressure = field_p_pv / field_pv;
197 }
198
205
210 Scalar
211 pressure(const RegionId r) const
212 {
213 const auto& ra = attr_.attributes(r);
214 return ra.pressure;
215 }
216
217
218 private:
222 const PhaseUsage phaseUsage_;
223
227 const RegionMapping<Region> rmap_;
228
232 struct Attributes {
233 Attributes()
234 : pressure(0.0)
235 , pv(0.0)
236
237 {}
238
240 Scalar pv;
241 };
242
243 RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
244
245 };
246 } // namespace RegionAverageCalculator
247} // namespace Opm
248
249#endif /* OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED */
#define OPM_END_PARALLEL_TRY_CATCH(prefix, comm)
Catch exception and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:192
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:158
Definition: RegionAverageCalculator.hpp:61
typename FluidSystem::Scalar Scalar
Definition: RegionAverageCalculator.hpp:63
AverageRegionalPressure(const PhaseUsage &phaseUsage, const Region &region)
Definition: RegionAverageCalculator.hpp:71
void defineState(const Simulator &simulator)
Definition: RegionAverageCalculator.hpp:84
RegionMapping< Region >::RegionId RegionId
Definition: RegionAverageCalculator.hpp:204
Scalar pressure(const RegionId r) const
Definition: RegionAverageCalculator.hpp:211
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:84
Vanguard & vanguard()
Return a reference to the grid manager of simulation.
Definition: simulator.hh:234
const GridView & gridView() const
Return the grid view for which the simulation is done.
Definition: simulator.hh:246
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:252
bool water(const PhaseUsage &pu)
Definition: RegionAttributeHelpers.hpp:309
bool oil(const PhaseUsage &pu)
Definition: RegionAttributeHelpers.hpp:322
bool gas(const PhaseUsage &pu)
Definition: RegionAttributeHelpers.hpp:335
Definition: blackoilboundaryratevector.hh:39
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
int RegionId
Definition: WellConstraints.hpp:37
Definition: BlackoilPhases.hpp:46