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
26
27#include <dune/grid/common/gridenums.hh>
28#include <dune/grid/common/rangegenerators.hh>
29
30#include <algorithm>
31#include <cassert>
32#include <numeric>
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;
70 explicit AverageRegionalPressure(const Region& region)
71 : rmap_ (region)
72 , attr_ (rmap_, Attributes())
73 {
74 }
75
76
80 template <typename ElementContext, class Simulator>
81 void defineState(const Simulator& simulator)
82 {
83 const auto& gridView = simulator.gridView();
84 const auto& comm = gridView.comm();
85 int numRegions =
86 std::accumulate(rmap_.activeRegions().begin(),
87 rmap_.activeRegions().end(), 0,
88 [](const auto acc, const auto& reg)
89 { return std::max(acc, reg); });
90 numRegions = comm.max(numRegions);
91 // reg = 0 is used for field
92 for (int reg = 0; reg <= numRegions ; ++ reg) {
93 if (!attr_.has(reg))
94 attr_.insert(reg, Attributes());
95 }
96
97 // quantities for pore volume average
98 std::unordered_map<RegionId, Attributes> attributes_pv;
99
100 // quantities for hydrocarbon volume average
101 std::unordered_map<RegionId, Attributes> attributes_hpv;
102
103 for (int reg = 1; reg <= numRegions ; ++ reg) {
104 attributes_pv.insert({reg, Attributes()});
105 attributes_hpv.insert({reg, Attributes()});
106 }
107
108 ElementContext elemCtx( simulator );
110 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
111 elemCtx.updatePrimaryStencil(elem);
112 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
113 const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
114 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
115 const auto& fs = intQuants.fluidState();
116 // use pore volume weighted averages.
117 const Scalar pv_cell =
118 simulator.model().dofTotalVolume(cellIdx)
119 * intQuants.porosity().value();
120
121 // only count oil and gas filled parts of the domain
122 Scalar hydrocarbon = 1.0;
123 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
124 hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
125 }
126
127 const int reg = rmap_.region(cellIdx);
128 assert(reg >= 0);
129
130 // sum p, rs, rv, and T.
131 const Scalar hydrocarbonPV = pv_cell*hydrocarbon;
132 if (hydrocarbonPV > 0.) {
133 auto& attr = attributes_hpv[reg];
134 attr.pv += hydrocarbonPV;
135 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
136 attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
137 } else {
138 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
139 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
140 }
141 }
142
143 if (pv_cell > 0.) {
144 auto& attr = attributes_pv[reg];
145 attr.pv += pv_cell;
146 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
147 attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
148 } else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
149 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
150 } else {
151 assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
152 attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
153 }
154 }
155 }
156 OPM_END_PARALLEL_TRY_CATCH("AverageRegionalPressure::defineState(): ", simulator.vanguard().grid().comm());
157
158 Scalar field_p_hpv = 0.0;
159 Scalar field_hpv = 0.0;
160 Scalar field_p_pv = 0.0;
161 Scalar field_pv = 0.0;
162 for (int reg = 1; reg <= numRegions ; ++ reg) {
163 auto& ra = attr_.attributes(reg);
164 const Scalar hpv_sum = comm.sum(attributes_hpv[reg].pv);
165 // TODO: should we have some epsilon here instead of zero?
166 if (hpv_sum > 0.) {
167 const auto& attri_hpv = attributes_hpv[reg];
168 const Scalar p_hpv_sum = comm.sum(attri_hpv.pressure);
169 ra.pressure = p_hpv_sum / hpv_sum;
170 ra.pv = hpv_sum;
171 field_p_hpv += p_hpv_sum;
172 field_hpv += hpv_sum;
173 } else {
174 // using the pore volume to do the averaging
175 const auto& attri_pv = attributes_pv[reg];
176 const Scalar pv_sum = comm.sum(attri_pv.pv);
177 // pore volums can be zero if a fipnum region is empty
178 if (pv_sum > 0) {
179 const Scalar p_pv_sum = comm.sum(attri_pv.pressure);
180 ra.pressure = p_pv_sum / pv_sum;
181 ra.pv = pv_sum;
182 field_p_pv += p_pv_sum;
183 field_pv += pv_sum;
184 }
185 }
186 }
187 // update field pressure
188 auto& ra_field = attr_.attributes(0);
189 if (field_hpv > 0)
190 ra_field.pressure = field_p_hpv / field_hpv;
191 else if (field_hpv > 0)
192 ra_field.pressure = field_p_pv / field_pv;
193 }
194
201
206 Scalar
207 pressure(const RegionId r) const
208 {
209 const auto& ra = attr_.attributes(r);
210 return ra.pressure;
211 }
212
213
214 private:
218 const RegionMapping<Region> rmap_;
219
223 struct Attributes {
224 Attributes()
225 : pressure(0.0)
226 , pv(0.0)
227
228 {}
229
231 Scalar pv;
232 };
233
234 RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
235
236 };
237 } // namespace RegionAverageCalculator
238} // namespace Opm
239
240#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:60
AverageRegionalPressure(const Region &region)
Definition: RegionAverageCalculator.hpp:70
typename FluidSystem::Scalar Scalar
Definition: RegionAverageCalculator.hpp:62
void defineState(const Simulator &simulator)
Definition: RegionAverageCalculator.hpp:81
RegionMapping< Region >::RegionId RegionId
Definition: RegionAverageCalculator.hpp:200
Scalar pressure(const RegionId r) const
Definition: RegionAverageCalculator.hpp:207
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
Definition: blackoilboundaryratevector.hh:39
int RegionId
Definition: WellConstraints.hpp:37