OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge remote branch 'remotes/opm/master' into relpermDiagnostics

Browse Source 
Conflicts:
	CMakeLists_files.cmake
This commit is contained in:
Liu Ming 2015-12-08 16:08:02 +08:00
commit ef203a1d04
20 changed files with 1247 additions and 159 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists_files.cmake
View File

@ -201,6 +201,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_anisotropiceikonal.cpp
tests/test_stoppedwells.cpp
tests/test_relpermdiagnostics.cpp
tests/test_norne_pvt.cpp
)
# originally generated with the command:
@ -235,7 +236,8 @@ list (APPEND TEST_DATA_FILES
tests/CORNERPOINT_ACTNUM.DATA
tests/wells_stopped.data
tests/relpermDiagnostics.DATA
)
tests/norne_pvt.data
)
# originally generated with the command:
# find tutorials examples -name '*.c*' -printf '\t%p\n' | sort

4
CTestConfig.cmake
View File

@ -2,6 +2,6 @@
set(CTEST_PROJECT_NAME "${${project}_NAME}")
set(CTEST_NIGHTLY_START_TIME "01:00:00 UTC")
set(CTEST_DROP_METHOD "http")
set(CTEST_DROP_SITE "opm-project.org")
set(CTEST_DROP_LOCATION "/CDash/submit.php?project=${${project}_NAME}")
set(CTEST_DROP_SITE "cdash.opm-project.org")
set(CTEST_DROP_LOCATION "/submit.php?project=${${project}_NAME}")
set(CTEST_DROP_SITE_CDASH TRUE)

5
opm/core/grid/GridHelpers.hpp
View File

@ -20,8 +20,13 @@
*/
#ifndef OPM_CORE_GRIDHELPERS_HEADER_INCLUDED
#define OPM_CORE_GRIDHELPERS_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/range/iterator_range.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
namespace Opm
{

6
opm/core/linalg/ParallelIstlInformation.hpp
View File

@ -144,8 +144,9 @@ public:
void copyOwnerToAll (const T& source, T& dest) const
{
typedef Dune::Combine<Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner>,Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::overlap>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> OwnerOverlapSet;
typedef Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner> OwnerSet;
typedef Dune::Combine<OwnerOverlapSet, Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::copy>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> AllSet;
OwnerOverlapSet sourceFlags;
OwnerSet sourceFlags;
AllSet destFlags;
Dune::Interface interface(communicator_);
if( !remoteIndices_->isSynced() )
@ -159,7 +160,7 @@ public:
communicator.free();
}
template<class T>
void updateOwnerMask(const T& container) const
const std::vector<double>& updateOwnerMask(const T& container) const
{
if( ! indexSet_ )
{
@ -176,6 +177,7 @@ public:
}
}
}
return ownerMask_;
}
/// \brief Compute one or more global reductions.
///

6
opm/core/props/BlackoilPropertiesFromDeck.hpp
View File

@ -233,6 +233,12 @@ namespace Opm
const double pcow,
double & swat);
// return a reference to the "raw" PVT fluid object for a phase.
const PvtInterface& pvt(int phaseIdx) const
{
return pvt_.pvt(phaseIdx);
}
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{

6
opm/core/props/pvt/BlackoilPvtProperties.hpp
View File

@ -117,6 +117,12 @@ namespace Opm
double* output_R,
double* output_dRdp) const;
// return a reference to the "raw" PVT fluid object for a phase.
const PvtInterface& pvt(int phaseIdx) const
{
return *props_[phaseIdx];
}
private:
// Disabling copying (just to avoid surprises, since we use shared_ptr).
BlackoilPvtProperties(const BlackoilPvtProperties&);

8
opm/core/props/rock/RockFromDeck.cpp
View File

@ -72,6 +72,14 @@ namespace Opm
{
}
RockFromDeck::RockFromDeck(std::size_t number_of_cells)
: porosity_(number_of_cells, 0),
// full permeability tensor in 3D stores 9 scalars
permeability_(number_of_cells*9, 0.0),
permfield_valid_(number_of_cells, false)
{
}
void RockFromDeck::init(Opm::EclipseStateConstPtr eclState,
int number_of_cells, const int* global_cell,
const int* cart_dims)

8
opm/core/props/rock/RockFromDeck.hpp
View File

