Merge pull request #1095 from andlaus/improve_FIP

Fixes for the fluid in place code
This commit is contained in:
Atgeirr Flø Rasmussen 2017-03-15 11:14:33 +01:00 committed by GitHub
commit 345f10bc18
2 changed files with 128 additions and 165 deletions

View File

@ -633,7 +633,6 @@ namespace Opm {
double step = dsMax()/maxVal;
step = std::min(step, 1.0);
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
if (active_[Water]) {
double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]];
@ -669,7 +668,6 @@ namespace Opm {
double& rs = reservoir_state.gasoilratio()[cell_idx];
double& rv = reservoir_state.rv()[cell_idx];
// phase translation sg <-> rs
const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx];
const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
@ -706,9 +704,6 @@ namespace Opm {
break;
}
const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
if (rs > ( rsSat * (1+epsilon) ) ) {
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
@ -727,6 +722,7 @@ namespace Opm {
//std::cout << "watonly rv -> sg" << cell_idx << std::endl;
break;
}
const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
if (rv > rvSat * (1+epsilon) ) {
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
@ -1016,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)
{
@ -1042,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 (->
@ -1058,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)
{
@ -1150,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
@ -1382,19 +1343,16 @@ namespace Opm {
return fip_;
}
const Simulator& ebosSimulator() const
{ return ebosSimulator_; }
protected:
protected:
const ISTLSolverType& istlSolver() const
{
assert( istlSolver_ );
return *istlSolver_;
}
// --------- Data members ---------
Simulator& ebosSimulator_;
@ -1428,10 +1386,6 @@ namespace Opm {
BVector dx_old_;
mutable FIPDataType fip_;
// --------- Protected methods ---------
public:
/// return the StandardWells object
@ -1714,10 +1668,8 @@ namespace Opm {
std::cout << "equation scaling not suported yet" << std::endl;
//updateEquationsScaling();
}
}
double dpMaxRel() const { return param_.dp_max_rel_; }
double dsMax() const { return param_.ds_max_; }
double drMaxRel() const { return param_.dr_max_rel_; }

View File

@ -250,6 +250,9 @@ public:
auto solver = createSolver(well_model);
std::vector<std::vector<double>> currentFluidInPlace;
std::vector<double> currentFluidInPlaceTotals;
// Compute orignal fluid in place if this has not been done yet
if (originalFluidInPlace.empty()) {
solver->model().convertInput(/*iterationIdx=*/0, state, ebosSimulator_ );
@ -259,6 +262,9 @@ public:
originalFluidInPlaceTotals = FIPTotals(originalFluidInPlace, state);
FIPUnitConvert(eclState().getUnits(), originalFluidInPlace);
FIPUnitConvert(eclState().getUnits(), originalFluidInPlaceTotals);
currentFluidInPlace = originalFluidInPlace;
currentFluidInPlaceTotals = originalFluidInPlaceTotals;
}
// write the inital state at the report stage
@ -266,6 +272,13 @@ public:
Dune::Timer perfTimer;
perfTimer.start();
if (terminal_output_) {
outputFluidInPlace(originalFluidInPlaceTotals, currentFluidInPlaceTotals,eclState().getUnits(), 0);
for (size_t reg = 0; reg < originalFluidInPlace.size(); ++reg) {
outputFluidInPlace(originalFluidInPlace[reg], currentFluidInPlace[reg], eclState().getUnits(), reg+1);
}
}
// No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating
output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
@ -335,9 +348,9 @@ public:
++timer;
// Compute current fluid in place.
std::vector<std::vector<double>> currentFluidInPlace;
currentFluidInPlace = solver->computeFluidInPlace(fipnum);
std::vector<double> currentFluidInPlaceTotals = FIPTotals(currentFluidInPlace, state);
currentFluidInPlaceTotals = FIPTotals(currentFluidInPlace, state);
const std::string version = moduleVersionName();
FIPUnitConvert(eclState().getUnits(), currentFluidInPlace);
@ -639,56 +652,54 @@ protected:
totals[i] += fip[reg][i];
}
}
const int numCells = Opm::AutoDiffGrid::numCells(grid());
const auto& pv = geo_.poreVolume();
const auto& gridView = ebosSimulator_.gridManager().gridView();
const auto& comm = gridView.comm();
double pv_hydrocarbon_sum = 0.0;
double p_pv_hydrocarbon_sum = 0.0;
if ( ! is_parallel_run_ )
ElementContext elemCtx(ebosSimulator_);
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const double& p = fs.pressure(FluidSystem::oilPhaseIdx).value();
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
totals[5] += pv[cellIdx];
pv_hydrocarbon_sum += pv[cellIdx] * hydrocarbon;
p_pv_hydrocarbon_sum += p * pv[cellIdx] * hydrocarbon;
}
}
else
{
#if HAVE_MPI
const auto & pinfo =
boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
// Mask with 1 for owned cell and 0 otherwise
const auto& mask = pinfo.updateOwnerMask(pv);
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const double& p = fs.pressure(FluidSystem::oilPhaseIdx).value();
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
if( mask[cellIdx] )
{
totals[5] += pv[cellIdx];
pv_hydrocarbon_sum += pv[cellIdx] * hydrocarbon;
p_pv_hydrocarbon_sum += p * pv[cellIdx] * hydrocarbon;
}
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
totals[5] = pinfo.communicator().sum(totals[5]);
pv_hydrocarbon_sum = pinfo.communicator().sum(pv_hydrocarbon_sum);
p_pv_hydrocarbon_sum= pinfo.communicator().sum(p_pv_hydrocarbon_sum);
#else
OPM_THROW(std::logic_error, "Requested a parallel run without MPI available!");
#endif
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();
const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
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();
totals[5] += pv;
pv_hydrocarbon_sum += pv*hydrocarbon;
p_pv_hydrocarbon_sum += p*pv*hydrocarbon;
}
pv_hydrocarbon_sum = comm.sum(pv_hydrocarbon_sum);
p_pv_hydrocarbon_sum = comm.sum(p_pv_hydrocarbon_sum);
totals[5] = comm.sum(totals[5]);
totals[6] = (p_pv_hydrocarbon_sum / pv_hydrocarbon_sum);
return totals;
}
@ -734,7 +745,7 @@ protected:
<< std::fixed << std::setprecision(0)
<< " : PORV =" << std::setw(14) << cip[5] << " RB :\n";
if (!reg) {
ss << " : Pressure is weighted by hydrocarbon pore voulme :\n"
ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
<< " : Pore volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";