20#ifndef OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED
21#define OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED
28#include <dune/grid/common/gridenums.hh>
29#include <dune/grid/common/rangegenerators.hh>
34#include <unordered_map>
47 namespace RegionAverageCalculator {
60 template <
class Flu
idSystem,
class Region>
63 using Scalar =
typename FluidSystem::Scalar;
75 , attr_ (rmap_, Attributes())
83 template <
typename ElementContext,
class Simulator>
86 const auto& gridView = simulator.
gridView();
87 const auto& comm = gridView.comm();
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);
95 for (
int reg = 0; reg <= numRegions ; ++ reg) {
97 attr_.insert(reg, Attributes());
101 std::unordered_map<RegionId, Attributes> attributes_pv;
104 std::unordered_map<RegionId, Attributes> attributes_hpv;
106 for (
int reg = 1; reg <= numRegions ; ++ reg) {
107 attributes_pv.insert({reg, Attributes()});
108 attributes_hpv.insert({reg, Attributes()});
111 ElementContext elemCtx( simulator );
113 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
114 elemCtx.updatePrimaryStencil(elem);
115 elemCtx.updatePrimaryIntensiveQuantities(0);
116 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
117 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
118 const auto& fs = intQuants.fluidState();
121 simulator.
model().dofTotalVolume(cellIdx)
122 * intQuants.porosity().value();
126 const auto& pu = phaseUsage_;
128 hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
131 const int reg = rmap_.region(cellIdx);
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;
143 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
148 auto& attr = attributes_pv[reg];
151 attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
153 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
156 attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
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);
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;
175 field_p_hpv += p_hpv_sum;
176 field_hpv += hpv_sum;
179 const auto& attri_pv = attributes_pv[reg];
180 const Scalar pv_sum = comm.sum(attri_pv.pv);
183 const Scalar p_pv_sum = comm.sum(attri_pv.pressure);
184 ra.pressure = p_pv_sum / pv_sum;
186 field_p_pv += p_pv_sum;
192 auto& ra_field = attr_.attributes(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;
213 const auto& ra = attr_.attributes(r);
227 const RegionMapping<Region> rmap_;
243 RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
#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 ®ion)
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