@ -31,10 +31,16 @@ namespace Opm
class RockFromDeck
{
// BlackoilPropsDataHandle needs mutable
// access to porosity and permeability
friend class BlackoilPropsDataHandle;
public:
/// Default constructor.
RockFromDeck();
/// Creates rock properties with zero porosity and permeability
/// \param number_of_cells The number of cells
explicit RockFromDeck(std::size_t number_of_cells);
/// Initialize from deck and cell mapping.
/// \param eclState The EclipseState (processed deck) produced by the opm-parser code
/// \param number_of_cells The number of cells in the grid.

1
opm/core/simulator/AdaptiveTimeStepping.hpp
View File

@ -43,7 +43,6 @@ namespace Opm {
//! \param pinfo The information about the data distribution
//! and communication for a parallel run.
AdaptiveTimeStepping( const parameter::ParameterGroup& param,
const boost::any& pinfo=boost::any(),
const bool terminal_output = true );
/** \brief step method that acts like the solver::step method

39
opm/core/simulator/AdaptiveTimeStepping_impl.hpp
View File

@ -31,11 +31,35 @@
namespace Opm {
namespace detail
{
template <class Solver, class State>
class SolutionTimeErrorSolverWrapper : public RelativeChangeInterface
{
const Solver& solver_;
const State& previous_;
const State& current_;
public:
SolutionTimeErrorSolverWrapper( const Solver& solver,
const State& previous,
const State& current )
: solver_( solver ),
previous_( previous ),
current_( current )
{}
/// return || u^n+1 - u^n || / || u^n+1 ||
double relativeChange() const
{
return solver_.model().relativeChange( previous_, current_ );
}
};
}
// AdaptiveTimeStepping
//---------------------
AdaptiveTimeStepping::AdaptiveTimeStepping( const parameter::ParameterGroup& param,
const boost::any& parallel_information,
const bool terminal_output )
: timeStepControl_()
, restart_factor_( param.getDefault("solver.restartfactor", double(0.33) ) )
@ -56,12 +80,12 @@ namespace Opm {
const double tol = param.getDefault("timestep.control.tol", double(1e-1) );
if( control == "pid" ) {
timeStepControl_ = TimeStepControlType( new PIDTimeStepControl( tol, parallel_information ) );
timeStepControl_ = TimeStepControlType( new PIDTimeStepControl( tol ) );
}
else if ( control == "pid+iteration" )
{
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations );
timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol, parallel_information ) );
timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
}
else if ( control == "iterationcount" )
{
@ -130,9 +154,6 @@ namespace Opm {
// get current delta t
const double dt = substepTimer.currentStepLength() ;
// initialize time step control in case current state is needed later
timeStepControl_->initialize( state );
if( timestep_verbose_ )
{
std::cout <<"Substep( " << substepTimer.currentStepNum() << " ), try with stepsize "
@ -172,9 +193,13 @@ namespace Opm {
// advance by current dt
++substepTimer;
// create object to compute the time error, simply forwards the call to the model
detail::SolutionTimeErrorSolverWrapper< Solver, State >
relativeChange( solver, last_state, state );
// compute new time step estimate
double dtEstimate =
timeStepControl_->computeTimeStepSize( dt, linearIterations, state );
timeStepControl_->computeTimeStepSize( dt, linearIterations, relativeChange );
// limit the growth of the timestep size by the growth factor
dtEstimate = std::min( dtEstimate, double(max_growth_ * dt) );

2
opm/core/simulator/EquilibrationHelpers.hpp
View File

@ -350,7 +350,7 @@ namespace Opm
const double temp,
const double sat_oil = 0.0 ) const
{
if (sat_oil > 0.0) {
if (std::abs(sat_oil) > 1e-16) {
return satRv(press, temp);
} else {
return std::min(satRv(press, temp), linearInterpolation(depth_, rv_, depth));

53
opm/core/simulator/TimeStepControl.cpp
View File

@ -55,7 +55,7 @@ namespace Opm
}
double SimpleIterationCountTimeStepControl::
computeTimeStepSize( const double dt, const int iterations, const SimulatorState& /* state */ ) const
computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */ ) const
{
double dtEstimate = dt ;
@ -82,58 +82,24 @@ namespace Opm
//
////////////////////////////////////////////////////////
PIDTimeStepControl::PIDTimeStepControl( const double tol, const boost::any& pinfo,
PIDTimeStepControl::PIDTimeStepControl( const double tol,
const bool verbose )
: p0_()
, sat0_()
, tol_( tol )
: tol_( tol )
, errors_( 3, tol_ )
, verbose_( verbose )
, parallel_information_(pinfo)
{}
void PIDTimeStepControl::initialize( const SimulatorState& state )
{
// store current state for later time step computation
p0_ = state.pressure();
sat0_ = state.saturation();
}
double PIDTimeStepControl::
computeTimeStepSize( const double dt, const int /* iterations */, const SimulatorState& state ) const
computeTimeStepSize( const double dt, const int /* iterations */, const RelativeChangeInterface& relChange ) const
{
const std::size_t pSize = p0_.size();
assert( state.pressure().size() == pSize );
const std::size_t satSize = sat0_.size();
assert( state.saturation().size() == satSize );
// compute u^n - u^n+1
for( std::size_t i=0; i<pSize; ++i ) {
p0_[ i ] -= state.pressure()[ i ];
}
for( std::size_t i=0; i<satSize; ++i ) {
sat0_[ i ] -= state.saturation()[ i ];
}
// compute || u^n - u^n+1 ||
const double stateOld = euclidianNormSquared( p0_.begin(), p0_.end() ) +
euclidianNormSquared( sat0_.begin(), sat0_.end(),
state.numPhases() );
// compute || u^n+1 ||
const double stateNew = euclidianNormSquared( state.pressure().begin(), state.pressure().end() ) +
euclidianNormSquared( state.saturation().begin(), state.saturation().end(),
state.numPhases() );
// shift errors
for( int i=0; i<2; ++i ) {
errors_[ i ] = errors_[i+1];
}
// store new error
const double error = stateOld / stateNew;
errors_[ 2 ] = error ;
const double error = relChange.relativeChange();
errors_[ 2 ] = error;
if( error > tol_ )
{
@ -169,16 +135,15 @@ namespace Opm
PIDAndIterationCountTimeStepControl::
PIDAndIterationCountTimeStepControl( const int target_iterations,
const double tol,
const boost::any& pinfo,
const bool verbose)
: BaseType( tol, pinfo, verbose )
: BaseType( tol, verbose )
, target_iterations_( target_iterations )
{}
double PIDAndIterationCountTimeStepControl::
computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const
computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relChange ) const
{
double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, state );
double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, relChange );
// further reduce step size if to many iterations were used
if( iterations > target_iterations_ )

57
opm/core/simulator/TimeStepControl.hpp
View File

@ -26,7 +26,6 @@
#include <boost/any.hpp>
#include <boost/range/iterator_range.hpp>
#include <opm/core/simulator/TimeStepControlInterface.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
namespace Opm
{
@ -49,7 +48,7 @@ namespace Opm
const bool verbose = false);
/// \brief \copydoc TimeStepControlInterface::computeTimeStepSize
double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const;
double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */ ) const;
protected:
const int target_iterations_;
@ -78,64 +77,18 @@ namespace Opm
/// \brief constructor
/// \param tol tolerance for the relative changes of the numerical solution to be accepted
/// in one time step (default is 1e-3)
/// \paramm pinfo The information about the parallel information. Needed to
/// compute parallel scalarproducts.
/// \param verbose if true get some output (default = false)
PIDTimeStepControl( const double tol = 1e-3,
const boost::any& pinfo = boost::any(),
const bool verbose = false );
/// \brief \copydoc TimeStepControlInterface::initialize
void initialize( const SimulatorState& state );
/// \brief \copydoc TimeStepControlInterface::computeTimeStepSize
double computeTimeStepSize( const double dt, const int /* iterations */, const SimulatorState& state ) const;
double computeTimeStepSize( const double dt, const int /* iterations */, const RelativeChangeInterface& relativeChange ) const;
protected:
template <class Iterator>
double euclidianNormSquared( Iterator it, const Iterator end, int num_components = 1 ) const
{
static_cast<void>(num_components); // Suppress warning in the serial case.
#if HAVE_MPI
if ( parallel_information_.type() == typeid(ParallelISTLInformation) )
{
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>(parallel_information_);
int size_per_component = (end - it) / num_components;
assert((end - it) == num_components * size_per_component);
double component_product = 0.0;
for( int i = 0; i < num_components; ++i )
{
auto component_container =
boost::make_iterator_range(it + i * size_per_component,
it + (i + 1) * size_per_component);
info.computeReduction(component_container,
Opm::Reduction::makeInnerProductFunctor<double>(),
component_product);
}
return component_product;
}
else
#endif
{
double product = 0.0 ;
for( ; it != end; ++it ) {
product += ( *it * *it );
}
return product;
}
}
protected:
mutable std::vector<double> p0_;
mutable std::vector<double> sat0_;
const double tol_;
mutable std::vector< double > errors_;
const bool verbose_;
private:
const boost::any parallel_information_;
};
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
@ -152,17 +105,13 @@ namespace Opm
/// \param target_iterations number of desired iterations per time step
/// \param tol tolerance for the relative changes of the numerical solution to be accepted
/// in one time step (default is 1e-3)
// \param maxgrowth max growth factor for new time step in relation of old time step (default = 3.0)
/// \paramm pinfo The information about the parallel information. Needed to
/// compute parallel scalarproducts.
/// \param verbose if true get some output (default = false)
PIDAndIterationCountTimeStepControl( const int target_iterations = 20,
const double tol = 1e-3,
const boost::any& = boost::any(),
const bool verbose = false);
/// \brief \copydoc TimeStepControlInterface::computeTimeStepSize
double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const;
double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relativeChange ) const;
protected:
const int target_iterations_;

38
opm/core/simulator/TimeStepControlInterface.hpp
View File

