mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-01-16 05:21:56 -06:00
Merge pull request #802 from atgeirr/qilicun-fip-support
Output total fluid in place
This commit is contained in:
Atgeirr Flø Rasmussen
GitHub
commit
5f1bcb2a59
@ -240,7 +240,7 @@ try
|
||||
<< std::flush;
|
||||
|
||||
Opm::BlackoilOutputWriter
|
||||
outputWriter(cGrid, param, eclipseState, Opm::NNC(), pu,
|
||||
outputWriter(cGrid, param, eclipseState, pu,
|
||||
new_props->permeability() );
|
||||
|
||||
SimulatorReport fullReport;
|
||||
|
||||
@ -260,6 +260,14 @@ namespace Opm {
|
||||
WellModel& wellModel() { return well_model_; }
|
||||
const WellModel& wellModel() const { return well_model_; }
|
||||
|
||||
/// Compute fluid in place.
|
||||
/// \param[in] ReservoirState
|
||||
/// \param[in] FIPNUM for active cells not global cells.
|
||||
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
|
||||
std::vector<V>
|
||||
computeFluidInPlace(const ReservoirState& x,
|
||||
const std::vector<int>& fipnum);
|
||||
|
||||
protected:
|
||||
|
||||
// --------- Types and enums ---------
|
||||
|
||||
@ -2295,6 +2295,83 @@ namespace detail {
|
||||
|
||||
|
||||
|
||||
|
||||
template <class Grid, class WellModel, class Implementation>
|
||||
std::vector<V>
|
||||
BlackoilModelBase<Grid, WellModel, Implementation>::
|
||||
computeFluidInPlace(const ReservoirState& x,
|
||||
const std::vector<int>& fipnum)
|
||||
{
|
||||
using namespace Opm::AutoDiffGrid;
|
||||
const int nc = numCells(grid_);
|
||||
const int np = x.numPhases();
|
||||
|
||||
SolutionState state(np);
|
||||
const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
|
||||
state.pressure = ADB::constant(Eigen::Map<const V>(& x.pressure()[0], nc, 1));
|
||||
state.temperature = ADB::constant(Eigen::Map<const V>(& x.temperature()[0], nc, 1));
|
||||
state.saturation[Water] = ADB::constant(s.col(Water));
|
||||
state.saturation[Oil] = ADB::constant(s.col(Oil));
|
||||
state.saturation[Gas] = ADB::constant(s.col(Gas));
|
||||
state.rs = ADB::constant(Eigen::Map<const V>(& x.gasoilratio()[0], nc, 1));
|
||||
state.rv = ADB::constant(Eigen::Map<const V>(& x.rv()[0], nc, 1));
|
||||
state.canonical_phase_pressures = computePressures(state.pressure,
|
||||
state.saturation[Water],
|
||||
state.saturation[Oil],
|
||||
state.saturation[Gas]);
|
||||
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
|
||||
const std::vector<PhasePresence> cond = phaseCondition();
|
||||
|
||||
|
||||
|
||||
const ADB pv_mult = poroMult(state.pressure);
|
||||
const V& pv = geo_.poreVolume();
|
||||
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
|
||||
std::vector<V> fip(5, V::Zero(nc));
|
||||
for (int phase = 0; phase < maxnp; ++phase) {
|
||||
if (active_[ phase ]) {
|
||||
const int pos = pu.phase_pos[ phase ];
|
||||
const auto& b = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], state.temperature, state.rs, state.rv, cond);
|
||||
fip[phase] = ((pv_mult * b * state.saturation[pos] * pv).value());
|
||||
}
|
||||
}
|
||||
|
||||
if (active_[ Oil ] && active_[ Gas ]) {
|
||||
// Account for gas dissolved in oil and vaporized oil
|
||||
const int po = pu.phase_pos[Oil];
|
||||
const int pg = pu.phase_pos[Gas];
|
||||
fip[3] = state.rs.value() * fip[po];
|
||||
fip[4] = state.rv.value() * fip[pg];
|
||||
}
|
||||
|
||||
const int dims = *std::max_element(fipnum.begin(), fipnum.end());
|
||||
std::vector<V> values(dims, V::Zero(7));
|
||||
for (int i = 0; i < 5; ++i) {
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
if (fipnum[c] != 0) {
|
||||
values[fipnum[c]-1][i] += fip[i][c];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// compute PAV and PORV or every regions.
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
if (fipnum[c] != 0) {
|
||||
values[fipnum[c]-1][5] += pv[c];
|
||||
values[fipnum[c]-1][6] += pv[c] * state.pressure.value()[c];
|
||||
}
|
||||
}
|
||||
|
||||
for (auto& val : values) {
|
||||
val[6] = val[6] / val[5];
|
||||
}
|
||||
|
||||
return values;
|
||||
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
|
||||
|
||||
@ -247,7 +247,17 @@ namespace Opm {
|
||||
}
|
||||
|
||||
|
||||
|
||||
/// Compute fluid in place.
|
||||
/// \param[in] ReservoirState
|
||||
/// \param[in] WellState
|
||||
/// \param[in] FIPNUM for active cells not global cells.
|
||||
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
|
||||
std::vector<V>
|
||||
computeFluidInPlace(const ReservoirState& x,
|
||||
const std::vector<int>& fipnum) const
|
||||
{
|
||||
return transport_solver_.computeFluidInPlace(x, fipnum);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
@ -160,8 +160,6 @@ namespace Opm {
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/// Solve the Jacobian system Jx = r where J is the Jacobian and
|
||||
/// r is the residual.
|
||||
V solveJacobianSystem() const
|
||||
@ -275,6 +273,10 @@ namespace Opm {
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void assembleMassBalanceEq(const SolutionState& state)
|
||||
{
|
||||
// Compute b_p and the accumulation term b_p*s_p for each phase,
|
||||
|
||||
@ -663,7 +663,6 @@ namespace Opm
|
||||
output_writer_.reset(new BlackoilOutputWriter(grid_init_->grid(),
|
||||
param_,
|
||||
eclipse_state_,
|
||||
geoprops_->nonCartesianConnections(),
|
||||
Opm::phaseUsageFromDeck(deck_),
|
||||
fluidprops_->permeability()));
|
||||
}
|
||||
|
||||
@ -125,6 +125,14 @@ namespace Opm {
|
||||
/// Number of well iterations used in all calls to step().
|
||||
int wellIterationsLastStep() const;
|
||||
|
||||
/// Compute fluid in place.
|
||||
/// \param[in] ReservoirState
|
||||
/// \param[in] FIPNUM for active cells not global cells.
|
||||
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
|
||||
std::vector<V>
|
||||
computeFluidInPlace(const ReservoirState& x,
|
||||
const std::vector<int>& fipnum) const;
|
||||
|
||||
/// Reference to physical model.
|
||||
const PhysicalModel& model() const;
|
||||
|
||||
|
||||
@ -99,6 +99,15 @@ namespace Opm
|
||||
return wellIterationsLast_;
|
||||
}
|
||||
|
||||
template <class PhysicalModel>
|
||||
std::vector<V>
|
||||
NonlinearSolver<PhysicalModel>::computeFluidInPlace(const ReservoirState& x,
|
||||
const std::vector<int>& fipnum) const
|
||||
{
|
||||
return model_->computeFluidInPlace(x, fipnum);
|
||||
}
|
||||
|
||||
|
||||
template <class PhysicalModel>
|
||||
int
|
||||
NonlinearSolver<PhysicalModel>::
|
||||
|
||||
@ -158,6 +158,16 @@ namespace Opm
|
||||
const BlackoilState& x,
|
||||
WellState& xw);
|
||||
|
||||
void
|
||||
FIPUnitConvert(const UnitSystem& units,
|
||||
std::vector<V>& fip);
|
||||
|
||||
V
|
||||
FIPTotals(const std::vector<V>& fip, const std::vector<double>& press);
|
||||
|
||||
void
|
||||
outputFluidInPlace(const V& oip, const V& cip, const UnitSystem& units, const int reg);
|
||||
|
||||
void computeWellPotentials(const Wells* wells,
|
||||
const WellState& xw,
|
||||
std::vector<double>& well_potentials);
|
||||
|
||||
@ -135,6 +135,18 @@ namespace Opm
|
||||
std::vector<double> well_potentials;
|
||||
DynamicListEconLimited dynamic_list_econ_limited;
|
||||
|
||||
bool ooip_computed = false;
|
||||
std::vector<int> fipnum_global = eclipse_state_->get3DProperties().getIntGridProperty("FIPNUM").getData();
|
||||
//Get compressed cell fipnum.
|
||||
std::vector<int> fipnum(AutoDiffGrid::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[AutoDiffGrid::globalCell(grid_)[c]];
|
||||
}
|
||||
}
|
||||
std::vector<V> OOIP;
|
||||
// Main simulation loop.
|
||||
while (!timer.done()) {
|
||||
// Report timestep.
|
||||
@ -185,6 +197,13 @@ namespace Opm
|
||||
|
||||
auto solver = asImpl().createSolver(well_model);
|
||||
|
||||
// Compute orignal FIP;
|
||||
if (!ooip_computed) {
|
||||
OOIP = solver->computeFluidInPlace(state, fipnum);
|
||||
FIPUnitConvert(eclipse_state_->getUnits(), OOIP);
|
||||
ooip_computed = true;
|
||||
}
|
||||
|
||||
if( terminal_output_ )
|
||||
{
|
||||
std::ostringstream step_msg;
|
||||
@ -248,10 +267,21 @@ namespace Opm
|
||||
|
||||
// Report timing.
|
||||
const double st = solver_timer.secsSinceStart();
|
||||
// Compute current FIP.
|
||||
std::vector<V> COIP;
|
||||
COIP = solver->computeFluidInPlace(state, fipnum);
|
||||
FIPUnitConvert(eclipse_state_->getUnits(), COIP);
|
||||
V OOIP_totals = FIPTotals(OOIP, state.pressure());
|
||||
V COIP_totals = FIPTotals(COIP, state.pressure());
|
||||
outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0);
|
||||
for (size_t reg = 0; reg < OOIP.size(); ++reg) {
|
||||
outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1);
|
||||
}
|
||||
|
||||
|
||||
// accumulate total time
|
||||
stime += st;
|
||||
|
||||
|
||||
if ( terminal_output_ )
|
||||
{
|
||||
std::string msg;
|
||||
@ -626,8 +656,89 @@ namespace Opm
|
||||
}
|
||||
|
||||
|
||||
template <class Implementation>
|
||||
void
|
||||
SimulatorBase<Implementation>::FIPUnitConvert(const UnitSystem& units,
|
||||
std::vector<V>& 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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <class Implementation>
|
||||
V
|
||||
SimulatorBase<Implementation>::FIPTotals(const std::vector<V>& fip, const std::vector<double>& press)
|
||||
{
|
||||
V totals(V::Zero(7));
|
||||
for (int i = 0; i < 5; ++i) {
|
||||
for (size_t reg = 0; reg < fip.size(); ++reg) {
|
||||
totals[i] += fip[reg][i];
|
||||
}
|
||||
}
|
||||
const V p = Eigen::Map<const V>(& press[0], press.size());
|
||||
totals[5] = geo_.poreVolume().sum();
|
||||
totals[6] = unit::convert::to((p * geo_.poreVolume()).sum() / totals[5], unit::barsa);
|
||||
|
||||
return totals;
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class Implementation>
|
||||
void
|
||||
SimulatorBase<Implementation>::outputFluidInPlace(const V& oip, const V& 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"
|
||||
<< " : Pressure is weighted by hydrocarbon pore voulme:\n"
|
||||
<< " : Porv volume are taken at reference conditions :\n"
|
||||
<< " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
|
||||
}
|
||||
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
|
||||
ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :"
|
||||
<< std::fixed << std::setprecision(0)
|
||||
<< " : PORV =" << std::setprecision(14) << cip[5] << " STB :"
|
||||
<< " : Pressure is weighted by hydrocarbon pore voulme :"
|
||||
<< " : Pore volume are taken at reference conditions :"
|
||||
<< " :--------------- 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());
|
||||
}
|
||||
|
||||
|
||||
template <class Implementation>
|
||||
void
|
||||
|
||||
@ -212,7 +212,6 @@ namespace Opm
|
||||
BlackoilOutputWriter(const Grid& grid,
|
||||
const parameter::ParameterGroup& param,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
const NNC&,
|
||||
const Opm::PhaseUsage &phaseUsage,
|
||||
const double* permeability );
|
||||
|
||||
@ -285,7 +284,6 @@ namespace Opm
|
||||
BlackoilOutputWriter(const Grid& grid,
|
||||
const parameter::ParameterGroup& param,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
const NNC& nnc,
|
||||
const Opm::PhaseUsage &phaseUsage,
|
||||
const double* permeability )
|
||||
: output_( param.getDefault("output", true) ),
|
||||
|
||||
@ -520,7 +520,7 @@ namespace {
|
||||
}
|
||||
rq_[0].accum[aix] = pv_mult * rq_[0].b * sat[0];
|
||||
rq_[1].accum[aix] = pv_mult * rq_[1].b * sat[1];
|
||||
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
|
||||
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
|
||||
const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
|
||||
const double rho_rock = polymer_props_ad_.rockDensity();
|
||||
const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], grid_.number_of_cells, 1);
|
||||
@ -532,6 +532,56 @@ namespace {
|
||||
|
||||
|
||||
|
||||
std::vector<V>
|
||||
FullyImplicitCompressiblePolymerSolver::computeFluidInPlace(const PolymerBlackoilState& x,
|
||||
const std::vector<int>& fipnum)
|
||||
{
|
||||
const int np = x.numPhases();
|
||||
const int nc = grid_.number_of_cells;
|
||||
|
||||
SolutionState state(np);
|
||||
state.pressure = ADB::constant(Eigen::Map<const V>(& x.pressure()[0], nc, 1));
|
||||
state.temperature = ADB::constant(Eigen::Map<const V>(& x.temperature()[0], nc, 1));
|
||||
const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
state.saturation[phase] = ADB::constant(s.col(phase));
|
||||
}
|
||||
|
||||
const ADB& press = state.pressure;
|
||||
const ADB& temp = state.temperature;
|
||||
const std::vector<ADB>& sat = state.saturation;
|
||||
|
||||
const std::vector<PhasePresence> cond = phaseCondition();
|
||||
std::vector<ADB> pressure = computePressures(state);
|
||||
|
||||
const ADB pv_mult = poroMult(press);
|
||||
const V& pv = geo_.poreVolume();
|
||||
std::vector<V> fip(5, V::Zero(nc));
|
||||
for (int phase = 0; phase < 2; ++phase) {
|
||||
const ADB& b = fluidReciprocFVF(phase, pressure[phase], temp, cond, cells_);
|
||||
fip[phase] = (pv_mult * b * sat[phase] * pv).value();
|
||||
}
|
||||
|
||||
|
||||
const int dims = *std::max_element(fipnum.begin(), fipnum.end());
|
||||
std::vector<V> values(dims, V::Zero(5));
|
||||
|
||||
for (int d = 0; d < dims; ++d) {
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
for (int i = 0; i < 5; ++i) {
|
||||
if (fipnum[c] == d) {
|
||||
values[d][i] += fip[c][i];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return values;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void
|
||||
FullyImplicitCompressiblePolymerSolver::
|
||||
|
||||
@ -54,6 +54,13 @@ namespace Opm {
|
||||
class FullyImplicitCompressiblePolymerSolver
|
||||
{
|
||||
public:
|
||||
typedef AutoDiffBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
typedef ADB::M M;
|
||||
typedef Eigen::Array<double,
|
||||
Eigen::Dynamic,
|
||||
Eigen::Dynamic,
|
||||
Eigen::RowMajor> DataBlock;
|
||||
/// Construct a solver. It will retain references to the
|
||||
/// arguments of this functions, and they are expected to
|
||||
/// remain in scope for the lifetime of the solver.
|
||||
@ -102,14 +109,17 @@ namespace Opm {
|
||||
double relativeChange(const PolymerBlackoilState& previous,
|
||||
const PolymerBlackoilState& current ) const;
|
||||
|
||||
/// Compute fluid in place.
|
||||
/// \param[in] ReservoirState
|
||||
/// \param[in] WellState
|
||||
/// \param[in] FIPNUM for active cells not global cells.
|
||||
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
|
||||
std::vector<V>
|
||||
computeFluidInPlace(const PolymerBlackoilState& x,
|
||||
const std::vector<int>& fipnum);
|
||||
|
||||
private:
|
||||
typedef AutoDiffBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
typedef ADB::M M;
|
||||
typedef Eigen::Array<double,
|
||||
Eigen::Dynamic,
|
||||
Eigen::Dynamic,
|
||||
Eigen::RowMajor> DataBlock;
|
||||
|
||||
|
||||
struct ReservoirResidualQuant {
|
||||
ReservoirResidualQuant();
|
||||
|
||||
Loading…
Reference in New Issue
Block a user