Output FIP in flow_ebos

This commit is contained in:
Tor Harald Sandve 2016-11-08 08:46:42 +01:00
parent 20353c7f94
commit 5a917a4828
5 changed files with 392 additions and 28 deletions

View File

@ -131,6 +131,19 @@ namespace Opm {
typedef ISTLSolver< MatrixBlockType, VectorBlockType > ISTLSolverType;
//typedef typename SolutionVector :: value_type PrimaryVariables ;
struct FIPData {
enum FipId {
FIP_AQUA = Opm::Water,
FIP_LIQUID = Opm::Oil,
FIP_VAPOUR = Opm::Gas,
FIP_DISSOLVED_GAS = 3,
FIP_VAPORIZED_OIL = 4,
FIP_PV = 5, //< Pore volume
FIP_WEIGHTED_PRESSURE = 6
};
std::array<std::vector<double>, 7> fip;
};
// --------- Public methods ---------
/// Construct the model. It will retain references to the
@ -188,6 +201,26 @@ namespace Opm {
}
}
bool
isParallel() const
{
#if HAVE_MPI
if ( linsolver_.parallelInformation().type() !=
typeid(ParallelISTLInformation) )
{
return false;
}
else
{
const auto& comm =boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation()).communicator();
return comm.size() > 1;
}
#else
return false;
#endif
}
const EclipseState& eclState() const
{ return *ebosSimulator_.gridManager().eclState(); }
@ -937,13 +970,191 @@ namespace Opm {
}
std::vector<std::vector<double> >
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum) const
computeFluidInPlace(const std::vector<int>& fipnum) const
{
OPM_THROW(std::logic_error,
"computeFluidInPlace() not implemented by BlackoilModelEbos!");
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
//const ADB pv_mult = poroMult(pressure);
const auto& pv = geo_.poreVolume();
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
for (int i = 0; i<7; i++) {
fip_.fip[i].resize(nc,0.0);
}
for (int c = 0; c < nc; ++c) {
const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(c, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
for (int phase = 0; phase < maxnp; ++phase) {
const double& b = fs.invB(flowPhaseToEbosPhaseIdx(phase)).value();
const double& s = fs.saturation(flowPhaseToEbosPhaseIdx(phase)).value();
const double pv_mult = 1.0; //todo
fip_.fip[phase][c] = pv_mult * b * s * pv[c];
}
if (active_[ Oil ] && active_[ Gas ]) {
// Account for gas dissolved in oil and vaporized oil
fip_.fip[FIPData::FIP_DISSOLVED_GAS][c] = fs.Rs().value() * fip_.fip[FIPData::FIP_LIQUID][c];
fip_.fip[FIPData::FIP_VAPORIZED_OIL][c] = fs.Rv().value() * fip_.fip[FIPData::FIP_VAPOUR][c];
}
}
// For a parallel run this is just a local maximum and needs to be updated later
int dims = *std::max_element(fipnum.begin(), fipnum.end());
std::vector<std::vector<double>> values(dims, std::vector<double>(7,0.0));
std::vector<double> hcpv(dims, 0.0);
std::vector<double> pres(dims, 0.0);
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) {
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];
}
}
}
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[FIPData::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[FIPData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[region];
} else {
fip_.fip[FIPData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
}
values[region][FIPData::FIP_PV] += fip_.fip[FIPData::FIP_PV][c];
values[region][FIPData::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPData::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&>(linsolver_.parallelInformation());
const auto& mask = pinfo.getOwnerMask();
auto comm = pinfo.communicator();
// Compute the global dims value and resize values accordingly.
dims = comm.max(dims);
values.resize(dims, std::vector<double>(7,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];
}
}
}
//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 && mask[c]) {
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];
}
}
}
hcpv = V::Zero(dims);
pres = V::Zero(dims);
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1 && mask[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();
hcpv[region] += pv[c] * hydrocarbon;
pres[region] += pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value();
}
}
comm.sum(hcpv.data(), hcpv.size());
comm.sum(pres.data(), pres.size());
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1 && mask[c]) {
fip_.fip[FIPData::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();
if (hcpv[region] != 0) {
fip_.fip[FIPData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[region];
} else {
fip_.fip[FIPData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
}
values[region][FIPData::FIP_PV] += fip_.fip[FIPData::FIP_PV][c];
values[region][FIPData::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPData::FIP_WEIGHTED_PRESSURE][c];
}
}
// For the frankenstein branch we hopefully can turn values into a vanilla
// std::vector<double>, use some index magic above, use one communication
// 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");
#endif
}
return values;
}
const FIPData& getFIPData() const {
return fip_;
}
const Simulator& ebosSimulator() const
{ return ebosSimulator_; }
@ -983,6 +1194,7 @@ namespace Opm {
std::vector<std::vector<double>> residual_norms_history_;
double current_relaxation_;
BVector dx_old_;
mutable FIPData fip_;

View File

@ -143,6 +143,12 @@ namespace Opm {
return model_->computeFluidInPlace(x, fipnum);
}
std::vector<std::vector<double>>
computeFluidInPlace(const std::vector<int>& fipnum) const
{
return model_->computeFluidInPlace(fipnum);
}
/// Reference to physical model.
const PhysicalModel& model() const;

View File

@ -192,6 +192,19 @@ public:
std::vector<double> well_potentials;
DynamicListEconLimited dynamic_list_econ_limited;
bool ooip_computed = false;
std::vector<int> fipnum_global = eclState().get3DProperties().getIntGridProperty("FIPNUM").getData();
//Get compressed cell fipnum.
std::vector<int> fipnum(Opm::UgGridHelpers::numCells(grid()));
if (fipnum_global.empty()) {
std::fill(fipnum.begin(), fipnum.end(), 0);
} else {
for (size_t c = 0; c < fipnum.size(); ++c) {
fipnum[c] = fipnum_global[Opm::UgGridHelpers::globalCell(grid())[c]];
}
}
std::vector<std::vector<double>> OOIP;
// Main simulation loop.
while (!timer.done()) {
// Report timestep.
@ -235,6 +248,13 @@ public:
auto solver = createSolver(well_model);
// Compute orignal FIP;
if (!timer.initialStep() && !ooip_computed) {
OOIP = solver->computeFluidInPlace(fipnum);
FIPUnitConvert(eclState().getUnits(), OOIP);
ooip_computed = true;
}
// write the inital state at the report stage
if (timer.initialStep()) {
//output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
@ -294,6 +314,21 @@ public:
// Report timing.
const double st = solver_timer.secsSinceStart();
// Compute current FIP.
std::vector<std::vector<double>> COIP;
COIP = solver->computeFluidInPlace(fipnum);
FIPUnitConvert(eclState().getUnits(), COIP);
std::vector<double> OOIP_totals = FIPTotals(OOIP, state);
std::vector<double> COIP_totals = FIPTotals(COIP, state);
if ( terminal_output_ )
{
outputFluidInPlace(OOIP_totals, COIP_totals,eclState().getUnits(), 0);
for (size_t reg = 0; reg < OOIP.size(); ++reg) {
outputFluidInPlace(OOIP[reg], COIP[reg], eclState().getUnits(), reg+1);
}
}
// accumulate total time
stime += st;
@ -559,6 +594,106 @@ protected:
well_state, list_econ_limited);
}
void FIPUnitConvert(const UnitSystem& units,
std::vector<std::vector<double>>& fip)
{
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
for (size_t i = 0; i < fip.size(); ++i) {
fip[i][0] = unit::convert::to(fip[i][0], unit::stb);
fip[i][1] = unit::convert::to(fip[i][1], unit::stb);
fip[i][2] = unit::convert::to(fip[i][2], 1000*unit::cubic(unit::feet));
fip[i][3] = unit::convert::to(fip[i][3], 1000*unit::cubic(unit::feet));
fip[i][4] = unit::convert::to(fip[i][4], unit::stb);
fip[i][5] = unit::convert::to(fip[i][5], unit::stb);
fip[i][6] = unit::convert::to(fip[i][6], unit::psia);
}
}
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
for (size_t i = 0; i < fip.size(); ++i) {
fip[i][6] = unit::convert::to(fip[i][6], unit::barsa);
}
}
}
std::vector<double> FIPTotals(const std::vector<std::vector<double>>& fip, const ReservoirState& state)
{
std::vector<double> totals(6,0.0);
for (int i = 0; i < 5; ++i) {
for (size_t reg = 0; reg < fip.size(); ++reg) {
totals[i] += fip[reg][i];
}
}
const int numCells = Opm::AutoDiffGrid::numCells(grid());
const auto& pv = geo_.poreVolume();
double pv_hydrocarbon_sum = 0.0;
double p_pv_hydrocarbon_sum = 0.0;
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 ( ! is_parallel_run_ )
{
totals[5] += pv[cellIdx];
pv_hydrocarbon_sum += pv[cellIdx] * hydrocarbon;
p_pv_hydrocarbon_sum += p * pv[cellIdx] * hydrocarbon;
}
else {
OPM_THROW(std::logic_error, "FIP not yet implemented for MPI");
}
}
totals[6] = unit::convert::to( (p_pv_hydrocarbon_sum / pv_hydrocarbon_sum), unit::barsa);
return totals;
}
void outputFluidInPlace(const std::vector<double>& oip, const std::vector<double>& cip, const UnitSystem& units, const int reg)
{
std::ostringstream ss;
if (!reg) {
ss << " ===================================================\n"
<< " : Field Totals :\n";
} else {
ss << " ===================================================\n"
<< " : FIPNUM report region "
<< std::setw(2) << reg << " :\n";
}
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n"
<< std::fixed << std::setprecision(0)
<< " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n";
if (!reg) {
ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
<< " : Porv volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
}
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :\n"
<< std::fixed << std::setprecision(0)
<< " : PORV =" << std::setw(14) << cip[5] << " RB :\n";
if (!reg) {
ss << " : Pressure is weighted by hydrocarbon pore voulme :\n"
<< " : Pore volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";
}
ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n"
<< ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n"
<< ":Currently in place :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) << ":"
<< std::setw(13) << cip[0] << " :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) << ":\n"
<< ":------------------------:------------------------------------------:----------------:------------------------------------------:\n"
<< ":Originally in place :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) << ":"
<< std::setw(13) << oip[0] << " :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) << ":\n"
<< ":========================:==========================================:================:==========================================:\n";
OpmLog::note(ss.str());
}
const EclipseState& eclState() const
{ return *ebosSimulator_.gridManager().eclState(); }