@ -1,5 +1,5 @@
/*
Copyright 2014 IRIS AS
Copyright 2014 IRIS AS
This file is part of the Open Porous Media project (OPM).
@ -19,31 +19,45 @@
#ifndef OPM_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED
#define OPM_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED
#include <opm/core/simulator/SimulatorState.hpp>
#include <opm/core/simulator/SimulatorState.hpp>
namespace Opm
{
///////////////////////////////////////////////////////////////////
///
/// TimeStepControlInterface
///
/// RelativeChangeInterface
///
///////////////////////////////////////////////////////////////////
class TimeStepControlInterface
class RelativeChangeInterface
{
protected:
protected:
RelativeChangeInterface() {}
public:
/// \return || u^n+1 - u^n || / || u^n+1 ||
virtual double relativeChange() const = 0;
/// virtual destructor (empty)
virtual ~RelativeChangeInterface() {}
};
///////////////////////////////////////////////////////////////////
///
/// TimeStepControlInterface
///
///////////////////////////////////////////////////////////////////
class TimeStepControlInterface
{
protected:
TimeStepControlInterface() {}
public:
/// \param state simulation state before computing update in the solver (default is empty)
virtual void initialize( const SimulatorState& /*state*/ ) {}
/// compute new time step size suggestions based on the PID controller
/// \param dt time step size used in the current step
/// \param iterations number of iterations used (linear/nonlinear)
/// \param state new solution state
/// \param iterations number of iterations used (linear/nonlinear)
/// \param timeError object to compute || u^n+1 - u^n || / || u^n+1 ||
///
/// \return suggested time step size for the next step
virtual double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& ) const = 0;
virtual double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relativeChange ) const = 0;
/// virtual destructor (empty)
virtual ~TimeStepControlInterface () {}

29
opm/core/simulator/initStateEquil_impl.hpp
View File

@ -313,7 +313,6 @@ namespace Opm
const CellRange& cells,
std::vector<double>& p )
{
const int nd = UgGridHelpers::dimensions(G);
enum { up = 0, down = 1 };
@ -324,7 +323,7 @@ namespace Opm
{
assert (c < p.size());
const double z = UgGridHelpers::cellCentroidCoordinate(G, *ci, nd-1);
const double z = UgGridHelpers::cellCenterDepth(G, *ci);
p[c] = (z < split) ? f[up](z) : f[down](z);
}
}
@ -366,9 +365,9 @@ namespace Opm
typedef Details::RK4IVP<ODE> WPress;
std::array<WPress,2> wpress = {
{
WPress(drho, up , p0, 100)
WPress(drho, up , p0, 2000)
,
WPress(drho, down, p0, 100)
WPress(drho, down, p0, 2000)
}
};
@ -426,9 +425,9 @@ namespace Opm
typedef Details::RK4IVP<ODE> OPress;
std::array<OPress,2> opress = {
{
OPress(drho, up , p0, 100)
OPress(drho, up , p0, 2000)
,
OPress(drho, down, p0, 100)
OPress(drho, down, p0, 2000)
}
};
@ -488,9 +487,9 @@ namespace Opm
typedef Details::RK4IVP<ODE> GPress;
std::array<GPress,2> gpress = {
{
GPress(drho, up , p0, 100)
GPress(drho, up , p0, 2000)
,
GPress(drho, down, p0, 100)
GPress(drho, down, p0, 2000)
}
};
@ -703,10 +702,8 @@ namespace Opm
double sw = 0.0;
if (water) {
if (isConstPc(props,waterpos,cell)){
const int nd = UgGridHelpers::dimensions(G);
const double cellDepth = UgGridHelpers::cellCentroidCoordinate(G,
cell,
nd-1);
const double cellDepth = UgGridHelpers::cellCenterDepth(G,
cell);
sw = satFromDepth(props,cellDepth,reg.zwoc(),waterpos,cell,false);
phase_saturations[waterpos][local_index] = sw;
}
@ -725,10 +722,8 @@ namespace Opm
double sg = 0.0;
if (gas) {
if (isConstPc(props,gaspos,cell)){
const int nd = UgGridHelpers::dimensions(G);
const double cellDepth = UgGridHelpers::cellCentroidCoordinate(G,
cell,
nd-1);
const double cellDepth = UgGridHelpers::cellCenterDepth(G,
cell);
sg = satFromDepth(props,cellDepth,reg.zgoc(),gaspos,cell,true);
phase_saturations[gaspos][local_index] = sg;
}
@ -829,7 +824,7 @@ namespace Opm
std::vector<double> rs(cells.size());
int count = 0;
for (auto it = cells.begin(); it != cells.end(); ++it, ++count) {
const double depth = UgGridHelpers::cellCentroidCoordinate(grid, *it, 2);
const double depth = UgGridHelpers::cellCenterDepth(grid, *it);
rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]);
}
return rs;

293
opm/core/utility/thresholdPressures.hpp
View File

