OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

BlackoilModelEbos: clean up and fix the FIP code

Browse Source 
mainly this should now work properly in parallel, because non-interior
cells are not counted multiple times anymore. also, the number of
loops over the global arrays has been reduced, some variables have
been renamed and some comments were added.

finally this fixes the average pressure for regions that do not
contain hydrocarbons (or at least it unifies it with the approach for
regions that contain hydrocarbons).
This commit is contained in:
Andreas Lauser 2017-03-13 10:04:31 +01:00
parent 3e1758f5c8
commit 5194cc5122
1 changed files with 70 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

175
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -1012,17 +1012,25 @@ namespace Opm {
std::vector<std::vector<double> >
computeFluidInPlace(const std::vector<int>& fipnum) const
{
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const auto& comm = grid_.comm();
const auto& gridView = grid_.leafGridView();
const int nc = gridView.size(/*codim=*/0);
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
int ntFip = *std::max_element(fipnum.begin(), fipnum.end());
ntFip = comm.max(ntFip);
std::vector<double> tpv(ntFip, 0.0);
std::vector<double> hcpv(ntFip, 0.0);
std::vector<std::vector<double> > regionValues(ntFip, std::vector<double>(FIPDataType::fipValues,0.0));
for (int i = 0; i<FIPDataType::fipValues; i++) {
fip_.fip[i].resize(nc,0.0);
}
ElementContext elemCtx(ebosSimulator_);
const auto& elemEndIt = elemCtx.gridView().template end</*codim=*/0>();
for (auto elemIt = elemCtx.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
@ -1038,6 +1046,12 @@ namespace Opm {
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const int regionIdx = fipnum[cellIdx] - 1;
if (regionIdx < 0) {
// the given cell is not attributed to any region
continue;
}
// calculate the pore volume of the current cell. Note that the porosity
// returned by the intensive quantities is defined as the ratio of pore
// space to total cell volume and includes all pressure dependent (->
@ -1054,86 +1068,33 @@ namespace Opm {
const double s = fs.saturation(flowPhaseToEbosPhaseIdx(phase)).value();
fip_.fip[phase][cellIdx] = b * s * pv;
if (active_[ phase ]) {
regionValues[regionIdx][phase] += fip_.fip[phase][cellIdx];
}
}
if (active_[ Oil ] && active_[ Gas ]) {
// Account for gas dissolved in oil and vaporized oil
fip_.fip[FIPDataType::FIP_DISSOLVED_GAS][cellIdx] = fs.Rs().value() * fip_.fip[FIPDataType::FIP_LIQUID][cellIdx];
fip_.fip[FIPDataType::FIP_VAPORIZED_OIL][cellIdx] = fs.Rv().value() * fip_.fip[FIPDataType::FIP_VAPOUR][cellIdx];
regionValues[regionIdx][FIPData::FIP_DISSOLVED_GAS] += fip_.fip[FIPData::FIP_DISSOLVED_GAS][cellIdx];
regionValues[regionIdx][FIPData::FIP_VAPORIZED_OIL] += fip_.fip[FIPData::FIP_VAPORIZED_OIL][cellIdx];
}
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
tpv[regionIdx] += pv;
hcpv[regionIdx] += pv * hydrocarbon;
}
// For a parallel run this is just a local maximum and needs to be updated later
const auto& comm = grid_.comm();
int dims = *std::max_element(fipnum.begin(), fipnum.end());
dims = comm.max(dims);
std::vector<std::vector<double>> values(dims, std::vector<double>(FIPDataType::fipValues,0.0));
//Accumulate phases for each region
for (int phase = 0; phase < maxnp; ++phase) {
if (active_[ phase ]) {
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1) {
values[region][phase] += fip_.fip[phase][c];
}
}
}
}
//Accumulate RS and RV-volumes for each region
if (active_[ Oil ] && active_[ Gas ]) {
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1) {
values[region][FIPData::FIP_DISSOLVED_GAS] += fip_.fip[FIPData::FIP_DISSOLVED_GAS][c];
values[region][FIPData::FIP_VAPORIZED_OIL] += fip_.fip[FIPData::FIP_VAPORIZED_OIL][c];
}
}
}
std::vector<double> hcpv(dims, 0.0);
std::vector<double> pres(dims, 0.0);
for (auto elemIt = elemCtx.gridView().template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const int region = fipnum[cellIdx] - 1;
if (region != -1) {
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
// calculate the pore volume of the current cell. Note that the
// porosity returned by the intensive quantities is defined as the
// ratio of pore space to total cell volume and includes all pressure
// dependent (-> rock compressibility) and static modifiers (MULTPV,
// MULTREGP, NTG, PORV, MINPV and friends). Also note that because of
// this, the porosity returned by the intensive quantities can be
// outside of the physical range [0, 1] in pathetic cases.
const double pv =
ebosSimulator_.model().dofTotalVolume(cellIdx)
* intQuants.porosity().value();
hcpv[region] += pv * hydrocarbon;
pres[region] += pv * fs.pressure(FluidSystem::oilPhaseIdx).value();
}
}
// sum tpv (-> total pore volume of the regions) and hcpv (-> pore volume of the
// the regions that is occupied by hydrocarbons)
comm.sum(tpv.data(), tpv.size());
comm.sum(hcpv.data(), hcpv.size());
comm.sum(pres.data(), pres.size());
for (auto elemIt = elemCtx.gridView().template begin</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
@ -1146,43 +1107,47 @@ namespace Opm {
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const int region = fipnum[cellIdx] - 1;
if (region != -1) {
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
// calculate the pore volume of the current cell. Note that the
// porosity returned by the intensive quantities is defined as the
// ratio of pore space to total cell volume and includes all pressure
// dependent (-> rock compressibility) and static modifiers (MULTPV,
// MULTREGP, NTG, PORV, MINPV and friends). Also note that because of
// this, the porosity returned by the intensive quantities can be
// outside of the physical range [0, 1] in pathetic cases.
const double pv =
ebosSimulator_.model().dofTotalVolume(cellIdx)
* intQuants.porosity().value();
fip_.fip[FIPDataType::FIP_PV][cellIdx] = pv;
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
//Compute hydrocarbon pore volume weighted average pressure.
//If we have no hydrocarbon in region, use pore volume weighted average pressure instead
if (hcpv[region] != 0.0) {
fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[region];
} else {
fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pres[region] / pv;
}
values[region][FIPDataType::FIP_PV] += fip_.fip[FIPDataType::FIP_PV][cellIdx];
values[region][FIPDataType::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx];
const int regionIdx = fipnum[cellIdx] - 1;
if (regionIdx < 0) {
// the cell is not attributed to any region. ignore it!
continue;
}
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
// calculate the pore volume of the current cell. Note that the
// porosity returned by the intensive quantities is defined as the
// ratio of pore space to total cell volume and includes all pressure
// dependent (-> rock compressibility) and static modifiers (MULTPV,
// MULTREGP, NTG, PORV, MINPV and friends). Also note that because of
// this, the porosity returned by the intensive quantities can be
// outside of the physical range [0, 1] in pathetic cases.
const double pv =
ebosSimulator_.model().dofTotalVolume(cellIdx)
* intQuants.porosity().value();
fip_.fip[FIPDataType::FIP_PV][cellIdx] = pv;
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
//Compute hydrocarbon pore volume weighted average pressure.
//If we have no hydrocarbon in region, use pore volume weighted average pressure instead
if (hcpv[regionIdx] > 1e-10) {
fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[regionIdx];
} else {
fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() / tpv[regionIdx];
}
regionValues[regionIdx][FIPDataType::FIP_PV] += fip_.fip[FIPDataType::FIP_PV][cellIdx];
regionValues[regionIdx][FIPDataType::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx];
}
for(int reg=0; reg < dims; ++reg) {
comm.sum(values[reg].data(), values[reg].size());
// sum the results over all processes
for(int regionIdx=0; regionIdx < ntFip; ++regionIdx) {
comm.sum(regionValues[regionIdx].data(), regionValues[regionIdx].size());
}
return values;
return regionValues;
}
SimulationDataContainer getSimulatorData () const