20#ifndef OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED
21#define OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED
27#include <dune/grid/common/gridenums.hh>
28#include <dune/grid/common/rangegenerators.hh>
33#include <unordered_map>
46 namespace RegionAverageCalculator {
59 template <
class Flu
idSystem,
class Region>
62 using Scalar =
typename FluidSystem::Scalar;
72 , attr_ (rmap_, Attributes())
80 template <
typename ElementContext,
class Simulator>
83 const auto& gridView = simulator.
gridView();
84 const auto& comm = gridView.comm();
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);
92 for (
int reg = 0; reg <= numRegions ; ++ reg) {
94 attr_.insert(reg, Attributes());
98 std::unordered_map<RegionId, Attributes> attributes_pv;
101 std::unordered_map<RegionId, Attributes> attributes_hpv;
103 for (
int reg = 1; reg <= numRegions ; ++ reg) {
104 attributes_pv.insert({reg, Attributes()});
105 attributes_hpv.insert({reg, Attributes()});
108 ElementContext elemCtx( simulator );
110 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
111 elemCtx.updatePrimaryStencil(elem);
112 elemCtx.updatePrimaryIntensiveQuantities(0);
113 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
114 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
115 const auto& fs = intQuants.fluidState();
118 simulator.
model().dofTotalVolume(cellIdx)
119 * intQuants.porosity().value();
123 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
124 hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
127 const int reg = rmap_.region(cellIdx);
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;
138 assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
139 attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
144 auto& attr = attributes_pv[reg];
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;
151 assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
152 attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
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);
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;
171 field_p_hpv += p_hpv_sum;
172 field_hpv += hpv_sum;
175 const auto& attri_pv = attributes_pv[reg];
176 const Scalar pv_sum = comm.sum(attri_pv.pv);
179 const Scalar p_pv_sum = comm.sum(attri_pv.pressure);
180 ra.pressure = p_pv_sum / pv_sum;
182 field_p_pv += p_pv_sum;
188 auto& ra_field = attr_.attributes(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;
209 const auto& ra = attr_.attributes(r);
218 const RegionMapping<Region> rmap_;
234 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:60
AverageRegionalPressure(const Region ®ion)
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