@ -20,14 +20,278 @@
#ifndef OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
#define OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <vector>
#include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
#include <opm/parser/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
#include <opm/parser/eclipse/Parser/ParseMode.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
namespace Opm
{
/// \brief Compute the maximum gravity corrected pressure difference of all
/// equilibration regions given a reservoir state.
/// \tparam Grid Type of grid object (UnstructuredGrid or CpGrid).
/// \param[out] maxDp The resulting pressure difference between equilibration regions
/// \param[in] deck Input deck, EQLOPTS and THPRES are accessed from it.
/// \param[in] eclipseState Processed eclipse state, EQLNUM is accessed from it.
/// \param[in] grid The grid to which the thresholds apply.
/// \param[in] initialState The state of the reservoir
/// \param[in] props The object which calculates fluid properties
/// \param[in] gravity The gravity constant
template <class Grid>
void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
const DeckConstPtr& deck,
EclipseStateConstPtr eclipseState,
const Grid& grid,
const BlackoilState& initialState,
const BlackoilPropertiesFromDeck& props,
const double gravity)
{
const PhaseUsage& pu = props.phaseUsage();
std::shared_ptr<GridProperty<int>> eqlnum = eclipseState->getIntGridProperty("EQLNUM");
const auto& eqlnumData = eqlnum->getData();
const int numPhases = initialState.numPhases();
const int numCells = UgGridHelpers::numCells(grid);
const int numPvtRegions = deck->getKeyword("TABDIMS")->getRecord(0)->getItem("NTPVT")->getInt(0);
// retrieve the minimum (residual!?) and the maximum saturations for all cells
std::vector<double> minSat(numPhases*numCells);
std::vector<double> maxSat(numPhases*numCells);
std::vector<int> allCells(numCells);
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
allCells[cellIdx] = cellIdx;
}
props.satRange(numCells, allCells.data(), minSat.data(), maxSat.data());
// retrieve the surface densities
std::vector<std::vector<double> > surfaceDensity(numPvtRegions);
Opm::DeckKeywordConstPtr densityKw = deck->getKeyword("DENSITY");
for (int regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx) {
surfaceDensity[regionIdx].resize(numPhases);
if (pu.phase_used[BlackoilPhases::Aqua]) {
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
surfaceDensity[regionIdx][wpos] =
densityKw->getRecord(regionIdx)->getItem("WATER")->getSIDouble(0);
}
if (pu.phase_used[BlackoilPhases::Liquid]) {
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
surfaceDensity[regionIdx][opos] =
densityKw->getRecord(regionIdx)->getItem("OIL")->getSIDouble(0);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
surfaceDensity[regionIdx][gpos] =
densityKw->getRecord(regionIdx)->getItem("GAS")->getSIDouble(0);
}
}
// retrieve the PVT region of each cell. note that we need c++ instead of
// Fortran indices.
const int* gc = UgGridHelpers::globalCell(grid);
std::vector<int> pvtRegion(numCells);
const auto& cartPvtRegion = eclipseState->getIntGridProperty("PVTNUM")->getData();
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
const int cartCellIdx = gc ? gc[cellIdx] : cellIdx;
pvtRegion[cellIdx] = std::max(0, cartPvtRegion[cartCellIdx] - 1);
}
// compute the initial "phase presence" of each cell (required to calculate
// the inverse formation volume factors
std::vector<PhasePresence> cond(numCells);
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
const double sw = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua]];
const double so = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]];
const double sg = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]];
if (pu.phase_used[BlackoilPhases::Aqua] && sw > 0.0) {
cond[cellIdx].setFreeWater();
}
if (pu.phase_used[BlackoilPhases::Liquid] && so > 0.0) {
cond[cellIdx].setFreeOil();
}
if (pu.phase_used[BlackoilPhases::Vapour] && sg > 0.0) {
cond[cellIdx].setFreeGas();
}
}
// calculate the initial fluid densities for the gravity correction.
std::vector<double> b(numCells);
std::vector<std::vector<double>> rho(numPhases);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
rho[phaseIdx].resize(numCells);
}
// compute the capillary pressures of the active phases
std::vector<double> capPress(numCells*numPhases);
std::vector<int> cellIdxArray(numCells);
for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
cellIdxArray[cellIdx] = cellIdx;
}
props.capPress(numCells, initialState.saturation().data(), cellIdxArray.data(), capPress.data(), NULL);
// compute the absolute pressure of each active phase: for some reason, E100
// defines the capillary pressure for the water phase as p_o - p_w while it
// uses p_g - p_o for the gas phase. (it would be more consistent to use the
// oil pressure as reference for both the other phases.) probably this is
// done to always have a positive number for the capillary pressure (as long
// as the medium is hydrophilic)
std::vector<std::vector<double> > phasePressure(numPhases);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
phasePressure[phaseIdx].resize(numCells);
}
for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
// we currently hard-code the oil phase as the reference phase!
assert(pu.phase_used[BlackoilPhases::Liquid]);
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
phasePressure[opos][cellIdx] = initialState.pressure()[cellIdx];
if (pu.phase_used[BlackoilPhases::Aqua]) {
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
phasePressure[wpos][cellIdx] =
initialState.pressure()[cellIdx]
+ (capPress[cellIdx*numPhases + opos] - capPress[cellIdx*numPhases + wpos]);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
phasePressure[gpos][cellIdx] =
initialState.pressure()[cellIdx]
+ (capPress[cellIdx*numPhases + gpos] - capPress[cellIdx*numPhases + opos]);
}
}
// calculate the inverse formation volume factors for the active phases and each cell
if (pu.phase_used[BlackoilPhases::Aqua]) {
std::vector<double> dummy(numCells*BlackoilPhases::MaxNumPhases);
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const PvtInterface& pvtw = props.pvt(wpos);
pvtw.b(numCells,
pvtRegion.data(),
phasePressure[wpos].data(),
initialState.temperature().data(),
initialState.gasoilratio().data(),
cond.data(),
b.data(),
dummy.data(),
dummy.data());
for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
rho[wpos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][wpos]*b[cellIdx];
}
}
if (pu.phase_used[BlackoilPhases::Liquid]) {
std::vector<double> dummy(numCells*BlackoilPhases::MaxNumPhases);
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
const PvtInterface& pvto = props.pvt(opos);
pvto.b(numCells,
pvtRegion.data(),
phasePressure[opos].data(),
initialState.temperature().data(),
initialState.gasoilratio().data(),
cond.data(),
b.data(),
dummy.data(),
dummy.data());
for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
rho[opos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][opos]*b[cellIdx];
if (pu.phase_used[BlackoilPhases::Vapour]) {
int gpos = pu.phase_pos[BlackoilPhases::Vapour];
rho[opos][cellIdx] +=
surfaceDensity[pvtRegion[cellIdx]][gpos]*initialState.gasoilratio()[cellIdx]*b[cellIdx];
}
}
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
std::vector<double> dummy(numCells*BlackoilPhases::MaxNumPhases);
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const PvtInterface& pvtg = props.pvt(gpos);
pvtg.b(numCells,
pvtRegion.data(),
phasePressure[gpos].data(),
initialState.temperature().data(),
initialState.rv().data(),
cond.data(),
b.data(),
dummy.data(),
dummy.data());
for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
rho[gpos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][gpos]*b[cellIdx];
if (pu.phase_used[BlackoilPhases::Liquid]) {
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
rho[gpos][cellIdx] +=
surfaceDensity[pvtRegion[cellIdx]][opos]*initialState.rv()[cellIdx]*b[cellIdx];
}
}
}
// Calculate the maximum pressure potential difference between all PVT region
// transitions of the initial solution.
const int num_faces = UgGridHelpers::numFaces(grid);
const auto& fc = UgGridHelpers::faceCells(grid);
for (int face = 0; face < num_faces; ++face) {
const int c1 = fc(face, 0);
const int c2 = fc(face, 1);
if (c1 < 0 || c2 < 0) {
// Boundary face, skip this.
continue;
}
const int gc1 = (gc == 0) ? c1 : gc[c1];
const int gc2 = (gc == 0) ? c2 : gc[c2];
const int eq1 = eqlnumData[gc1];
const int eq2 = eqlnumData[gc2];
if (eq1 == eq2) {
// not an equilibration region boundary. skip this.
continue;
}
// update the maximum pressure potential difference between the two
// regions
const auto barrierId = std::make_pair(eq1, eq2);
if (maxDp.count(barrierId) == 0) {
maxDp[barrierId] = 0.0;
}
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const double z1 = UgGridHelpers::cellCenterDepth(grid, c1);
const double z2 = UgGridHelpers::cellCenterDepth(grid, c2);
const double zAvg = (z1 + z2)/2; // average depth
const double rhoAvg = (rho[phaseIdx][c1] + rho[phaseIdx][c2])/2;
const double s1 = initialState.saturation()[numPhases*c1 + phaseIdx];
const double s2 = initialState.saturation()[numPhases*c2 + phaseIdx];
const double sResid1 = minSat[numPhases*c1 + phaseIdx];
const double sResid2 = minSat[numPhases*c2 + phaseIdx];
// compute gravity corrected pressure potentials at the average depth
const double p1 = phasePressure[phaseIdx][c1] + rhoAvg*gravity*(zAvg - z1);
const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(zAvg - z2);
if ((p1 > p2 && s1 > sResid1) || (p2 > p1 && s2 > sResid2))
maxDp[barrierId] = std::max(maxDp[barrierId], std::abs(p1 - p2));
}
}
}
/// \brief Get a vector of pressure thresholds from EclipseState.
/// This function looks at EQLOPTS, THPRES and EQLNUM to determine
@ -37,6 +301,8 @@ namespace Opm
/// \tparam Grid Type of grid object (UnstructuredGrid or CpGrid).
/// \param[in] deck Input deck, EQLOPTS and THPRES are accessed from it.
/// \param[in] eclipseState Processed eclipse state, EQLNUM is accessed from it.
/// \param[in] maxDp The maximum gravity corrected pressure differences between
/// the equilibration regions.
/// \param[in] grid The grid to which the thresholds apply.
/// \return A vector of pressure thresholds, one
/// for each face in the grid. A value
@ -48,25 +314,28 @@ namespace Opm
template <class Grid>
std::vector<double> thresholdPressures(const ParseMode& parseMode ,EclipseStateConstPtr eclipseState, const Grid& grid)
std::vector<double> thresholdPressures(const DeckConstPtr& /* deck */,
EclipseStateConstPtr eclipseState,
const Grid& grid,
const std::map<std::pair<int, int>, double>& maxDp)
{
SimulationConfigConstPtr simulationConfig = eclipseState->getSimulationConfig();
std::vector<double> thpres_vals;
if (simulationConfig->hasThresholdPressure()) {
std::shared_ptr<const ThresholdPressure> thresholdPressure = simulationConfig->getThresholdPressure();
std::shared_ptr<GridProperty<int>> eqlnum = eclipseState->getIntGridProperty("EQLNUM");
auto eqlnumData = eqlnum->getData();
const auto& eqlnumData = eqlnum->getData();
// Set values for each face.
// Set threshold pressure values for each cell face.
const int num_faces = UgGridHelpers::numFaces(grid);
thpres_vals.resize(num_faces, 0.0);
const auto& fc = UgGridHelpers::faceCells(grid);
const int* gc = UgGridHelpers::globalCell(grid);
auto fc = UgGridHelpers::faceCells(grid);
thpres_vals.resize(num_faces, 0.0);
for (int face = 0; face < num_faces; ++face) {
const int c1 = fc(face, 0);
const int c2 = fc(face, 1);
if (c1 < 0 || c2 < 0) {
// Boundary face, skip this.
// Boundary face, skip it.
continue;
}
const int gc1 = (gc == 0) ? c1 : gc[c1];
@ -77,9 +346,13 @@ namespace Opm
if (thresholdPressure->hasRegionBarrier(eq1,eq2)) {
if (thresholdPressure->hasThresholdPressure(eq1,eq2)) {
thpres_vals[face] = thresholdPressure->getThresholdPressure(eq1,eq2);
} else {
std::string msg = "Initializing the THPRES pressure values from the initial state is not supported - using 0.0";
parseMode.handleError( ParseMode::UNSUPPORTED_INITIAL_THPRES , msg );
}
else {
// set the threshold pressure for faces of PVT regions where the third item
// has been defaulted to the maximum pressure potential difference between
// these regions
const auto barrierId = std::make_pair(eq1, eq2);
thpres_vals[face] = maxDp.at(barrierId);
}
}
}

17
opm/core/wells/WellsManager_impl.hpp
View File

@ -128,7 +128,11 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
wellperf_data.resize(wells.size());
wells_on_proc.resize(wells.size(), 1);
int well_index = 0;
// The well index on the current process.
// Note that some wells are deactivated as they live on the interior
// domain of another proccess. Therefore this might different from
// the index of the well according to the eclipse state
int well_index_on_proc = 0;
for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
WellConstPtr well = (*wellIter);
@ -198,7 +202,7 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
}
pd.well_index *= wellPi;
}
wellperf_data[well_index].push_back(pd);
wellperf_data[well_index_on_proc].push_back(pd);
}
} else {
++shut_completions_number;
@ -226,18 +230,19 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
// Check that the complete well is on this process
if ( sum_completions_on_proc < completionSet->size() )
{
std::cout<< "Well "<< well->name() << " semms not be in "
std::cout<< "Well "<< well->name() << " does not seem to be"
<< "completely in the disjoint partition of "
<< "process deactivating here." << std::endl;
<< "process. Therefore we deactivate it here." << std::endl;
// Mark well as not existent on this process
wells_on_proc[wellIter-wells.begin()] = 0;
wellperf_data[well_index_on_proc].clear();
continue;
}
}
}
}
{ // WELSPECS handling
well_names_to_index[well->name()] = well_index;
well_names_to_index[well->name()] = well_index_on_proc;
well_names.push_back(well->name());
{
WellData wd;
@ -253,7 +258,7 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
}
}
well_index++;
well_index_on_proc++;
}
// Set up reference depths that were defaulted. Count perfs.

