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

Merge pull request #1025 from andlaus/flow_ebos-fix-FIP

Browse Source 
Flow ebos fix fip
This commit is contained in:
Atgeirr Flø Rasmussen 2017-02-21 14:31:21 +01:00 committed by GitHub
commit 2797eca463
2 changed files with 131 additions and 145 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

239
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -1006,8 +1006,6 @@ namespace Opm {
{ {
using namespace Opm::AutoDiffGrid; using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_); const int nc = numCells(grid_);
//const ADB pv_mult = poroMult(pressure);
const auto& pv = geo_.poreVolume();
const int maxnp = Opm::BlackoilPhases::MaxNumPhases; const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
for (int i = 0; i<FIPDataType::fipValues; i++) { for (int i = 0; i<FIPDataType::fipValues; i++) {
@ -1015,9 +1013,11 @@ namespace Opm {
} }
ElementContext elemCtx(ebosSimulator_); ElementContext elemCtx(ebosSimulator_);
auto elemIt = elemCtx.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = elemCtx.gridView().template end</*codim=*/0>(); const auto& elemEndIt = elemCtx.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) { for (auto elemIt = elemCtx.gridView().template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
const auto& elem = *elemIt; const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity) { if (elem.partitionType() != Dune::InteriorEntity) {
continue; continue;
@ -1030,11 +1030,22 @@ namespace Opm {
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState(); 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();
for (int phase = 0; phase < maxnp; ++phase) { for (int phase = 0; phase < maxnp; ++phase) {
const double b = fs.invB(flowPhaseToEbosPhaseIdx(phase)).value(); const double b = fs.invB(flowPhaseToEbosPhaseIdx(phase)).value();
const double s = fs.saturation(flowPhaseToEbosPhaseIdx(phase)).value(); const double s = fs.saturation(flowPhaseToEbosPhaseIdx(phase)).value();
const double pv_mult = 1.0; //todo
fip_.fip[phase][cellIdx] = pv_mult * b * s * pv[cellIdx]; fip_.fip[phase][cellIdx] = b * s * pv;
} }
if (active_[ Oil ] && active_[ Gas ]) { if (active_[ Oil ] && active_[ Gas ]) {
@ -1045,150 +1056,122 @@ namespace Opm {
} }
// For a parallel run this is just a local maximum and needs to be updated later // 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()); 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)); std::vector<std::vector<double>> values(dims, std::vector<double>(FIPDataType::fipValues,0.0));
std::vector<double> hcpv(dims, 0.0); //Accumulate phases for each region
std::vector<double> pres(dims, 0.0); for (int phase = 0; phase < maxnp; ++phase) {
if (active_[ phase ]) {
if ( !isParallel() )
{
//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) { for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1; const int region = fipnum[c] - 1;
if (region != -1) { if (region != -1) {
values[region][FIPDataType::FIP_DISSOLVED_GAS] += fip_.fip[FIPDataType::FIP_DISSOLVED_GAS][c];
values[region][FIPDataType::FIP_VAPORIZED_OIL] += fip_.fip[FIPDataType::FIP_VAPORIZED_OIL][c];
}
}
}
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1) {
const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(c, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
hcpv[region] += pv[c] * hydrocarbon;
pres[region] += pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value();
}
}
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1) {
fip_.fip[FIPDataType::FIP_PV][c] = pv[c];
const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(c, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
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) {
fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][c] = pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[region];
} else {
fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
}
values[region][FIPDataType::FIP_PV] += fip_.fip[FIPDataType::FIP_PV][c];
values[region][FIPDataType::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][c];
}
}
}
else
{
#if HAVE_MPI
// mask[c] is 1 if we need to compute something in parallel
const auto & pinfo =
boost::any_cast<const ParallelISTLInformation&>(istlSolver().parallelInformation());
const auto& mask = pinfo.updateOwnerMask( fipnum );
auto comm = pinfo.communicator();
// Compute the global dims value and resize values accordingly.
dims = comm.max(dims);
values.resize(dims, std::vector<double>(FIPDataType::fipValues,0.0));
//Accumulate phases for each region
for (int phase = 0; phase < maxnp; ++phase) {
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1 && mask[c]) {
values[region][phase] += fip_.fip[phase][c]; values[region][phase] += fip_.fip[phase][c];
} }
} }
} }
}
//Accumulate RS and RV-volumes for each region //Accumulate RS and RV-volumes for each region
if (active_[ Oil ] && active_[ Gas ]) { if (active_[ Oil ] && active_[ Gas ]) {
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1 && mask[c]) {
values[region][FIPDataType::FIP_DISSOLVED_GAS] += fip_.fip[FIPDataType::FIP_DISSOLVED_GAS][c];
values[region][FIPDataType::FIP_VAPORIZED_OIL] += fip_.fip[FIPDataType::FIP_VAPORIZED_OIL][c];
}
}
}
hcpv = std::vector<double>(dims, 0.0);
pres = std::vector<double>(dims, 0.0);
for (int c = 0; c < nc; ++c) { for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1; const int region = fipnum[c] - 1;
if (region != -1 && mask[c]) { if (region != -1) {
const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(c, /*timeIdx=*/0); values[region][FIPData::FIP_DISSOLVED_GAS] += fip_.fip[FIPData::FIP_DISSOLVED_GAS][c];
const auto& fs = intQuants.fluidState(); values[region][FIPData::FIP_VAPORIZED_OIL] += fip_.fip[FIPData::FIP_VAPORIZED_OIL][c];
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
hcpv[region] += pv[c] * hydrocarbon;
pres[region] += pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value();
} }
} }
}
comm.sum(hcpv.data(), hcpv.size()); std::vector<double> hcpv(dims, 0.0);
comm.sum(pres.data(), pres.size()); std::vector<double> pres(dims, 0.0);
for (int c = 0; c < nc; ++c) { for (auto elemIt = elemCtx.gridView().template begin</*codim=*/0>();
const int region = fipnum[c] - 1; elemIt != elemEndIt;
if (region != -1 && mask[c]) { ++elemIt)
fip_.fip[FIPDataType::FIP_PV][c] = pv[c]; {
const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(c, /*timeIdx=*/0); const auto& elem = *elemIt;
const auto& fs = intQuants.fluidState(); if (elem.partitionType() != Dune::InteriorEntity) {
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value(); continue;
}
if (hcpv[region] != 0) { elemCtx.updatePrimaryStencil(elem);
fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][c] = pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[region]; elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
} else {
fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
}
values[region][FIPDataType::FIP_PV] += fip_.fip[FIPDataType::FIP_PV][c]; unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
values[region][FIPDataType::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][c]; 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();
}
}
comm.sum(hcpv.data(), hcpv.size());
comm.sum(pres.data(), pres.size());
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();
// 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;
} }
}
// For the frankenstein branch we hopefully can turn values into a vanilla values[region][FIPDataType::FIP_PV] += fip_.fip[FIPDataType::FIP_PV][cellIdx];
// std::vector<double>, use some index magic above, use one communication values[region][FIPDataType::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx];
// to sum up the vector entries instead of looping over the regions.
for(int reg=0; reg < dims; ++reg)
{
comm.sum(values[reg].data(), values[reg].size());
} }
#else }
// This should never happen!
OPM_THROW(std::logic_error, "HAVE_MPI should be defined if we are running in parallel"); for(int reg=0; reg < dims; ++reg) {
#endif comm.sum(values[reg].data(), values[reg].size());
} }
return values; return values;

