Loading...
Searching...
No Matches
RegionAverageCalculator.hpp
Go to the documentation of this file.
/*
Copyright 2021, Equinor
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED
#define OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED
#include <opm/simulators/wells/RegionAttributeHelpers.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <dune/grid/common/gridenums.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <algorithm>
#include <cassert>
#include <numeric>
#include <unordered_map>
namespace Opm {
namespace RegionAverageCalculator {
template <class FluidSystem, class Region>
class AverageRegionalPressure {
public:
using Scalar = typename FluidSystem::Scalar;
explicit AverageRegionalPressure(const Region& region)
: rmap_ (region)
, attr_ (rmap_, Attributes())
{
}
template <typename ElementContext, class Simulator>
void defineState(const Simulator& simulator)
{
const auto& gridView = simulator.gridView();
const auto& comm = gridView.comm();
int numRegions =
std::accumulate(rmap_.activeRegions().begin(),
rmap_.activeRegions().end(), 0,
[](const auto acc, const auto& reg)
{ return std::max(acc, reg); });
numRegions = comm.max(numRegions);
// reg = 0 is used for field
for (int reg = 0; reg <= numRegions ; ++ reg) {
if (!attr_.has(reg))
attr_.insert(reg, Attributes());
}
// quantities for pore volume average
std::unordered_map<RegionId, Attributes> attributes_pv;
// quantities for hydrocarbon volume average
std::unordered_map<RegionId, Attributes> attributes_hpv;
for (int reg = 1; reg <= numRegions ; ++ reg) {
attributes_pv.insert({reg, Attributes()});
attributes_hpv.insert({reg, Attributes()});
}
ElementContext elemCtx( simulator );
OPM_BEGIN_PARALLEL_TRY_CATCH();
for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
// use pore volume weighted averages.
const Scalar pv_cell =
simulator.model().dofTotalVolume(cellIdx)
* intQuants.porosity().value();
// only count oil and gas filled parts of the domain
Scalar hydrocarbon = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
}
const int reg = rmap_.region(cellIdx);
assert(reg >= 0);
// sum p, rs, rv, and T.
const Scalar hydrocarbonPV = pv_cell*hydrocarbon;
if (hydrocarbonPV > 0.) {
auto& attr = attributes_hpv[reg];
attr.pv += hydrocarbonPV;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
} else {
assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
}
}
if (pv_cell > 0.) {
auto& attr = attributes_pv[reg];
attr.pv += pv_cell;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
} else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
} else {
assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
}
}
}
OPM_END_PARALLEL_TRY_CATCH("AverageRegionalPressure::defineState(): ", simulator.vanguard().grid().comm());
Scalar field_p_hpv = 0.0;
Scalar field_hpv = 0.0;
Scalar field_p_pv = 0.0;
Scalar field_pv = 0.0;
for (int reg = 1; reg <= numRegions ; ++ reg) {
auto& ra = attr_.attributes(reg);
const Scalar hpv_sum = comm.sum(attributes_hpv[reg].pv);
// TODO: should we have some epsilon here instead of zero?
if (hpv_sum > 0.) {
const auto& attri_hpv = attributes_hpv[reg];
const Scalar p_hpv_sum = comm.sum(attri_hpv.pressure);
ra.pressure = p_hpv_sum / hpv_sum;
ra.pv = hpv_sum;
field_p_hpv += p_hpv_sum;
field_hpv += hpv_sum;
} else {
// using the pore volume to do the averaging
const auto& attri_pv = attributes_pv[reg];
const Scalar pv_sum = comm.sum(attri_pv.pv);
// pore volums can be zero if a fipnum region is empty
if (pv_sum > 0) {
const Scalar p_pv_sum = comm.sum(attri_pv.pressure);
ra.pressure = p_pv_sum / pv_sum;
ra.pv = pv_sum;
field_p_pv += p_pv_sum;
field_pv += pv_sum;
}
}
}
// update field pressure
auto& ra_field = attr_.attributes(0);
if (field_hpv > 0)
ra_field.pressure = field_p_hpv / field_hpv;
else if (field_hpv > 0)
ra_field.pressure = field_p_pv / field_pv;
}
typedef typename RegionMapping<Region>::RegionId RegionId;
Scalar
pressure(const RegionId r) const
{
const auto& ra = attr_.attributes(r);
return ra.pressure;
}
private:
const RegionMapping<Region> rmap_;
struct Attributes {
Attributes()
: pressure(0.0)
, pv(0.0)
{}
Scalar pressure;
Scalar pv;
};
RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
};
} // namespace RegionAverageCalculator
} // namespace Opm
#endif /* OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED */
Opm::RegionAverageCalculator::AverageRegionalPressure::AverageRegionalPressure
AverageRegionalPressure(const Region &region)
Constructor.
Definition RegionAverageCalculator.hpp:70
Opm::RegionAverageCalculator::AverageRegionalPressure::defineState
void defineState(const Simulator &simulator)
Compute pore volume averaged hydrocarbon state pressure.
Definition RegionAverageCalculator.hpp:81
Opm::RegionAverageCalculator::AverageRegionalPressure::RegionId
RegionMapping< Region >::RegionId RegionId
Region identifier.
Definition RegionAverageCalculator.hpp:200
Opm::RegionAverageCalculator::AverageRegionalPressure::pressure
Scalar pressure(const RegionId r) const
Average pressure.
Definition RegionAverageCalculator.hpp:207
Opm::Simulator
Manages the initializing and running of time dependent problems.
Definition simulator.hh:84
Opm::Simulator::vanguard
Vanguard & vanguard()
Return a reference to the grid manager of simulation.
Definition simulator.hh:234
Opm::Simulator::gridView
const GridView & gridView() const
Return the grid view for which the simulation is done.
Definition simulator.hh:246
Opm::Simulator::model
Model & model()
Return the physical model used in the simulation.
Definition simulator.hh:252
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilbioeffectsmodules.hh:43