2
tests/DuneIstlTestHelpers.hpp
View File

@ -61,7 +61,7 @@ struct MPIFixture {
};
BOOST_GLOBAL_FIXTURE(MPIFixture)
BOOST_GLOBAL_FIXTURE(MPIFixture);
struct MyMatrix
{

527
tests/norne_pvt.data Normal file
View File

@ -0,0 +1,527 @@
-- This reservoir simulation deck is made available under the Open Database
-- License: http://opendatacommons.org/licenses/odbl/1.0/. Any rights in
-- individual contents of the database are licensed under the Database Contents
-- License: http://opendatacommons.org/licenses/dbcl/1.0/
-- Copyright (C) 2015 Statoil
TABDIMS
--ntsfun ntpvt nssfun nppvt ntfip nrpvt ntendp
2 2 33 60 16 60 /
PVTG
-- PRESSURE RSG B-GAS VISCOSITY
-- BAR (CP)
50.00 0.00000497 0.024958 0.01441
0.00000248 0.024958 0.01440
0.00000000 0.024958 0.01440 /
70.00 0.00000521 0.017639 0.01491
0.00000261 0.017641 0.01490
0.00000000 0.017643 0.01490 /
90.00 0.00000627 0.013608 0.01547
0.00000313 0.013611 0.01546
0.00000000 0.013615 0.01544 /
110.00 0.00000798 0.011072 0.01609
0.00000399 0.011076 0.01607
0.00000000 0.011081 0.01605 /
130.00 0.00001041 0.009340 0.01677
0.00000520 0.009346 0.01674
0.00000000 0.009352 0.01671 /
150.00 0.00001365 0.008092 0.01752
0.00000683 0.008099 0.01748
0.00000000 0.008106 0.01743 /
170.00 0.00001786 0.007156 0.01834
0.00000893 0.007164 0.01827
0.00000000 0.007172 0.01819 /
190.00 0.00002316 0.006433 0.01923
0.00001158 0.006442 0.01912
0.00000000 0.006451 0.01900 /
210.00 0.00002972 0.005861 0.02019
0.00001486 0.005871 0.02001
0.00000000 0.005881 0.01984 /
230.00 0.00003767 0.005402 0.02121
0.00001883 0.005412 0.02095
0.00000000 0.005422 0.02071 /
250.80 0.00004756 0.005013 0.02234
0.00002378 0.005022 0.02197
0.00000000 0.005032 0.02162 /
268.42 0.00005757 0.004737 0.02335
0.00002878 0.004746 0.02287
0.00000000 0.004754 0.02240 /
285.33 0.00006853 0.004511 0.02438
0.00003427 0.004518 0.02375
0.00000000 0.004525 0.02315 /
301.59 0.00008041 0.004323 0.02542
0.00004020 0.004327 0.02463
0.00000000 0.004332 0.02387 /
317.23 0.00009313 0.004165 0.02648
0.00004657 0.004166 0.02549
0.00000000 0.004169 0.02456 /
332.29 0.00010668 0.004031 0.02755
0.00005334 0.004029 0.02634
0.00000000 0.004028 0.02522 /
346.80 0.00012100 0.003917 0.02863
0.00006050 0.003911 0.02719
0.00000000 0.003906 0.02585 /
360.80 0.00013607 0.003819 0.02974
0.00006803 0.003808 0.02803
0.00000000 0.003799 0.02645 /
374.31 0.00015188 0.003735 0.03087
0.00007594 0.003718 0.02887
0.00000000 0.003705 0.02703 /
387.36 0.00016843 0.003662 0.03202
0.00008421 0.003639 0.02970
0.00000000 0.003621 0.02758 /
399.99 0.00018571 0.003598 0.03320
0.00009286 0.003570 0.03053
0.00000000 0.003545 0.02810 /
412.21 0.00020375 0.003543 0.03442
0.00010188 0.003508 0.03137
0.00000000 0.003477 0.02861 /
424.05 0.00022256 0.003496 0.03566
0.00011128 0.003453 0.03220
0.00000000 0.003416 0.02909 /
435.53 0.00024218 0.003454 0.03695
0.00012109 0.003404 0.03305
0.00000000 0.003360 0.02956 /
446.68 0.00026266 0.003419 0.03828
0.00013133 0.003360 0.03390
0.00000000 0.003309 0.03000 /
457.51 0.00028404 0.003388 0.03967
0.00014202 0.003320 0.03477
0.00000000 0.003262 0.03043 /
468.04 0.00030639 0.003362 0.04110
0.00015319 0.003285 0.03566
0.00000000 0.003218 0.03085 /
478.30 0.00032980 0.003341 0.04261
0.00016490 0.003253 0.03656
0.00000000 0.003178 0.03125 /
488.30 0.00035436 0.003323 0.04418
0.00017718 0.003225 0.03749
0.00000000 0.003140 0.03164 /
498.06 0.00038020 0.003310 0.04583
0.00019010 0.003200 0.03845
0.00000000 0.003105 0.03202 /
507.59 0.00040745 0.003300 0.04758
0.00020373 0.003178 0.03944
0.00000000 0.003073 0.03238 /
516.92 0.00043630 0.003293 0.04943
0.00021815 0.003158 0.04048
0.00000000 0.003042 0.03273 /
526.06 0.00046694 0.003290 0.05141
0.00023347 0.003141 0.04156
0.00000000 0.003013 0.03308 /
535.02 0.00049963 0.003291 0.05353
0.00024981 0.003126 0.04271
0.00000000 0.002986 0.03342 /
543.83 0.00053469 0.003295 0.05582
0.00026734 0.003114 0.04393
0.00000000 0.002960 0.03374 /
552.49 0.00057251 0.003303 0.05830
0.00028625 0.003105 0.04523
0.00000000 0.002935 0.03407 /
561.04 0.00061359 0.003315 0.06103
0.00030679 0.003098 0.04664
0.00000000 0.002912 0.03438 /
569.48 0.00065855 0.003332 0.06405
0.00032928 0.003093 0.04818
0.00000000 0.002890 0.03469 /
577.82 0.00070820 0.003354 0.06744
0.00035410 0.003092 0.04988
0.00000000 0.002868 0.03500 /
586.09 0.00076355 0.003382 0.07127
0.00038178 0.003094 0.05178
0.00000000 0.002847 0.03530 /
594.29 0.00082592 0.003418 0.07567
0.00041296 0.003099 0.05394
0.00000000 0.002828 0.03560 /
/
-- PVT region 2 --
80.00 0.00000485 0.015151 0.01506
0.00000242 0.015154 0.01505
0.00000000 0.015157 0.01504 /
100.00 0.00000621 0.012032 0.01566
0.00000310 0.012036 0.01564
0.00000000 0.012040 0.01563 /
120.00 0.00000821 0.009980 0.01632
0.00000411 0.009985 0.01630
0.00000000 0.009990 0.01628 /
140.00 0.00001096 0.008537 0.01706
0.00000548 0.008544 0.01702
0.00000000 0.008550 0.01698 /
160.00 0.00001457 0.007476 0.01786
0.00000728 0.007484 0.01780
0.00000000 0.007492 0.01774 /
180.00 0.00001918 0.006669 0.01873
0.00000959 0.006678 0.01864
0.00000000 0.006687 0.01854 /
200.00 0.00002493 0.006038 0.01967
0.00001247 0.006048 0.01952
0.00000000 0.006058 0.01939 /
216.50 0.00003061 0.005616 0.02049
0.00001530 0.005626 0.02029
0.00000000 0.005636 0.02010 /
/
PVTO
-- RSO PRESSURE B-OIL VISCOSITY
-- (BAR) (CP)
20.59 50.00 1.10615 1.180
75.00 1.10164 1.247
100.00 1.09744 1.315
125.00 1.09351 1.384
150.00 1.08984 1.453 /
28.19 70.00 1.12522 1.066
95.00 1.12047 1.124
120.00 1.11604 1.182
145.00 1.11191 1.241
170.00 1.10804 1.300 /
36.01 90.00 1.14458 0.964
115.00 1.13959 1.014
140.00 1.13494 1.064
165.00 1.13060 1.115
190.00 1.12653 1.166 /
44.09 110.00 1.16437 0.880
135.00 1.15915 0.924
160.00 1.15428 0.968
185.00 1.14973 1.012
210.00 1.14547 1.056 /
52.46 130.00 1.18467 0.805
155.00 1.17921 0.843
180.00 1.17413 0.882
205.00 1.16937 0.920
230.00 1.16491 0.959 /
61.13 150.00 1.20555 0.746
175.00 1.19985 0.780
200.00 1.19454 0.814
225.00 1.18958 0.849
250.00 1.18492 0.883 /
70.14 170.00 1.22704 0.698
195.00 1.22111 0.729
220.00 1.21558 0.759
245.00 1.21040 0.790
270.00 1.20555 0.821 /
79.50 190.00 1.24922 0.658
215.00 1.24305 0.686
240.00 1.23729 0.714
265.00 1.23190 0.742
290.00 1.22685 0.770 /
89.24 210.00 1.27214 0.637
235.00 1.26573 0.664
260.00 1.25974 0.693
285.00 1.25414 0.725
310.00 1.24888 0.760 /
99.39 230.00 1.29586 0.622
255.00 1.28921 0.641
280.00 1.28300 0.661
305.00 1.27718 0.680
330.00 1.27171 0.699 /
110.41 250.80 1.32148 0.610
275.80 1.31457 0.628
300.80 1.30812 0.647
325.80 1.30207 0.665
350.80 1.29638 0.682 /
120.32 268.42 1.34449 0.576
293.42 1.33735 0.593
318.42 1.33068 0.609
343.42 1.32442 0.626
368.42 1.31853 0.642 /
130.23 285.33 1.36737 0.5335
310.33 1.36001 0.5487
335.33 1.35313 0.5638
360.33 1.34667 0.5787
385.33 1.34059 0.5934 /
140.12 301.59 1.39015 0.4956
326.59 1.38257 0.5094
351.59 1.37548 0.5230
376.59 1.36882 0.5365
401.59 1.36255 0.5498 /
150.01 317.23 1.41282 0.4614
342.23 1.40503 0.4739
367.23 1.39773 0.4863
392.23 1.39088 0.4986
417.23 1.38443 0.5107 /
159.89 332.29 1.43539 0.43042
357.29 1.42739 0.44183
382.29 1.41990 0.45312
407.29 1.41286 0.46430
432.29 1.40622 0.47537 /
169.76 346.80 1.45788 0.41191
371.80 1.44967 0.42260
396.80 1.44198 0.43318
421.80 1.43475 0.44365
446.80 1.42794 0.45402 /
179.63 360.80 1.48028 0.39503
385.80 1.47187 0.40508
410.80 1.46398 0.41502
435.80 1.45657 0.42487
460.80 1.44958 0.43461 /
189.48 374.31 1.50260 0.37959
399.31 1.49399 0.38907
424.31 1.48591 0.39845
449.31 1.47832 0.40773
474.31 1.47116 0.41692 /
199.34 387.36 1.52484 0.36543
412.36 1.51603 0.37439
437.36 1.50777 0.38326
462.36 1.50000 0.39203
487.36 1.49267 0.40072 /
209.18 399.99 1.54700 0.35239
424.99 1.53800 0.36089
449.99 1.52956 0.36929
474.99 1.52161 0.37762
499.99 1.51411 0.38585 /
219.02 412.21 1.56910 0.34035
437.21 1.55991 0.34843
462.21 1.55128 0.35642
487.21 1.54316 0.36433
512.21 1.53549 0.37216 /
228.85 424.05 1.59112 0.32921
449.05 1.58174 0.33691
474.05 1.57294 0.34453
499.05 1.56464 0.35206
524.05 1.55681 0.35952 /
238.67 435.53 1.61307 0.31888
460.53 1.60351 0.32623
485.53 1.59453 0.33350
510.53 1.58606 0.34070
535.53 1.57807 0.34782 /
248.48 446.68 1.63496 0.30927
471.68 1.62522 0.31630
496.68 1.61606 0.32326
521.68 1.60743 0.33014
546.68 1.59927 0.33695 /
258.29 457.51 1.65678 0.30032
482.51 1.64686 0.30706
507.51 1.63753 0.31373
532.51 1.62873 0.32032
557.51 1.62042 0.32685 /
268.09 468.04 1.67853 0.29196
493.04 1.66843 0.29843
518.04 1.65893 0.30483
543.04 1.64997 0.31117
568.04 1.64150 0.31743 /
277.89 478.30 1.70022 0.28414
503.30 1.68994 0.29037
528.30 1.68028 0.29652
553.30 1.67116 0.30261
578.30 1.66253 0.30864 /
287.68 488.30 1.72184 0.27681
513.30 1.71139 0.28281
538.30 1.70156 0.28874
563.30 1.69228 0.29460
588.30 1.68350 0.30040 /
297.46 498.06 1.74339 0.26994
523.06 1.73277 0.27572
548.06 1.72278 0.28144
573.06 1.71334 0.28709
598.06 1.70442 0.29269 /
307.23 507.59 1.76487 0.26347
532.59 1.75409 0.26906
557.59 1.74393 0.27458
582.59 1.73434 0.28004
607.59 1.72527 0.28544 /
317.00 516.92 1.78628 0.25738
541.92 1.77533 0.26279
566.92 1.76502 0.26812
591.92 1.75528 0.27340
616.92 1.74606 0.27863 /
326.76 526.06 1.80761 0.25165
551.06 1.79651 0.25688
576.06 1.78604 0.26204
601.06 1.77615 0.26716
626.06 1.76679 0.27221 /
336.51 535.02 1.82887 0.24623
560.02 1.81761 0.25130
585.02 1.80699 0.25631
610.02 1.79696 0.26126
635.02 1.78746 0.26616 /
346.26 543.83 1.85005 0.24112
568.83 1.83864 0.24603
593.83 1.82787 0.25089
618.83 1.81770 0.25570
643.83 1.80806 0.26045 /
356.00 552.49 1.87115 0.23628
577.49 1.85959 0.24105
602.49 1.84868 0.24577
627.49 1.83836 0.25043
652.49 1.82858 0.25505 /
365.73 561.04 1.89217 0.23170
586.04 1.88046 0.23634
611.04 1.86940 0.24092
636.04 1.85895 0.24546
661.04 1.84904 0.24994 /
375.46 569.48 1.91309 0.22736
594.48 1.90124 0.23187
619.48 1.89004 0.23633
644.48 1.87946 0.24074
669.48 1.86942 0.24510 /
385.18 577.82 1.93391 0.22325
602.82 1.92192 0.22764
627.82 1.91060 0.23198
652.82 1.89988 0.23627
677.82 1.88971 0.24052 /
394.89 586.09 1.95464 0.21934
611.09 1.94252 0.22362
636.09 1.93106 0.22785
661.09 1.92021 0.23204
686.09 1.90993 0.23617 /
404.60 594.29 1.97527 0.21564
619.29 1.96301 0.21981
644.29 1.95143 0.22393
669.29 1.94046 0.22801
694.29 1.93005 0.23204 /
/
-- PVT region 2
32.91 80.00 1.13304 1.04537
114.00 1.12837 1.10009
148.00 1.12401 1.15521
182.00 1.11994 1.21063 /
40.99 100.00 1.15276 0.97219
134.00 1.14786 1.02086
168.00 1.14328 1.06983
202.00 1.13900 1.11901 /
49.36 120.00 1.17297 0.91124
154.00 1.16783 0.95496
188.00 1.16303 0.99891
222.00 1.15854 1.04301 /
58.04 140.00 1.19374 0.85942
174.00 1.18836 0.89902
208.00 1.18334 0.93878
242.00 1.17864 0.97864 /
67.04 160.00 1.21512 0.81456
194.00 1.20951 0.85065
228.00 1.20426 0.88686
262.00 1.19935 0.92313 /
76.39 180.00 1.23718 0.77508
214.00 1.23132 0.80815
248.00 1.22585 0.84130
282.00 1.22073 0.87448 /
86.11 200.00 1.25996 0.73928
234.00 1.25386 0.76989
268.00 1.24816 0.80050
302.00 1.24283 0.83108 /
94.44 216.50 1.27934 0.67686
250.50 1.27304 0.70995
284.50 1.26716 0.74427
318.50 1.26164 0.77857 /
/

301
tests/test_norne_pvt.cpp Normal file
View File

@ -0,0 +1,301 @@
/*
Copyright 2015 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 NORNE_PVT_TESTS
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/test/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <sstream>
#include <iostream>
#include <opm/parser/eclipse/Units/ConversionFactors.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseMode.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/core/props/pvt/PvtLiveOil.hpp>
using namespace Opm;
/*
This file contains a regression test of the LiveOilPvt class for the
Norne data. The test is created by calling the PvtLiveOil::mu( )
function for a list of P and Rs values and storing the results in
the XX_expected vectors. No actual validation of the data has been
done, so the test will only serve to further cement possible bugs.
The columns of input data and expected data is organized in two
columns, that is mainly because two (P,Rs) values were selected from
each of the 'inner tables' in the PVTO keyword, and carries no
further semantic meaning.
*/
void check_vectors( const std::vector<double>& v1 , const std::vector<double>& v2) {
double tol = 1e-5;
for (decltype(v1.size()) i = 0; i < v1.size(); i++) {
BOOST_CHECK_CLOSE( v1[i] , v2[i] , tol );
}
}
void verify_norne_oil_pvt_region2(const TableManager& tableManager) {
auto pvtoTables = tableManager.getPvtoTables();
PvtLiveOil oilPvt(pvtoTables);
std::vector<double> rs = {33, 33,
43, 43,
53, 53,
61, 61,
70, 70,
80, 80,
100, 100, 100};
std::vector<double> P = {114, 148,
134, 168,
154, 188,
174, 208,
194, 228,
214, 248,
234, 268,
270 , 300, 330};
std::vector<double> mu_expected = {0.00106736588, 0.00113961037,
0.00093801366, 0.00099871729,
0.00083529743, 0.00088728769,
0.00077986989, 0.00082627508,
0.00072883113, 0.00076988665,
0.00068250424, 0.00072040786,
0.00062347677, 0.00064963306,
0.00065122911, 0.00409946846, 0.00472871311};
std::vector<double> b_expected = {0.88421444595, 0.88893909117,
0.86493342861, 0.86978957420,
0.84676402016, 0.85171762998,
0.83354279748, 0.83851861429,
0.81904041272, 0.82404719615,
0.80341044483, 0.80845950744,
0.77131381726, 0.77661604334,
0.77691738473, 0.98268067196, 0.98617572121};
{
std::vector<int> tableIndex(P.size() , 0);
std::vector<double> mu(P.size());
std::vector<double> dmudp(P.size());
std::vector<double> dmudr(P.size());
std::vector<double> b(P.size());
std::vector<double> dbdp(P.size());
std::vector<double> dbdr(P.size());
for (auto& value : P)
value *= Metric::Pressure;
for (auto& value : rs)
value *= Metric::GasDissolutionFactor;
oilPvt.mu( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , mu.data() , dmudp.data() , dmudr.data());
oilPvt.b( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , b.data() , dbdp.data() , dbdr.data());
check_vectors( mu , mu_expected );
check_vectors( b , b_expected );
}
}
void verify_norne_oil_pvt_region1(const TableManager& tableManager) {
auto pvtoTables = tableManager.getPvtoTables();
PvtLiveOil oilPvt(pvtoTables);
std::vector<double> rs = {21 , 21,
30 , 30,
38 , 38,
48 , 48,
55 , 55,
65 , 65,
75 , 75,
85 , 85,
95 , 95,
105 , 105,
115 , 115,
125 , 125,
135 , 135,
145 , 145,
155 , 155,
165 , 165,
175 , 175,
185 , 185,
195 , 195,
205 , 205,
215 , 215,
225 , 225,
234 , 234,
240 , 240,
252 , 252,
262 , 262,
272 , 272,
280 , 280,
410, 410, 410};
std::vector<double> P = {70, 110,
95, 145,
115, 165,
135, 185,
155, 205,
195, 245,
215, 265,
235, 285,
255, 305,
275, 325,
293, 343,
310, 360,
326, 376,
342, 392,
357, 407,
371, 420,
385, 435,
399, 450,
420, 480,
437, 487,
449, 499,
460, 510,
471, 521,
482, 532,
503, 553,
650, 680, 710};
std::vector<double> mu_expected = {0.00120767750, 0.00129077352,
0.00111063039, 0.00119627038,
0.00103118116, 0.00110633521,
0.00094413471, 0.00100998373,
0.00090320931, 0.00096374536,
0.00086714481, 0.00092142974,
0.00081811098, 0.00086735227,
0.00077704364, 0.00082229010,
0.00070975205, 0.00076029164,
0.00065679329, 0.00071124175,
0.00061496175, 0.00067213642,
0.00058000381, 0.00064115346,
0.00055124739, 0.00061633274,
0.00052840888, 0.00059781928,
0.00050926184, 0.00058323394,
0.00049295739, 0.00056996321,
0.00048026810, 0.00056474486,
0.00047088998, 0.00056427878,
0.00047649659, 0.00060774836,
0.00048006188, 0.00059909192,
0.00026623648, 0.00060915386,
0.00025670489, 0.00062157315,
0.00024760210, 0.00064290735,
0.00023889979, 0.00067946283,
0.00022330662, 0.00077837223,
0.01142273040, -0.00351292519, -0.00129867195};
std::vector<double> b_expected = {0.90699449462, 0.91120449633,
0.89040695696, 0.89551008140,
0.87548859167, 0.88062965205,
0.85697013389, 0.86224235632,
0.84533618728, 0.85061301709,
0.83069819286, 0.83585867335,
0.81473536808, 0.81994107210,
0.79955491390, 0.80479144821,
0.78507711370, 0.79032915313,
0.77073097762, 0.77596189361,
0.75627401890, 0.76141290296,
0.74161331648, 0.74678198081,
0.72686889575, 0.73206734035,
0.71214353439, 0.71737175926,
0.69733207231, 0.70259007745,
0.68243272267, 0.68761475238,
0.66755004999, 0.67286761567,
0.65268405426, 0.65813834713,
0.63858753316, 0.64504008462,
0.62408347496, 0.62949038145,
0.61223874629, 0.61449268543,
0.60422344638, 0.59939995459,
0.59620814647, 0.58594855211,
0.58819284656, 0.57739165219,
0.57289091037, 0.56019050084,
0.55474601877, 0.55809201119, 0.54526832277};
{
std::vector<int> tableIndex(P.size() , 1);
std::vector<double> mu(P.size());
std::vector<double> dmudp(P.size());
std::vector<double> dmudr(P.size());
std::vector<double> b(P.size());
std::vector<double> dbdp(P.size());
std::vector<double> dbdr(P.size());
for (auto& value : P)
value *= Metric::Pressure;
for (auto& value : rs)
value *= Metric::GasDissolutionFactor;
oilPvt.mu( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , mu.data() , dmudp.data() , dmudr.data());
oilPvt.b( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , b.data() , dbdp.data() , dbdr.data());
check_vectors( mu , mu_expected );
check_vectors( b , b_expected );
}
}
TableManager loadTables( const std::string& deck_file) {
Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
Opm::ParserPtr parser(new Parser());
std::shared_ptr<const Deck> deck;
deck = parser->parseFile(deck_file, parseMode);
return TableManager(*deck);
}
BOOST_AUTO_TEST_CASE( Test_Norne_PVT) {
TableManager tableManager = loadTables( "norne_pvt.data" );
verify_norne_oil_pvt_region1( tableManager );
verify_norne_oil_pvt_region2( tableManager );
}