View File

@ -162,10 +162,12 @@ namespace Opm
if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) {
std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl;
} else {
data::Solution combined_sol = simToSolution(state, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...)
combined_sol.insert(sol.begin(), sol.end()); // ... insert "extra" data (KR, VISC, ...)
eclWriter_->writeTimeStep(timer.reportStepNum(),
substep,
timer.simulationTimeElapsed(),
simToSolution( state, phaseUsage_ ),
combined_sol,
wellState.report(phaseUsage_));
}
}

View File

@ -538,6 +538,7 @@ namespace Opm
|| summaryConfig.hasKeyword(block_kw);
}
/**
* Returns the data as asked for in the summaryConfig
*/
@ -545,35 +546,34 @@ namespace Opm
void getSummaryData(data::Solution& output,
const Opm::PhaseUsage& phaseUsage,
const Model& physicalModel,
const SummaryConfig& summaryConfig)
{
const SummaryConfig& summaryConfig) {
const typename Model::FIPData& fip = physicalModel.getFIPData();
//Get shorthands for water, oil, gas
const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
#warning "TODO: fluid in place"
#if 0
#warning "TODO"
const int numFip = 0;
/**
* Now process all of the summary config files
*/
// Water in place
if (aqua_active && hasFRBKeyword(summaryConfig, "WIP")) {
#warning "TODO"
const std::vector<double> wip(numFip, 0.0);
output.insert("WIP",
Opm::UnitSystem::measure::volume,
wip,
fip.fip[Model::FIPData::FIP_AQUA],
data::TargetType::SUMMARY );
}
if (liquid_active) {
#warning "TODO"
const std::vector<double> oipl(numFip, 0.0);
const std::vector<double> oipg(numFip, 0.0);
const std::vector<double> oip(numFip, 0.0);
const std::vector<double>& oipl = fip.fip[Model::FIPData::FIP_LIQUID];
const int size = oipl.size();
const std::vector<double>& oipg = vapour_active ? fip.fip[Model::FIPData::FIP_VAPORIZED_OIL] : std::vector<double>(size,0.0);
std::vector<double> oip = oipl;
if (vapour_active) {
oip.insert(oip.end(), oipg.begin(), oipg.end());
}
//Oil in place (liquid phase only)
if (hasFRBKeyword(summaryConfig, "OIPL")) {
@ -582,6 +582,7 @@ namespace Opm
oipl,
data::TargetType::SUMMARY );
}
//Oil in place (gas phase only)
if (hasFRBKeyword(summaryConfig, "OIPG")) {
output.insert("OIPG",
@ -596,12 +597,19 @@ namespace Opm
oip,
data::TargetType::SUMMARY );
}
}
if (vapour_active) {
#warning "TODO"
const std::vector<double> gipl(numFip, 0.0);
const std::vector<double> gipg(numFip, 0.0);
const std::vector<double> gip(numFip, 0.0);
const std::vector<double>& gipg = fip.fip[Model::FIPData::FIP_VAPOUR];
const int size = gipg.size();
const std::vector<double>& gipl= liquid_active ? fip.fip[Model::FIPData::FIP_DISSOLVED_GAS] : std::vector<double>(size,0.0);
std::vector<double> gip = gipg;
if (liquid_active) {
gip.insert(gip.end(), gipl.begin(), gipl.end());
}
// Gas in place (gas phase only)
if (hasFRBKeyword(summaryConfig, "GIPG")) {
@ -627,24 +635,25 @@ namespace Opm
}
// Cell pore volume in reservoir conditions
if (hasFRBKeyword(summaryConfig, "RPV")) {
const std::vector<double> pv(numFip, 0.0);
output.insert("RPV",
Opm::UnitSystem::measure::volume,
pv,
fip.fip[Model::FIPData::FIP_PV],
data::TargetType::SUMMARY );
}
// Pressure averaged value (hydrocarbon pore volume weighted)
if (summaryConfig.hasKeyword("FPRH") || summaryConfig.hasKeyword("RPRH")) {
const std::vector<double> prh(numFip, 0.0);
output.insert("PRH",
Opm::UnitSystem::measure::pressure,
prh,
fip.fip[Model::FIPData::FIP_WEIGHTED_PRESSURE],
data::TargetType::SUMMARY );
}
#endif
}
}
template<class Model>
inline void
BlackoilOutputWriterEbos::