OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-10 02:25:33 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #188 from OPM/support-resv

Browse Source 
Support RESV controls
This commit is contained in:
Bård Skaflestad 2014-08-25 13:54:30 +02:00
commit 077b01b767
8 changed files with 1067 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -46,6 +46,7 @@ list (APPEND MAIN_SOURCE_FILES
list (APPEND TEST_SOURCE_FILES
tests/test_block.cpp
tests/test_boprops_ad.cpp
tests/test_rateconverter.cpp
tests/test_span.cpp
tests/test_syntax.cpp
tests/test_scalar_mult.cpp
@ -104,6 +105,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/NewtonIterationBlackoilInterface.hpp
opm/autodiff/NewtonIterationBlackoilSimple.hpp
opm/autodiff/LinearisedBlackoilResidual.hpp
opm/autodiff/RateConverter.hpp
opm/autodiff/SimulatorCompressibleAd.hpp
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

3
examples/sim_fibo_ad_cp.cpp
View File

@ -217,9 +217,6 @@ try
param.writeParam(output_dir + "/simulation.param");
}
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
SimulatorTimer simtimer;

2
opm/autodiff/FullyImplicitBlackoilSolver.hpp
View File

@ -140,7 +140,7 @@ namespace Opm {
const NewtonIterationBlackoilInterface& linsolver_;
// For each canonical phase -> true if active
const std::vector<bool> active_;
// Size = # active faces. Maps active -> canonical phase indices.
// Size = # active phases. Maps active -> canonical phase indices.
const std::vector<int> canph_;
const std::vector<int> cells_; // All grid cells
HelperOps ops_;

139
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
View File

@ -417,14 +417,18 @@ namespace {
if (active_[ Gas ]){
for (int c = 0; c < nc ; c++ ) {
if ( primalVariable_[c] == PrimalVariables::RS ) {
switch (primalVariable_[c]) {
case PrimalVariables::RS:
isRs[c] = 1;
}
else if ( primalVariable_[c] == PrimalVariables::RV ) {
break;
case PrimalVariables::RV:
isRv[c] = 1;
}
else {
break;
default:
isSg[c] = 1;
break;
}
}
@ -955,37 +959,50 @@ namespace {
const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
const double target = well_controls_iget_target(wc, ctrl_index);
const double* distr = well_controls_iget_distr(wc, ctrl_index);
bool broken = false;
switch (well_type) {
case INJECTOR:
{
switch (ctrl_type) {
case BHP:
return bhp.value()[well] > target;
broken = bhp.value()[well] > target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
return rateToCompare(well_phase_flow_rate, well, num_phases, distr) > target;
case RESERVOIR_RATE:
// Intentional fall-through.
default:
OPM_THROW(std::logic_error, "Can only handle BHP and SURFACE_RATE controls.");
broken = rateToCompare(well_phase_flow_rate,
well, num_phases, distr) > target;
break;
}
break;
}
break;
case PRODUCER:
{
switch (ctrl_type) {
case BHP:
return bhp.value()[well] < target;
broken = bhp.value()[well] < target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
// Note that the rates compared below are negative,
// so breaking the constraints means: too high flow rate
// (as for injection).
return rateToCompare(well_phase_flow_rate, well, num_phases, distr) < target;
case RESERVOIR_RATE:
// Intentional fall-through.
default:
OPM_THROW(std::logic_error, "Can only handle BHP and SURFACE_RATE controls.");
broken = rateToCompare(well_phase_flow_rate,
well, num_phases, distr) < target;
break;
}
break;
}
break;
default:
OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells.");
}
return broken;
}
} // anonymous namespace
@ -1022,14 +1039,6 @@ namespace {
// inequality constraint, and therefore skipped.
continue;
}
if (well_controls_iget_type(wc, ctrl_index) == RESERVOIR_RATE) {
// We cannot handle this yet.
#ifdef OPM_VERBOSE
std::cout << "Warning: a RESERVOIR_RATE well control exists for well "
<< wells_.name[w] << ", but will never be checked." << std::endl;
#endif
continue;
}
if (constraintBroken(bhp, well_phase_flow_rate, w, np, wells_.type[w], wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
@ -1051,6 +1060,18 @@ namespace {
xw.bhp()[w] = target;
bhp_changed = true;
break;
case RESERVOIR_RATE:
// No direct change to any observable quantity at
// surface condition. In this case, use existing
// flow rates as initial conditions as reservoir
// rate acts only in aggregate.
//
// Just record the fact that we need to recompute
// the 'well_phase_flow_rate'.
rates_changed = true;
break;
case SURFACE_RATE:
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0) {
@ -1059,9 +1080,6 @@ namespace {
}
rates_changed = true;
break;
default:
OPM_THROW(std::logic_error, "Programmer error: should not have switched "
"to anything but BHP or SURFACE_RATE.");
}
}
}
@ -1089,13 +1107,11 @@ namespace {
const WellStateFullyImplicitBlackoil& xw,
const V& aliveWells)
{
// Handling BHP and SURFACE_RATE wells.
const int np = wells_.number_of_phases;
const int nw = wells_.number_of_wells;
V bhp_targets(nw);
V rate_targets(nw);
V bhp_targets = V::Zero(nw);
V rate_targets = V::Zero(nw);
M rate_distr(nw, np*nw);
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells_.ctrls[w];
@ -1103,20 +1119,33 @@ namespace {
// the current control set in the Wells struct, which
// is instead treated as a default.
const int current = xw.currentControls()[w];
if (well_controls_iget_type(wc, current) == BHP) {
bhp_targets[w] = well_controls_iget_target(wc, current);
rate_targets[w] = -1e100;
} else if (well_controls_iget_type(wc, current) == SURFACE_RATE) {
bhp_targets[w] = -1e100;
rate_targets[w] = well_controls_iget_target(wc, current);
{
const double * distr = well_controls_iget_distr(wc, current);
for (int phase = 0; phase < np; ++phase) {
rate_distr.insert(w, phase*nw + w) = distr[phase];
}
switch (well_controls_iget_type(wc, current)) {
case BHP:
{
bhp_targets (w) = well_controls_iget_target(wc, current);
rate_targets(w) = -1e100;
}
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
{
// RESERVOIR and SURFACE rates look the same, from a
// high-level point of view, in the system of
// simultaneous linear equations.
const double* const distr =
well_controls_iget_distr(wc, current);
for (int p = 0; p < np; ++p) {
rate_distr.insert(w, p*nw + w) = distr[p];
}
} else {
OPM_THROW(std::runtime_error, "Can only handle BHP and SURFACE_RATE type controls.");
bhp_targets (w) = -1.0e100;
rate_targets(w) = well_controls_iget_target(wc, current);
}
break;
}
}
const ADB bhp_residual = state.bhp - bhp_targets;
@ -1181,15 +1210,19 @@ namespace {
V isRv = V::Zero(nc,1);
V isSg = V::Zero(nc,1);
if (active_[Gas]) {
for (int c = 0; c < nc ; c++ ) {
if ( primalVariable_[c] == PrimalVariables::RS ) {
for (int c = 0; c < nc; ++c) {
switch (primalVariable_[c]) {
case PrimalVariables::RS:
isRs[c] = 1;
}
else if ( primalVariable_[c] == PrimalVariables::RV ) {
break;
case PrimalVariables::RV:
isRv[c] = 1;
}
else {
break;
default:
isSg[c] = 1;
break;
}
}
}

612
opm/autodiff/RateConverter.hpp Normal file
View File

@ -0,0 +1,612 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_RATECONVERTER_HPP_HEADER_INCLUDED
#define OPM_RATECONVERTER_HPP_HEADER_INCLUDED
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <Eigen/Core>
#include <algorithm>
#include <cmath>
#include <vector>
/**
* \file
* Facility for converting component rates at surface conditions to
* phase (voidage) rates at reservoir conditions.
*
* This uses the average hydrocarbon pressure to define fluid
* properties. The facility is intended to support Reservoir Voidage
* rates only ('RESV').
*/
namespace Opm {
namespace RateConverter {
/**
* Convenience tools for implementing the rate conversion
* facility.
*/
namespace Details {
/**
* Count number of cells in all regions.
*
* This value is needed to compute the average (arithmetic
* mean) hydrocarbon pressure in each region.
*
* \tparam RMap Region mapping. Typically an instance of
* class Opm::RegionMapping<> from module opm-core.
*
* \param[in] m Specific region mapping.
*
* \return Array containing number of cells in each region
* defined by the region mapping.
*/
template <class RMap>
Eigen::ArrayXd
countCells(const RMap& m)
{
// Note: Floating point type (double) to represent
// cell counts is intentional. The count will be
// subsequently used to compute average (pressure)
// values only, and that operation is safer if we
// guarantee a floating point type here.
Eigen::ArrayXd n(m.numRegions());
for (typename RMap::RegionId
r = 0, nr = m.numRegions(); r < nr; ++r)
{
n(r) = double(m.cells(r).size());
}
return n;
}
/**
* Extract representative cell in each region.
*
* These are the cells for which fluid properties will be
* computed.
*
* \tparam Cells Type of cell container. Must be
* commensurable with the properties object's input
* requirements and support a single (integer) argument
* constructor that specifies the number of regions.
* Typically \code std::vector<int> \endcode.
*
* \tparam RMap Region mapping. Typically an instance of
* class Opm::RegionMapping<> from module opm-core.
*
* \param[in] m Specific region mapping.
*
* \return Array of representative cells, one cell in each
* region defined by @c m.
*/
template <class Cells, class RMap>
Cells
representative(const RMap& m)
{
Cells c(m.numRegions());
for (typename RMap::RegionId
r = 0, nr = m.numRegions(); r < nr; ++r)
{
c[r] = *m.cells(r).begin();
}
return c;
}
/**
* Convenience functions for querying presence/absence of
* active phases.
*/
namespace PhaseUsed {
/**
* Active water predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not water is an active phase.
*/
inline bool
water(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Aqua ] != 0;
}
/**
* Active oil predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not oil is an active phase.
*/
inline bool
oil(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Liquid ] != 0;
}
/**
* Active gas predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not gas is an active phase.
*/
inline bool
gas(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Vapour ] != 0;
}
} // namespace PhaseUsed
/**
* Convenience functions for querying numerical IDs
* ("positions") of active phases.
*/
namespace PhasePos {
/**
* Numerical ID of active water phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of water if
* active, -1 if not.
*/
inline int
water(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::water(pu)) {
p = pu.phase_pos[ BlackoilPhases::Aqua ];
}
return p;
}
/**
* Numerical ID of active oil phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of oil if
* active, -1 if not.
*/
inline int
oil(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::oil(pu)) {
p = pu.phase_pos[ BlackoilPhases::Liquid ];
}
return p;
}
/**
* Numerical ID of active gas phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of gas if
* active, -1 if not.
*/
inline int
gas(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::gas(pu)) {
p = pu.phase_pos[ BlackoilPhases::Vapour ];
}
return p;
}
} // namespace PhasePos
} // namespace Details
/**
* Convert component rates at surface conditions to phase
* (voidage) rates at reservoir conditions.
*
* The conversion uses fluid properties evaluated at average
* hydrocarbon pressure in regions or field.
*
* \tparam Property Fluid property object. Expected to
* feature the formation volume factor functions of the
* BlackoilPropsAdInterface.
*
* \tparam Region Type of a forward region mapping. Expected
* to provide indexed access through \code operator[]()
* \endcode as well as inner types \c value_type, \c
* size_type, and \c const_iterator. Typically \code
* std::vector<int> \endcode.
*/
template <class Property, class Region>
class SurfaceToReservoirVoidage {
public:
/**
* Constructor.
*
* \param[in] props Fluid property object.
*
* \param[in] region Forward region mapping. Often
* corresponds to the "FIPNUM" mapping of an ECLIPSE input
* deck.
*/
SurfaceToReservoirVoidage(const Property& props,
const Region& region)
: props_ (props)
, rmap_ (region)
, repcells_(Details::representative<typename Property::Cells>(rmap_))
, ncells_ (Details::countCells(rmap_))
, p_avg_ (rmap_.numRegions())
, Rmax_ (rmap_.numRegions(), props.numPhases())
{}
/**
* Compute average hydrocarbon pressure and maximum
* dissolution and evaporation at average hydrocarbon
* pressure in all regions in field.
*
* Fluid properties are evaluated at average hydrocarbon
* pressure for purpose of conversion from surface rate to
* reservoir voidage rate.
*
* \param[in] state Dynamic reservoir state.
*/
void
defineState(const BlackoilState& state)
{
averagePressure(state);
calcRmax();
}
/**
* Region identifier.
*
* Integral type.
*/
typedef typename RegionMapping<Region>::RegionId RegionId;
/**
* Compute coefficients for surface-to-reservoir voidage
* conversion.
*
* \tparam Input Type representing contiguous collection
* of component rates at surface conditions. Must support
* direct indexing through \code operator[]()\endcode.
*
* \tparam Coeff Type representing contiguous collection
* of surface-to-reservoir conversion coefficients. Must
* support direct indexing through \code operator[]()
* \endcode.
*
* \param[in] in Single tuple of active component rates at
* surface conditions.
*
* \param[in] r Fluid-in-place region to which the
* component rates correspond.
*
* \param[out] coeff Surface-to-reservoir conversion
* coefficients for all active phases, corresponding to
* input rates \c in in region \c r.
*/
template <class Input,
class Coeff>
void
calcCoeff(const Input& in, const RegionId r, Coeff& coeff)
{
typedef typename Property::V V;
const PhaseUsage& pu = props_.phaseUsage();
const V& p = getRegPress(r);
const typename Property::Cells& c = getRegCell (r);
const int iw = Details::PhasePos::water(pu);
const int io = Details::PhasePos::oil (pu);
const int ig = Details::PhasePos::gas (pu);
std::fill(& coeff[0], & coeff[0] + props_.numPhases(), 0.0);
if (Details::PhaseUsed::water(pu)) {
// q[w]_r = q[w]_s / bw
const V& bw = props_.bWat(p, c);
coeff[iw] = 1.0 / bw(0);
}
const Miscibility& m = calcMiscibility(in, r);
// Determinant of 'R' matrix
const double detR = 1.0 - (m.rs(0) * m.rv(0));
if (Details::PhaseUsed::oil(pu)) {
// q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
const V& bo = props_.bOil(p, m.rs, m.cond, c);
const double den = bo(0) * detR;
coeff[io] += 1.0 / den;
if (Details::PhaseUsed::gas(pu)) {
coeff[ig] -= m.rv(0) / den;
}
}
if (Details::PhaseUsed::gas(pu)) {
// q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
const V& bg = props_.bGas(p, m.rv, m.cond, c);
const double den = bg(0) * detR;
coeff[ig] += 1.0 / den;
if (Details::PhaseUsed::oil(pu)) {
coeff[io] -= m.rs(0) / den;
}
}
}
private:
/**
* Fluid property object.
*/
const Property& props_;
/**
* "Fluid-in-place" region mapping (forward and reverse).
*/
const RegionMapping<Region> rmap_;
/**
* Representative cells in each FIP region.
*/
const typename Property::Cells repcells_;
/**
* Number of cells in each region.
*
* Floating-point type (double) for purpose of average
* pressure calculation.
*/
const Eigen::ArrayXd ncells_;
/**
* Average hydrocarbon pressure in each FIP region.
*/
Eigen::ArrayXd p_avg_;
/**
* Maximum dissolution and evaporation ratios at average
* hydrocarbon pressure.
*
* Size (number of regions)-by-(number of fluid phases).
* Water value is, strictly speaking, wasted if present.
*/
Eigen::ArrayXXd Rmax_;
/**
* Aggregate structure defining fluid miscibility
* conditions in single region with particular input
* surface rates.
*/
struct Miscibility {
Miscibility()
: rs (1)
, rv (1)
, cond(1)
{
rs << 0.0;
rv << 0.0;
}
/**
* Dissolved gas-oil ratio at particular component oil
* and gas rates at surface conditions.
*
* Limited by "RSmax" at average hydrocarbon pressure
* in region.
*/
typename Property::V rs;
/**
* Evaporated oil-gas ratio at particular component oil
* and gas rates at surface conditions.
*
* Limited by "RVmax" at average hydrocarbon pressure
* in region.
*/
typename Property::V rv;
/**
* Fluid condition in representative region cell.
*
* Needed for purpose of FVF evaluation.
*/
std::vector<PhasePresence> cond;
};
/**
* Compute average hydrocarbon pressure in all regions.
*
* \param[in] state Dynamic reservoir state.
*/
void
averagePressure(const BlackoilState& state)
{
p_avg_.setZero();
const std::vector<double>& p = state.pressure();
for (std::vector<double>::size_type
i = 0, n = p.size(); i < n; ++i)
{
p_avg_(rmap_.region(i)) += p[i];
}
p_avg_ /= ncells_;
}
/**
* Compute maximum dissolution and evaporation ratios at
* average hydrocarbon pressure.
*
* Uses the pressure value computed by averagePressure()
* and must therefore be called *after* that method.
*/
void
calcRmax()
{
Rmax_.setZero();
const PhaseUsage& pu = props_.phaseUsage();
if (Details::PhaseUsed::oil(pu) &&
Details::PhaseUsed::gas(pu))
{
const Eigen::ArrayXXd::Index
io = Details::PhasePos::oil(pu),
ig = Details::PhasePos::gas(pu);
// Note: Intentionally does not take capillary
// pressure into account. This facility uses the
// average *hydrocarbon* pressure rather than
// average phase pressure.
Rmax_.col(io) = props_.rsSat(p_avg_, repcells_);
Rmax_.col(ig) = props_.rvSat(p_avg_, repcells_);
}
}
/**
* Compute fluid conditions in particular region for a
* given set of component rates at surface conditions.
*
* \tparam Input Type representing collection of (active)
* component rates at surface conditions. Must support
* direct indexing through \code operator[]()\endcode.
*
* \param[in] in Single tuple of active component rates at
* surface conditions.
*
* \param[in] r Fluid-in-place region to which the
* component rates correspond.
*
* \return Fluid conditions in region \c r corresponding
* to surface component rates \c in.
*/
template <class Input>
Miscibility
calcMiscibility(const Input& in, const RegionId r) const
{
const PhaseUsage& pu = props_.phaseUsage();
const int io = Details::PhasePos::oil(pu);
const int ig = Details::PhasePos::gas(pu);
Miscibility m;
PhasePresence& cond = m.cond[0];
if (Details::PhaseUsed::water(pu)) {
cond.setFreeWater();
}
if (Details::PhaseUsed::oil(pu)) {
cond.setFreeOil();
if (Details::PhaseUsed::gas(pu)) {
const double rsmax = Rmax_(r, io);
const double rs =
(0.0 < std::abs(in[io]))
? in[ig] / in[io]
: (0.0 < std::abs(in[ig])) ? rsmax : 0.0;
if (rsmax < rs) {
cond.setFreeGas();
}
m.rs(0) = std::min(rs, rsmax);
}
}
if (Details::PhaseUsed::gas(pu)) {
if (! Details::PhaseUsed::oil(pu)) {
// Oil *NOT* active -- not really supported.
cond.setFreeGas();
}
if (Details::PhaseUsed::oil(pu)) {
const double rvmax = Rmax_(r, ig);
const double rv =
(0.0 < std::abs(in[ig]))
? (in[io] / in[ig])
: (0.0 < std::abs(in[io])) ? rvmax : 0.0;
m.rv(0) = std::min(rv, rvmax);
}
}
return m;
}
/**
* Retrieve average hydrocarbon pressure in region.
*
* \param[in] r Particular region.
*
* \return Average hydrocarbon pressure in region \c r.
*/
typename Property::V
getRegPress(const RegionId r) const
{
typename Property::V p(1);
p << p_avg_(r);
return p;
}
/**
* 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, repcells_[r]);
return c;
}
};
} // namespace RateConverter
} // namespace Opm
#endif /* OPM_RATECONVERTER_HPP_HEADER_INCLUDED */

5
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
View File

@ -33,7 +33,6 @@ namespace Opm
class BlackoilPropsAdInterface;
class RockCompressibility;
class DerivedGeology;
class WellsManager;
class NewtonIterationBlackoilInterface;
class SimulatorTimer;
class BlackoilState;
@ -66,11 +65,13 @@ namespace Opm
/// segregation is ignored).
///
/// \param[in] grid grid data structure
/// \param[in] geo derived geological properties
/// \param[in] props fluid and rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] well_manager well manager, may manage no (null) wells
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
/// \param[in] disgas true for dissolved gas option
/// \param[in] vapoil true for vaporized oil option
/// \param[in] eclipse_state
/// \param[in] output_writer
SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,

251
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
View File

@ -26,8 +26,11 @@
#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/io/eclipse/EclipseWriter.hpp>
@ -38,22 +41,31 @@
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/grid/ColumnExtract.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
#include <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>
#include <algorithm>
#include <cstddef>
#include <cassert>
#include <functional>
#include <memory>
#include <numeric>
#include <fstream>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
namespace Opm
{
@ -78,15 +90,17 @@ namespace Opm
private:
// Data.
typedef RateConverter::
SurfaceToReservoirVoidage< BlackoilPropsAdInterface,
std::vector<int> > RateConverterType;
const parameter::ParameterGroup param_;
// Parameters for output.
bool output_;
bool output_vtk_;
std::string output_dir_;
int output_interval_;
// Parameters for well control
bool check_well_controls_;
int max_well_control_iterations_;
// Observed objects.
const Grid& grid_;
BlackoilPropsAdInterface& props_;
@ -103,6 +117,13 @@ namespace Opm
std::shared_ptr<EclipseState> eclipse_state_;
// output_writer
EclipseWriter& output_writer_;
RateConverterType rateConverter_;
void
computeRESV(const std::size_t step,
const Wells* wells,
const BlackoilState& x,
WellStateFullyImplicitBlackoil& xw);
};
@ -219,7 +240,8 @@ namespace Opm
has_disgas_(has_disgas),
has_vapoil_(has_vapoil),
eclipse_state_(eclipse_state),
output_writer_(output_writer)
output_writer_(output_writer),
rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0))
{
// For output.
output_ = param.getDefault("output", true);
@ -237,10 +259,6 @@ namespace Opm
output_interval_ = param.getDefault("output_interval", 1);
}
// Well control related init.
check_well_controls_ = param.getDefault("check_well_controls", false);
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
// Misc init.
const int num_cells = AutoDiffGrid::numCells(grid);
allcells_.resize(num_cells);
@ -308,6 +326,9 @@ namespace Opm
props_.updateSatOilMax(state.saturation());
props_.updateSatHyst(state.saturation(), allcells_);
// Compute reservoir volumes for RESV controls.
computeRESV(timer.currentStepNum(), wells, state, well_state);
// Run a single step of the solver.
solver_timer.start();
FullyImplicitBlackoilSolver<T> solver(param_, grid_, props_, geo_, rock_comp_props_, *wells, solver_, has_disgas_, has_vapoil_);
@ -349,5 +370,211 @@ namespace Opm
return report;
}
namespace SimFIBODetails {
typedef std::unordered_map<std::string, WellConstPtr> WellMap;
inline WellMap
mapWells(const std::vector<WellConstPtr>& wells)
{
WellMap wmap;
for (std::vector<WellConstPtr>::const_iterator
w = wells.begin(), e = wells.end();
w != e; ++w)
{
wmap.insert(std::make_pair((*w)->name(), *w));
}
return wmap;
}
inline int
resv_control(const WellControls* ctrl)
{
int i, n = well_controls_get_num(ctrl);
bool match = false;
for (i = 0; (! match) && (i < n); ++i) {
match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE;
}
if (! match) { i = 0; }
return i - 1; // -1 if no match, undo final "++" otherwise
}
inline bool
is_resv_prod(const Wells& wells,
const int w)
{
return ((wells.type[w] == PRODUCER) &&
(0 <= resv_control(wells.ctrls[w])));
}
inline bool
is_resv_prod(const WellMap& wmap,
const std::string& name,
const std::size_t step)
{
bool match = false;
WellMap::const_iterator i = wmap.find(name);
if (i != wmap.end()) {
WellConstPtr wp = i->second;
match = (wp->isProducer(step) &&
wp->getProductionProperties(step)
.hasProductionControl(WellProducer::RESV));
}
return match;
}
inline std::vector<int>
resvProducers(const Wells& wells,
const std::size_t step,
const WellMap& wmap)
{
std::vector<int> resv_prod;
for (int w = 0, nw = wells.number_of_wells; w < nw; ++w) {
if (is_resv_prod(wells, w) ||
((wells.name[w] != 0) &&
is_resv_prod(wmap, wells.name[w], step)))
{
resv_prod.push_back(w);
}
}
return resv_prod;
}
inline void
historyRates(const PhaseUsage& pu,
const WellProductionProperties& p,
std::vector<double>& rates)
{
assert (! p.predictionMode);
assert (rates.size() ==
std::vector<double>::size_type(pu.num_phases));
if (pu.phase_used[ BlackoilPhases::Aqua ]) {
const std::vector<double>::size_type
i = pu.phase_pos[ BlackoilPhases::Aqua ];
rates[i] = p.WaterRate;
}
if (pu.phase_used[ BlackoilPhases::Liquid ]) {
const std::vector<double>::size_type
i = pu.phase_pos[ BlackoilPhases::Liquid ];
rates[i] = p.OilRate;
}
if (pu.phase_used[ BlackoilPhases::Vapour ]) {
const std::vector<double>::size_type
i = pu.phase_pos[ BlackoilPhases::Vapour ];
rates[i] = p.GasRate;
}
}
} // namespace SimFIBODetails
template <class T>
void
SimulatorFullyImplicitBlackoil<T>::
Impl::computeRESV(const std::size_t step,
const Wells* wells,
const BlackoilState& x,
WellStateFullyImplicitBlackoil& xw)
{
typedef SimFIBODetails::WellMap WellMap;
const std::vector<WellConstPtr>& w_ecl = eclipse_state_->getSchedule()->getWells(step);
const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
const std::vector<int>& resv_prod =
SimFIBODetails::resvProducers(*wells, step, wmap);
if (! resv_prod.empty()) {
const PhaseUsage& pu = props_.phaseUsage();
const std::vector<double>::size_type np = props_.numPhases();
rateConverter_.defineState(x);
std::vector<double> distr (np);
std::vector<double> hrates(np);
std::vector<double> prates(np);
for (std::vector<int>::const_iterator
rp = resv_prod.begin(), e = resv_prod.end();
rp != e; ++rp)
{
WellControls* ctrl = wells->ctrls[*rp];
// RESV control mode, all wells
{
const int rctrl = SimFIBODetails::resv_control(ctrl);
if (0 <= rctrl) {
const std::vector<double>::size_type off = (*rp) * np;
// Convert to positive rates to avoid issues
// in coefficient calculations.
std::transform(xw.wellRates().begin() + (off + 0*np),
xw.wellRates().begin() + (off + 1*np),
prates.begin(), std::negate<double>());
const int fipreg = 0; // Hack. Ignore FIP regions.
rateConverter_.calcCoeff(prates, fipreg, distr);
well_controls_iset_distr(ctrl, rctrl, & distr[0]);
}
}
// RESV control, WCONHIST wells. A bit of duplicate
// work, regrettably.
if (wells->name[*rp] != 0) {
WellMap::const_iterator i = wmap.find(wells->name[*rp]);
if (i != wmap.end()) {
WellConstPtr wp = i->second;
const WellProductionProperties& p =
wp->getProductionProperties(step);
if (! p.predictionMode) {
// History matching (WCONHIST/RESV)
SimFIBODetails::historyRates(pu, p, hrates);
const int fipreg = 0; // Hack. Ignore FIP regions.
rateConverter_.calcCoeff(hrates, fipreg, distr);
// WCONHIST/RESV target is sum of all
// observed phase rates translated to
// reservoir conditions. Recall sign
// convention: Negative for producers.
const double target =
- std::inner_product(distr.begin(), distr.end(),
hrates.begin(), 0.0);
well_controls_clear(ctrl);
well_controls_assert_number_of_phases(ctrl, int(np));
const int ok =
well_controls_add_new(RESERVOIR_RATE, target,
& distr[0], ctrl);
if (ok != 0) {
xw.currentControls()[*rp] = 0;
well_controls_set_current(ctrl, 0);
}
}
}
}
}
}
}
} // namespace Opm

124
tests/test_rateconverter.cpp Normal file
View File

@ -0,0 +1,124 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define BOOST_TEST_MODULE VoidageRateConversionTest
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <boost/test/unit_test.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
struct SetupSimple {
SetupSimple()
{
Opm::ParserPtr parser(new Opm::Parser());
deck = parser->parseFile("fluid.data");
eclState.reset(new Opm::EclipseState(deck));
param.disableOutput();
param.insertParameter("init_rock" , "false" );
param.insertParameter("threephase_model", "simple");
param.insertParameter("pvt_tab_size" , "0" );
param.insertParameter("sat_tab_size" , "0" );
}
Opm::parameter::ParameterGroup param;
Opm::DeckConstPtr deck;
Opm::EclipseStateConstPtr eclState;
};
template <class Setup>
struct TestFixture : public Setup
{
TestFixture()
: Setup()
, grid (deck)
, props(deck, eclState, *grid.c_grid(), param,
param.getDefault("init_rock", false))
{
}
using Setup::param;
using Setup::deck;
using Setup::eclState;
Opm::GridManager grid;
Opm::BlackoilPropertiesFromDeck props;
};
BOOST_FIXTURE_TEST_CASE(Construction, TestFixture<SetupSimple>)
{
typedef std::vector<int> Region;
typedef Opm::BlackoilPropsAd Props;
typedef Opm::RateConverter::
SurfaceToReservoirVoidage<Props, Region> RCvrt;
Region reg{ 0 };
Props ad_props(props);
RCvrt cvrt(ad_props, reg);
}
BOOST_FIXTURE_TEST_CASE(TwoPhaseII, TestFixture<SetupSimple>)
{
// Immiscible and incompressible two-phase fluid
typedef std::vector<int> Region;
typedef Opm::BlackoilPropsAd Props;
typedef Opm::RateConverter::
SurfaceToReservoirVoidage<Props, Region> RCvrt;
Region reg{ 0 };
Props ad_props(props);
RCvrt cvrt(ad_props, reg);
Opm::BlackoilState x;
x.init(*grid.c_grid(), 2);
cvrt.defineState(x);
std::vector<double> qs{1.0e3, 1.0e1};
std::vector<double> coeff(qs.size(), 0.0);
// Immiscible and incompressible: All coefficients are one (1),
// irrespective of actual surface rates.
cvrt.calcCoeff(qs, 0, coeff);
BOOST_CHECK_CLOSE(coeff[0], 1.0, 1.0e-6);
BOOST_CHECK_CLOSE(coeff[1], 1.0, 1.0e-6);
}