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

port the RateConverter to use the FluidSystem instead of the old fluid property API

Browse Source 
this makes the RateConverter stuff independent of Eigen and it
simplifies some things because the the old PVT API is designed as a
"bulk-with-derivatives" API while the rate converter code used it in
"single shot" mode without derivatives.
This commit is contained in:
Andreas Lauser 2016-12-29 16:56:31 +01:00
parent 66decb4bda
commit 5fd83985a9
7 changed files with 83 additions and 113 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
opm/autodiff/BlackoilModelBase.hpp
View File

@ -117,7 +117,7 @@ namespace Opm {
// For the conversion between the surface volume rate and resrevoir voidage rate // For the conversion between the surface volume rate and resrevoir voidage rate
using RateConverterType = RateConverter:: using RateConverterType = RateConverter::
SurfaceToReservoirVoidage<BlackoilPropsAdFromDeck, std::vector<int> >; SurfaceToReservoirVoidage<BlackoilPropsAdFromDeck::FluidSystem, std::vector<int> >;
// --------- Public methods --------- // --------- Public methods ---------

2
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -144,7 +144,7 @@ typedef Eigen::Array<double,
// the field will be calculated. // the field will be calculated.
// TODO: more delicate implementation will be required if we want to handle different // TODO: more delicate implementation will be required if we want to handle different
// FIP regions specified from the well specifications. // FIP regions specified from the well specifications.
, rate_converter_(fluid_, std::vector<int>(AutoDiffGrid::numCells(grid_),0)) , rate_converter_(fluid_.phaseUsage(), fluid_.cellPvtRegionIndex(), AutoDiffGrid::numCells(grid_), std::vector<int>(AutoDiffGrid::numCells(grid_),0))
{ {
if (active_[Water]) { if (active_[Water]) {
material_name_.push_back("Water"); material_name_.push_back("Water");

143
opm/autodiff/RateConverter.hpp
View File

@ -28,8 +28,6 @@
#include <opm/core/utility/RegionMapping.hpp> #include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp> #include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <Eigen/Core>
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
#include <memory> #include <memory>
@ -385,9 +383,7 @@ namespace Opm {
* The conversion uses fluid properties evaluated at average * The conversion uses fluid properties evaluated at average
* hydrocarbon pressure in regions or field. * hydrocarbon pressure in regions or field.
* *
* \tparam Property Fluid property object. Expected to * \tparam FluidSystem Fluid system class. Expected to be a BlackOilFluidSystem
* feature the formation volume factor functions of the
* BlackoilPropsAdInterface.
* *
* \tparam Region Type of a forward region mapping. Expected * \tparam Region Type of a forward region mapping. Expected
* to provide indexed access through \code operator[]() * to provide indexed access through \code operator[]()
@ -395,24 +391,29 @@ namespace Opm {
* size_type, and \c const_iterator. Typically \code * size_type, and \c const_iterator. Typically \code
* std::vector<int> \endcode. * std::vector<int> \endcode.
*/ */
template <class Property, class Region> template <class FluidSystem, class Region>
class SurfaceToReservoirVoidage { class SurfaceToReservoirVoidage {
public: public:
/** /**
* Constructor. * Constructor.
* *
* \param[in] props Fluid property object.
*
* \param[in] region Forward region mapping. Often * \param[in] region Forward region mapping. Often
* corresponds to the "FIPNUM" mapping of an ECLIPSE input * corresponds to the "FIPNUM" mapping of an ECLIPSE input
* deck. * deck.
*/ */
SurfaceToReservoirVoidage(const Property& props, SurfaceToReservoirVoidage(const PhaseUsage& phaseUsage,
const int* cellPvtRegionIdx,
const int numCells,
const Region& region) const Region& region)
: props_(props) : phaseUsage_(phaseUsage)
, rmap_ (region) , rmap_ (region)
, attr_ (rmap_, Attributes(props_.numPhases())) , attr_ (rmap_, Attributes(phaseUsage_.num_phases))
{} {
cellPvtIdx_.resize(numCells, 0);
if (cellPvtRegionIdx) {
std::copy_n(cellPvtRegionIdx, numCells, cellPvtIdx_.begin());
}
}
/** /**
* Compute average hydrocarbon pressure and maximum * Compute average hydrocarbon pressure and maximum
@ -484,59 +485,58 @@ namespace Opm {
void void
calcCoeff(const Input& in, const RegionId r, Coeff& coeff) const calcCoeff(const Input& in, const RegionId r, Coeff& coeff) const
{ {
using V = typename Property::V; const auto& pu = phaseUsage_;
const auto& pu = props_.phaseUsage();
const auto& ra = attr_.attributes(r); const auto& ra = attr_.attributes(r);
const auto p = this->constant(ra.pressure); const double p = ra.pressure;
const auto T = this->constant(ra.temperature); const double T = ra.temperature;
const auto c = this->getRegCell(r); const int cellIdx = attr_.cell(r);
const int pvtRegionIdx = cellPvtIdx_[cellIdx];
const int iw = Details::PhasePos::water(pu); const int iw = Details::PhasePos::water(pu);
const int io = Details::PhasePos::oil (pu); const int io = Details::PhasePos::oil (pu);
const int ig = Details::PhasePos::gas (pu); const int ig = Details::PhasePos::gas (pu);
std::fill(& coeff[0], & coeff[0] + props_.numPhases(), 0.0); std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
if (Details::PhaseUsed::water(pu)) { if (Details::PhaseUsed::water(pu)) {
// q[w]_r = q[w]_s / bw // q[w]_r = q[w]_s / bw
const V bw = props_.bWat(p, T, c).value(); const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p);
coeff[iw] = 1.0 / bw(0); coeff[iw] = 1.0 / bw;
} }
const Miscibility& m = calcMiscibility(in, r); const Miscibility& m = calcMiscibility(in, r);
// Determinant of 'R' matrix // Determinant of 'R' matrix
const double detR = 1.0 - (m.rs(0) * m.rv(0)); const double detR = 1.0 - (m.rs * m.rv);
if (Details::PhaseUsed::oil(pu)) { if (Details::PhaseUsed::oil(pu)) {
// q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s) // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
const auto rs = this->constant(m.rs); const double Rs = m.rs;
const V bo = props_.bOil(p, T, rs, m.cond, c).value(); const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
const double den = bo(0) * detR; const double den = bo * detR;
coeff[io] += 1.0 / den; coeff[io] += 1.0 / den;
if (Details::PhaseUsed::gas(pu)) { if (Details::PhaseUsed::gas(pu)) {
coeff[ig] -= m.rv(0) / den; coeff[ig] -= m.rv / den;
} }
} }
if (Details::PhaseUsed::gas(pu)) { if (Details::PhaseUsed::gas(pu)) {
// q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s) // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
const auto rv = this->constant(m.rv); const double Rv = m.rv;
const V bg = props_.bGas(p, T, rv, m.cond, c).value(); const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
const double den = bg(0) * detR; const double den = bg * detR;
coeff[ig] += 1.0 / den; coeff[ig] += 1.0 / den;
if (Details::PhaseUsed::oil(pu)) { if (Details::PhaseUsed::oil(pu)) {
coeff[io] -= m.rs(0) / den; coeff[io] -= m.rs / den;
} }
} }
} }
@ -545,7 +545,8 @@ namespace Opm {
/** /**
* Fluid property object. * Fluid property object.
*/ */
const Property& props_; const PhaseUsage phaseUsage_;
std::vector<int> cellPvtIdx_;
/** /**
* "Fluid-in-place" region mapping (forward and reverse). * "Fluid-in-place" region mapping (forward and reverse).
@ -559,12 +560,12 @@ namespace Opm {
Attributes(const int np) Attributes(const int np)
: pressure (0.0) : pressure (0.0)
, temperature(0.0) , temperature(0.0)
, Rmax (Eigen::ArrayXd::Zero(np, 1)) , Rmax(np, 0.0)
{} {}
double pressure; double pressure;
double temperature; double temperature;
Eigen::ArrayXd Rmax; std::vector<double> Rmax;
}; };
Details::RegionAttributes<RegionId, Attributes> attr_; Details::RegionAttributes<RegionId, Attributes> attr_;
@ -580,8 +581,8 @@ namespace Opm {
, rv (1) , rv (1)
, cond(1) , cond(1)
{ {
rs << 0.0; rs = 0.0;
rv << 0.0; rv = 0.0;
} }
/** /**
@ -591,7 +592,7 @@ namespace Opm {
* Limited by "RSmax" at average hydrocarbon pressure * Limited by "RSmax" at average hydrocarbon pressure
* in region. * in region.
*/ */
typename Property::V rs; double rs;
/** /**
* Evaporated oil-gas ratio at particular component oil * Evaporated oil-gas ratio at particular component oil
@ -600,7 +601,7 @@ namespace Opm {
* Limited by "RVmax" at average hydrocarbon pressure * Limited by "RVmax" at average hydrocarbon pressure
* in region. * in region.
*/ */
typename Property::V rv; double rv;
/** /**
* Fluid condition in representative region cell. * Fluid condition in representative region cell.
@ -674,14 +675,13 @@ namespace Opm {
void void
calcRmax() calcRmax()
{ {
const PhaseUsage& pu = props_.phaseUsage(); const PhaseUsage& pu = phaseUsage_;
if (Details::PhaseUsed::oil(pu) && if (Details::PhaseUsed::oil(pu) &&
Details::PhaseUsed::gas(pu)) Details::PhaseUsed::gas(pu))
{ {
const Eigen::ArrayXd::Index const int io = Details::PhasePos::oil(pu);
io = Details::PhasePos::oil(pu), const int ig = Details::PhasePos::gas(pu);
ig = Details::PhasePos::gas(pu);
// Note: Intentionally does not take capillary // Note: Intentionally does not take capillary
// pressure into account. This facility uses the // pressure into account. This facility uses the
@ -691,13 +691,14 @@ namespace Opm {
for (const auto& reg : rmap_.activeRegions()) { for (const auto& reg : rmap_.activeRegions()) {
auto& ra = attr_.attributes(reg); auto& ra = attr_.attributes(reg);
const auto c = this->getRegCell(reg); const double T = ra.temperature;
const auto p = this->constant(ra.pressure); const double p = ra.pressure;
const auto T = this->constant(ra.temperature); const int cellIdx = attr_.cell(reg);
const int pvtRegionIdx = cellPvtIdx_[cellIdx];
auto& Rmax = ra.Rmax; std::vector<double>& Rmax = ra.Rmax;
Rmax.row(io) = props_.rsSat(p, T, c).value(); Rmax[io] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx, T, p);
Rmax.row(ig) = props_.rvSat(p, T, c).value(); Rmax[ig] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx, T, p);
} }
} }
} }
@ -723,7 +724,7 @@ namespace Opm {
Miscibility Miscibility
calcMiscibility(const Input& in, const RegionId r) const calcMiscibility(const Input& in, const RegionId r) const
{ {
const auto& pu = props_.phaseUsage(); const auto& pu = phaseUsage_;
const auto& attr = attr_.attributes(r); const auto& attr = attr_.attributes(r);
const int io = Details::PhasePos::oil(pu); const int io = Details::PhasePos::oil(pu);
@ -740,7 +741,7 @@ namespace Opm {
cond.setFreeOil(); cond.setFreeOil();
if (Details::PhaseUsed::gas(pu)) { if (Details::PhaseUsed::gas(pu)) {
const double rsmax = attr.Rmax(io); const double rsmax = attr.Rmax[io];
const double rs = const double rs =
(0.0 < std::abs(in[io])) (0.0 < std::abs(in[io]))
? in[ig] / in[io] ? in[ig] / in[io]
@ -750,7 +751,7 @@ namespace Opm {
cond.setFreeGas(); cond.setFreeGas();
} }
m.rs(0) = std::min(rs, rsmax); m.rs = std::min(rs, rsmax);
} }
} }
@ -761,56 +762,18 @@ namespace Opm {
} }
if (Details::PhaseUsed::oil(pu)) { if (Details::PhaseUsed::oil(pu)) {
const double rvmax = attr.Rmax(ig); const double rvmax = attr.Rmax[ig];
const double rv = const double rv =
(0.0 < std::abs(in[ig])) (0.0 < std::abs(in[ig]))
? (in[io] / in[ig]) ? (in[io] / in[ig])
: (0.0 < std::abs(in[io])) ? rvmax : 0.0; : (0.0 < std::abs(in[io])) ? rvmax : 0.0;
m.rv(0) = std::min(rv, rvmax); m.rv = std::min(rv, rvmax);
} }
} }
return m; return m;
} }
/**
* Conversion from \code Property::V \endcode to (constant)
* \code Property::ADB \endcode (zero derivatives).
*/
typename Property::ADB
constant(const typename Property::V& x) const
{
return Property::ADB::constant(x);
}
/**
* Conversion from \c double to (constant) \code Property::ADB
* \endcode (zero derivatives).
*/
typename Property::ADB
constant(const double x) const
{
typename Property::V y(1);
y << x;
return this->constant(y);
}
/**
* Retrieve representative cell of region
*
* \param[in] r Particular region.
*
* \return Representative cell of region \c r.
*/
typename Property::Cells
getRegCell(const RegionId r) const
{
typename Property::Cells c(1, this->attr_.cell(r));
return c;
}
}; };
} // namespace RateConverter } // namespace RateConverter
} // namespace Opm } // namespace Opm

2
opm/autodiff/SimulatorBase.hpp
View File

@ -189,7 +189,7 @@ namespace Opm
// Data. // Data.
typedef RateConverter:: typedef RateConverter::
SurfaceToReservoirVoidage< BlackoilPropsAdFromDeck, SurfaceToReservoirVoidage< BlackoilPropsAdFromDeck::FluidSystem,
std::vector<int> > RateConverterType; std::vector<int> > RateConverterType;
typedef typename Traits::Model Model; typedef typename Traits::Model Model;
typedef typename Model::ModelParameters ModelParameters; typedef typename Model::ModelParameters ModelParameters;

2
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -60,7 +60,7 @@ namespace Opm
terminal_output_(param.getDefault("output_terminal", true)), terminal_output_(param.getDefault("output_terminal", true)),
eclipse_state_(eclipse_state), eclipse_state_(eclipse_state),
output_writer_(output_writer), output_writer_(output_writer),
rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0)), rateConverter_(props_.phaseUsage(), props.cellPvtRegionIndex(), AutoDiffGrid::numCells(grid_), std::vector<int>(AutoDiffGrid::numCells(grid_), 0)),
threshold_pressures_by_face_(threshold_pressures_by_face), threshold_pressures_by_face_(threshold_pressures_by_face),
is_parallel_run_( false ), is_parallel_run_( false ),
defunct_well_names_(defunct_well_names) defunct_well_names_(defunct_well_names)

33
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -117,15 +117,11 @@ public:
defunct_well_names_( defunct_well_names ), defunct_well_names_( defunct_well_names ),
is_parallel_run_( false ) is_parallel_run_( false )
{ {
DUNE_UNUSED_PARAMETER(eclState); extractLegacyCellPvtRegionIndex_();
// Misc init. rateConverter_.reset(new RateConverterType(props.phaseUsage(),
const int num_cells = AutoDiffGrid::numCells(grid()); legacyCellPvtRegionIdx_.data(),
allcells_.resize(num_cells); AutoDiffGrid::numCells(grid()),
for (int cell = 0; cell < num_cells; ++cell) { std::vector<int>(AutoDiffGrid::numCells(grid()), 0)));
allcells_[cell] = cell;
}
rateConverter_.reset(new RateConverterType(props_, std::vector<int>(AutoDiffGrid::numCells(grid()), 0)));
#if HAVE_MPI #if HAVE_MPI
if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
@ -708,15 +704,24 @@ protected:
const EclipseState& eclState() const const EclipseState& eclState() const
{ return ebosSimulator_.gridManager().eclState(); } { return ebosSimulator_.gridManager().eclState(); }
EclipseState& eclState() void extractLegacyCellPvtRegionIndex_()
{ return ebosSimulator_.gridManager().eclState(); } {
const auto& grid = ebosSimulator_.gridManager().grid();
const auto& eclProblem = ebosSimulator_.problem();
const unsigned numCells = grid.size(/*codim=*/0);
legacyCellPvtRegionIdx_.resize(numCells);
for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
legacyCellPvtRegionIdx_[cellIdx] =
eclProblem.pvtRegionIndex(cellIdx);
}
}
// Data. // Data.
Simulator& ebosSimulator_; Simulator& ebosSimulator_;
typedef RateConverter:: std::vector<int> legacyCellPvtRegionIdx_;
SurfaceToReservoirVoidage< BlackoilPropsAdFromDeck, typedef RateConverter::SurfaceToReservoirVoidage<FluidSystem, std::vector<int> > RateConverterType;
std::vector<int> > RateConverterType;
typedef typename Solver::SolverParameters SolverParameters; typedef typename Solver::SolverParameters SolverParameters;
const parameter::ParameterGroup param_; const parameter::ParameterGroup param_;

12
tests/test_rateconverter.cpp
View File

@ -87,10 +87,11 @@ BOOST_FIXTURE_TEST_CASE(Construction, TestFixture<SetupSimple>)
typedef std::vector<int> Region; typedef std::vector<int> Region;
typedef Opm::BlackoilPropsAdFromDeck Props; typedef Opm::BlackoilPropsAdFromDeck Props;
typedef Opm::RateConverter:: typedef Opm::RateConverter::
SurfaceToReservoirVoidage<Props, Region> RCvrt; SurfaceToReservoirVoidage<Props::FluidSystem, Region> RCvrt;
Region reg{ 0 }; Region reg{ 0 };
RCvrt cvrt(ad_props, reg); int numCells = Opm::UgGridHelpers::numCells(*grid.c_grid());
RCvrt cvrt(ad_props.phaseUsage(), ad_props.cellPvtRegionIndex(), numCells, reg);
} }
@ -100,12 +101,13 @@ BOOST_FIXTURE_TEST_CASE(ThreePhase, TestFixture<SetupSimple>)
typedef std::vector<int> Region; typedef std::vector<int> Region;
typedef Opm::BlackoilPropsAdFromDeck Props; typedef Opm::BlackoilPropsAdFromDeck Props;
typedef Opm::RateConverter:: typedef Opm::RateConverter::
SurfaceToReservoirVoidage<Props, Region> RCvrt; SurfaceToReservoirVoidage<Props::FluidSystem, Region> RCvrt;
Region reg{ 0 }; Region reg{ 0 };
RCvrt cvrt(ad_props, reg); int numCells = Opm::UgGridHelpers::numCells(*grid.c_grid());
RCvrt cvrt(ad_props.phaseUsage(), ad_props.cellPvtRegionIndex(), numCells, reg);
Opm::BlackoilState x( Opm::UgGridHelpers::numCells( *grid.c_grid()) , Opm::UgGridHelpers::numFaces( *grid.c_grid()) , 3); Opm::BlackoilState x(numCells, Opm::UgGridHelpers::numFaces( *grid.c_grid()) , 3);
cvrt.defineState(x); cvrt.defineState(x);