37
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -191,7 +191,6 @@ public:
SimulatorReport report; SimulatorReport report;
SimulatorReport stepReport; SimulatorReport stepReport;
bool ooip_computed = false;
std::vector<int> fipnum_global = eclState().get3DProperties().getIntGridProperty("FIPNUM").getData(); std::vector<int> fipnum_global = eclState().get3DProperties().getIntGridProperty("FIPNUM").getData();
//Get compressed cell fipnum. //Get compressed cell fipnum.
std::vector<int> fipnum(Opm::UgGridHelpers::numCells(grid())); std::vector<int> fipnum(Opm::UgGridHelpers::numCells(grid()));
@ -202,7 +201,9 @@ public:
fipnum[c] = fipnum_global[Opm::UgGridHelpers::globalCell(grid())[c]]; fipnum[c] = fipnum_global[Opm::UgGridHelpers::globalCell(grid())[c]];
} }
} }
std::vector<std::vector<double>> OOIP; std::vector<std::vector<double>> originalFluidInPlace;
std::vector<double> originalFluidInPlaceTotals;
// Main simulation loop. // Main simulation loop.
while (!timer.done()) { while (!timer.done()) {
// Report timestep. // Report timestep.
@ -256,11 +257,15 @@ public:
auto solver = createSolver(well_model); auto solver = createSolver(well_model);
// Compute orignal FIP; // Compute orignal fluid in place if this has not been done yet
if (!ooip_computed) { if (originalFluidInPlace.empty()) {
OOIP = solver->computeFluidInPlace(fipnum); solver->model().convertInput(/*iterationIdx=*/0, state, ebosSimulator_ );
FIPUnitConvert(eclState().getUnits(), OOIP); ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
ooip_computed = true;
originalFluidInPlace = solver->computeFluidInPlace(fipnum);
originalFluidInPlaceTotals = FIPTotals(originalFluidInPlace, state);
FIPUnitConvert(eclState().getUnits(), originalFluidInPlace);
FIPUnitConvert(eclState().getUnits(), originalFluidInPlaceTotals);
} }
if( terminal_output_ ) if( terminal_output_ )
@ -315,20 +320,18 @@ public:
report.solver_time += solver_timer.secsSinceStart(); report.solver_time += solver_timer.secsSinceStart();
// Compute current fluid in place. // Compute current fluid in place.
std::vector<std::vector<double>> COIP; std::vector<std::vector<double>> currentFluidInPlace;
COIP = solver->computeFluidInPlace(fipnum); currentFluidInPlace = solver->computeFluidInPlace(fipnum);
std::vector<double> OOIP_totals = FIPTotals(OOIP, state); std::vector<double> currentFluidInPlaceTotals = FIPTotals(currentFluidInPlace, state);
std::vector<double> COIP_totals = FIPTotals(COIP, state);
FIPUnitConvert(eclState().getUnits(), COIP); FIPUnitConvert(eclState().getUnits(), currentFluidInPlace);
FIPUnitConvert(eclState().getUnits(), OOIP_totals); FIPUnitConvert(eclState().getUnits(), currentFluidInPlaceTotals);
FIPUnitConvert(eclState().getUnits(), COIP_totals);
if (terminal_output_ ) if (terminal_output_ )
{ {
outputFluidInPlace(OOIP_totals, COIP_totals,eclState().getUnits(), 0); outputFluidInPlace(originalFluidInPlaceTotals, currentFluidInPlaceTotals,eclState().getUnits(), 0);
for (size_t reg = 0; reg < OOIP.size(); ++reg) { for (size_t reg = 0; reg < originalFluidInPlace.size(); ++reg) {
outputFluidInPlace(OOIP[reg], COIP[reg], eclState().getUnits(), reg+1); outputFluidInPlace(originalFluidInPlace[reg], currentFluidInPlace[reg], eclState().getUnits(), reg+1);
} }
std::string msg; std::string msg;