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

Merge pull request #1630 from akva2/split_vfp

Browse Source 
Split VFP classes in ebos and legacy
This commit is contained in:
Atgeirr Flø Rasmussen
2018-11-14 11:04:05 +01:00
committed by GitHub
parent e1e39b4bd9 d3fa23bc2e
commit 49a7773b30
23 changed files with 1427 additions and 1128 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists_files.cmake
View File

@@ -41,9 +41,10 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/WellDensitySegmented.cpp
opm/autodiff/LinearisedBlackoilResidual.cpp
opm/autodiff/MPIUtilities.cpp
opm/autodiff/VFPProperties.cpp
opm/autodiff/VFPProdProperties.cpp
opm/autodiff/VFPProdPropertiesLegacy.cpp
opm/autodiff/VFPInjProperties.cpp
opm/autodiff/VFPInjPropertiesLegacy.cpp
opm/autodiff/MissingFeatures.cpp
opm/core/flowdiagnostics/AnisotropicEikonal.cpp
opm/core/flowdiagnostics/DGBasis.cpp
@@ -142,6 +143,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_transmissibilitymultipliers.cpp
tests/test_welldensitysegmented.cpp
tests/test_vfpproperties.cpp
tests/test_vfpproperties_legacy.cpp
tests/test_singlecellsolves.cpp
tests/test_milu.cpp
tests/test_multmatrixtransposed.cpp

5
opm/autodiff/BlackoilModelBase.hpp
View File

@@ -31,6 +31,8 @@
#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/autodiff/BlackoilModelEnums.hpp>
#include <opm/autodiff/VFPInjPropertiesLegacy.hpp>
#include <opm/autodiff/VFPProdPropertiesLegacy.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/IterationReport.hpp>
@@ -51,7 +53,6 @@ namespace Opm {
class DerivedGeology;
class RockCompressibility;
class NewtonIterationBlackoilInterface;
class VFPProperties;
/// Traits to encapsulate the types used by classes using or
/// extending this model. Forward declared here, must be
@@ -324,7 +325,7 @@ namespace Opm {
const BlackoilPropsAdFromDeck& fluid_;
const DerivedGeology& geo_;
const RockCompressibility* rock_comp_props_;
VFPProperties vfp_properties_;
VFPProperties<VFPInjPropertiesLegacy,VFPProdPropertiesLegacy> vfp_properties_;
const NewtonIterationBlackoilInterface& linsolver_;
// For each canonical phase -> true if active
const std::vector<bool> active_;

4
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@@ -36,8 +36,8 @@
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/autodiff/VFPProdPropertiesLegacy.hpp>
#include <opm/autodiff/VFPInjPropertiesLegacy.hpp>
#include <opm/grid/UnstructuredGrid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>

6
opm/autodiff/BlackoilWellModel.hpp
View File

@@ -40,11 +40,11 @@
#include <opm/core/wells/DynamicListEconLimited.hpp>
#include <opm/core/wells/WellCollection.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/WellHelpers.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/autodiff/BlackoilDetails.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/WellInterface.hpp>
@@ -302,7 +302,7 @@ namespace Opm {
bool initial_step_;
std::unique_ptr<RateConverterType> rateConverter_;
std::unique_ptr<VFPProperties> vfp_properties_;
std::unique_ptr<VFPProperties<VFPInjProperties,VFPProdProperties>> vfp_properties_;
SimulatorReport last_report_;

2
opm/autodiff/BlackoilWellModel_impl.hpp
View File

@@ -196,7 +196,7 @@ namespace Opm {
computeRESV(timeStepIdx);
// update VFP properties
vfp_properties_.reset (new VFPProperties (
vfp_properties_.reset (new VFPProperties<VFPInjProperties,VFPProdProperties> (
schedule().getVFPInjTables(timeStepIdx),
schedule().getVFPProdTables(timeStepIdx)) );

4
opm/autodiff/StandardWells.hpp
View File

@@ -76,7 +76,7 @@ namespace Opm {
void init(const BlackoilPropsAdFromDeck* fluid_arg,
const std::vector<bool>* active_arg,
const std::vector<PhasePresence>* pc_arg,
const VFPProperties* vfp_properties_arg,
const VFPProperties<VFPInjPropertiesLegacy,VFPProdPropertiesLegacy>* vfp_properties_arg,
const double gravity_arg,
const Vector& depth_arg);
@@ -223,7 +223,7 @@ namespace Opm {
const BlackoilPropsAdFromDeck* fluid_;
const std::vector<bool>* active_;
const std::vector<PhasePresence>* phase_condition_;
const VFPProperties* vfp_properties_;
const VFPProperties<VFPInjPropertiesLegacy,VFPProdPropertiesLegacy>* vfp_properties_;
double gravity_;
// the depth of the all the cell centers
// for standard Wells, it the same with the perforation depth

6
opm/autodiff/StandardWells_impl.hpp
View File

@@ -22,8 +22,8 @@
#include <opm/autodiff/StandardWells.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/VFPInjPropertiesLegacy.hpp>
#include <opm/autodiff/VFPProdPropertiesLegacy.hpp>
#include <opm/autodiff/WellHelpers.hpp>
@@ -96,7 +96,7 @@ namespace Opm
StandardWells::init(const BlackoilPropsAdFromDeck* fluid_arg,
const std::vector<bool>* active_arg,
const std::vector<PhasePresence>* pc_arg,
const VFPProperties* vfp_properties_arg,
const VFPProperties<VFPInjPropertiesLegacy,VFPProdPropertiesLegacy>* vfp_properties_arg,
const double gravity_arg,
const Vector& depth_arg)
{

215
opm/autodiff/VFPHelpers.hpp
View File

@@ -22,9 +22,10 @@
#define OPM_AUTODIFF_VFPHELPERS_HPP_
#include <cmath>
#include <opm/common/ErrorMacros.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/VFPInjTable.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
@@ -35,9 +36,6 @@ namespace Opm {
namespace detail {
typedef AutoDiffBlock<double> ADB;
/**
* Returns zero if input value is NaN
*/
@@ -45,6 +43,7 @@ inline double zeroIfNan(const double& value) {
return (std::isnan(value)) ? 0.0 : value;
}
/**
* Returns zero if input value is NaN
*/
@@ -54,24 +53,6 @@ inline EvalWell zeroIfNan(const EvalWell& value) {
}
/**
* Returns zero for every entry in the ADB which is NaN
*/
inline ADB zeroIfNan(const ADB& values) {
Selector<ADB::V::Scalar> not_nan_selector(values.value(), Selector<ADB::V::Scalar>::NotNaN);
const ADB::V z = ADB::V::Zero(values.size());
const ADB zero = ADB::constant(z, values.blockPattern());
ADB retval = not_nan_selector.select(values, zero);
return retval;
}
/**
* Computes the flo parameter according to the flo_type_
* for production tables
@@ -564,77 +545,6 @@ const T* getTable(const std::map<int, T*> tables, int table_id) {
}
/**
* Sets block_pattern to be the "union of x.blockPattern() and block_pattern".
*/
inline void extendBlockPattern(const ADB& x, std::vector<int>& block_pattern) {
std::vector<int> x_block_pattern = x.blockPattern();
if (x_block_pattern.empty()) {
return;
}
else {
if (block_pattern.empty()) {
block_pattern = x_block_pattern;
return;
}
else {
if (x_block_pattern != block_pattern) {
OPM_THROW(std::logic_error, "Block patterns do not match");
}
}
}
}
/**
* Finds the common block pattern for all inputs
*/
inline std::vector<int> commonBlockPattern(
const ADB& x1,
const ADB& x2,
const ADB& x3,
const ADB& x4) {
std::vector<int> block_pattern;
extendBlockPattern(x1, block_pattern);
extendBlockPattern(x2, block_pattern);
extendBlockPattern(x3, block_pattern);
extendBlockPattern(x4, block_pattern);
return block_pattern;
}
inline std::vector<int> commonBlockPattern(
const ADB& x1,
const ADB& x2,
const ADB& x3,
const ADB& x4,
const ADB& x5) {
std::vector<int> block_pattern = commonBlockPattern(x1, x2, x3, x4);
extendBlockPattern(x5, block_pattern);
return block_pattern;
}
/**
* Returns the type variable for FLO/GFR/WFR for production tables
*/
@@ -670,125 +580,6 @@ VFPInjTable::FLO_TYPE getType(const VFPInjTable* table) {
}
/**
* Returns the actual ADB for the type of FLO/GFR/WFR type
*/
template <typename TYPE>
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour, TYPE type);
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPProdTable::FLO_TYPE type) {
return detail::getFlo(aqua, liquid, vapour, type);
}
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPProdTable::WFR_TYPE type) {
return detail::getWFR(aqua, liquid, vapour, type);
}
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPProdTable::GFR_TYPE type) {
return detail::getGFR(aqua, liquid, vapour, type);
}
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPInjTable::FLO_TYPE type) {
return detail::getFlo(aqua, liquid, vapour, type);
}
/**
* Given m wells and n types of VFP variables (e.g., FLO = {FLO_OIL, FLO_LIQ}
* this function combines the n types of ADB objects, so that each of the
* m wells gets the right ADB.
* @param TYPE Type of variable to return, e.g., FLO_TYPE, WFR_TYPE, GFR_TYPE
* @param TABLE Type of table to use, e.g., VFPInjTable, VFPProdTable.
*/
template <typename TYPE, typename TABLE>
ADB combineADBVars(const std::vector<const TABLE*>& well_tables,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour) {
const int num_wells = static_cast<int>(well_tables.size());
assert(aqua.size() == num_wells);
assert(liquid.size() == num_wells);
assert(vapour.size() == num_wells);
//Caching variable for flo/wfr/gfr
std::map<TYPE, ADB> map;
//Indexing variable used when combining the different ADB types
std::map<TYPE, std::vector<int> > elems;
//Compute all of the different ADB types,
//and record which wells use which types
for (int i=0; i<num_wells; ++i) {
const TABLE* table = well_tables[i];
//Only do something if this well is under THP control
if (table != NULL) {
TYPE type = getType<TYPE>(table);
//"Caching" of flo_type etc: Only calculate used types
//Create type if it does not exist
if (map.find(type) == map.end()) {
map.insert(std::pair<TYPE, ADB>(
type,
detail::getValue<TYPE>(aqua, liquid, vapour, type)
));
}
//Add the index for assembly later in gather_vars
elems[type].push_back(i);
}
}
//Loop over all types of ADB variables, and combine them
//so that each well gets the proper variable
ADB retval = ADB::constant(ADB::V::Zero(num_wells));
for (const auto& entry : elems) {
const auto& key = entry.first;
const auto& value = entry.second;
//Get the ADB for this type of variable
assert(map.find(key) != map.end());
const ADB& values = map.find(key)->second;
//Get indices to all elements that should use this ADB
const std::vector<int>& current = value;
//Add these elements to retval
retval = retval + superset(subset(values, current), current, values.size());
}
return retval;
}
/**
* Helper function that finds x for a given value of y for a line
* *NOTE ORDER OF ARGUMENTS*

238
opm/autodiff/VFPHelpersLegacy.hpp Normal file
View File

@@ -0,0 +1,238 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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_AUTODIFF_VFPHELPERSLEGACY_HPP_
#define OPM_AUTODIFF_VFPHELPERSLEGACY_HPP_
#include <opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/VFPInjTable.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
/**
* This file contains a set of helper functions used by VFPProd / VFPInj.
*/
namespace Opm {
namespace detail {
typedef AutoDiffBlock<double> ADB;
/**
* Returns zero for every entry in the ADB which is NaN
*/
inline ADB zeroIfNan(const ADB& values) {
Selector<ADB::V::Scalar> not_nan_selector(values.value(), Selector<ADB::V::Scalar>::NotNaN);
const ADB::V z = ADB::V::Zero(values.size());
const ADB zero = ADB::constant(z, values.blockPattern());
ADB retval = not_nan_selector.select(values, zero);
return retval;
}
/**
* Sets block_pattern to be the "union of x.blockPattern() and block_pattern".
*/
inline void extendBlockPattern(const ADB& x, std::vector<int>& block_pattern) {
std::vector<int> x_block_pattern = x.blockPattern();
if (x_block_pattern.empty()) {
return;
}
else {
if (block_pattern.empty()) {
block_pattern = x_block_pattern;
return;
}
else {
if (x_block_pattern != block_pattern) {
OPM_THROW(std::logic_error, "Block patterns do not match");
}
}
}
}
}
}
#include "VFPHelpers.hpp"
namespace Opm {
namespace detail {
/**
* Finds the common block pattern for all inputs
*/
inline std::vector<int> commonBlockPattern(
const ADB& x1,
const ADB& x2,
const ADB& x3,
const ADB& x4) {
std::vector<int> block_pattern;
extendBlockPattern(x1, block_pattern);
extendBlockPattern(x2, block_pattern);
extendBlockPattern(x3, block_pattern);
extendBlockPattern(x4, block_pattern);
return block_pattern;
}
inline std::vector<int> commonBlockPattern(
const ADB& x1,
const ADB& x2,
const ADB& x3,
const ADB& x4,
const ADB& x5) {
std::vector<int> block_pattern = commonBlockPattern(x1, x2, x3, x4);
extendBlockPattern(x5, block_pattern);
return block_pattern;
}
/**
* Returns the actual ADB for the type of FLO/GFR/WFR type
*/
template <typename TYPE>
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour, TYPE type);
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPProdTable::FLO_TYPE type) {
return detail::getFlo(aqua, liquid, vapour, type);
}
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPProdTable::WFR_TYPE type) {
return detail::getWFR(aqua, liquid, vapour, type);
}
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPProdTable::GFR_TYPE type) {
return detail::getGFR(aqua, liquid, vapour, type);
}
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPInjTable::FLO_TYPE type) {
return detail::getFlo(aqua, liquid, vapour, type);
}
/**
* Given m wells and n types of VFP variables (e.g., FLO = {FLO_OIL, FLO_LIQ}
* this function combines the n types of ADB objects, so that each of the
* m wells gets the right ADB.
* @param TYPE Type of variable to return, e.g., FLO_TYPE, WFR_TYPE, GFR_TYPE
* @param TABLE Type of table to use, e.g., VFPInjTable, VFPProdTable.
*/
template <typename TYPE, typename TABLE>
ADB combineADBVars(const std::vector<const TABLE*>& well_tables,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour) {
const int num_wells = static_cast<int>(well_tables.size());
assert(aqua.size() == num_wells);
assert(liquid.size() == num_wells);
assert(vapour.size() == num_wells);
//Caching variable for flo/wfr/gfr
std::map<TYPE, ADB> map;
//Indexing variable used when combining the different ADB types
std::map<TYPE, std::vector<int> > elems;
//Compute all of the different ADB types,
//and record which wells use which types
for (int i=0; i<num_wells; ++i) {
const TABLE* table = well_tables[i];
//Only do something if this well is under THP control
if (table != NULL) {
TYPE type = getType<TYPE>(table);
//"Caching" of flo_type etc: Only calculate used types
//Create type if it does not exist
if (map.find(type) == map.end()) {
map.insert(std::pair<TYPE, ADB>(
type,
detail::getValue<TYPE>(aqua, liquid, vapour, type)
));
}
//Add the index for assembly later in gather_vars
elems[type].push_back(i);
}
}
//Loop over all types of ADB variables, and combine them
//so that each well gets the proper variable
ADB retval = ADB::constant(ADB::V::Zero(num_wells));
for (const auto& entry : elems) {
const auto& key = entry.first;
const auto& value = entry.second;
//Get the ADB for this type of variable
assert(map.find(key) != map.end());
const ADB& values = map.find(key)->second;
//Get indices to all elements that should use this ADB
const std::vector<int>& current = value;
//Add these elements to retval
retval = retval + superset(subset(values, current), current, values.size());
}
return retval;
}
} // namespace detail
} // namespace
#endif /* OPM_AUTODIFF_VFPHELPERSLEGACY_HPP_ */

146
opm/autodiff/VFPInjProperties.cpp
View File

@@ -26,146 +26,34 @@
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/VFPHelpers.hpp>
namespace Opm {
namespace detail {
} //Namespace detail
VFPInjProperties::VFPInjProperties() {
}
VFPInjProperties::VFPInjProperties(const VFPInjTable* table){
m_tables[table->getTableNum()] = table;
}
VFPInjProperties::VFPInjProperties(const std::map<int, std::shared_ptr<const VFPInjTable> >& tables) {
VFPInjProperties::VFPInjProperties(const VFPInjProperties::InjTable& tables) {
for (const auto& table : tables) {
m_tables[table.first] = table.second.get();
}
}
VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp_val) const {
const int nw = wells.number_of_wells;
//Short-hands for water / oil / gas phases
//TODO enable support for two-phase.
assert(wells.number_of_phases == 3);
const ADB& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
return bhp(table_id, w, o, g, thp_val);
}
VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp_arg) const {
const int nw = thp_arg.size();
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp_arg);
assert(static_cast<int>(table_id.size()) == nw);
assert(aqua.size() == nw);
assert(liquid.size() == nw);
assert(vapour.size() == nw);
assert(thp_arg.size() == nw);
//Allocate data for bhp's and partial derivatives
ADB::V value = ADB::V::Zero(nw);
ADB::V dthp = ADB::V::Zero(nw);
ADB::V dflo = ADB::V::Zero(nw);
//Get the table for each well
std::vector<const VFPInjTable*> well_tables(nw, nullptr);
for (int i=0; i<nw; ++i) {
if (table_id[i] > 0) {
well_tables[i] = detail::getTable(m_tables, table_id[i]);
}
}
//Get the right FLO variable for each well as a single ADB
const ADB flo = detail::combineADBVars<VFPInjTable::FLO_TYPE>(well_tables, aqua, liquid, vapour);
//Compute the BHP for each well independently
for (int i=0; i<nw; ++i) {
const VFPInjTable* table = well_tables[i];
if (table != nullptr) {
//First, find the values to interpolate between
auto flo_i = detail::findInterpData(flo.value()[i], table->getFloAxis());
auto thp_i = detail::findInterpData(thp_arg.value()[i], table->getTHPAxis());
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i);
value[i] = bhp_val.value;
dthp[i] = bhp_val.dthp;
dflo[i] = bhp_val.dflo;
}
else {
value[i] = -1e100; //Signal that this value has not been calculated properly, due to "missing" table
}
}
//Create diagonal matrices from ADB::Vs
ADB::M dthp_diag(dthp.matrix().asDiagonal());
ADB::M dflo_diag(dflo.matrix().asDiagonal());
//Calculate the Jacobians
const int num_blocks = block_pattern.size();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
//Could have used fastSparseProduct and temporary variables
//but may not save too much on that.
jacs[block] = ADB::M(nw, block_pattern[block]);
if (!thp_arg.derivative().empty()) {
jacs[block] += dthp_diag * thp_arg.derivative()[block];
}
if (!flo.derivative().empty()) {
jacs[block] += dflo_diag * flo.derivative()[block];
}
}
ADB retval = ADB::function(std::move(value), std::move(jacs));
return retval;
}
double VFPInjProperties::bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp_arg) const {
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp_arg) const {
const VFPInjTable* table = detail::getTable(m_tables, table_id);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg);
@@ -173,17 +61,11 @@ double VFPInjProperties::bhp(int table_id,
}
double VFPInjProperties::thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp_arg) const {
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp_arg) const {
const VFPInjTable* table = detail::getTable(m_tables, table_id);
const VFPInjTable::array_type& data = table->getTable();
@@ -210,18 +92,8 @@ double VFPInjProperties::thp(int table_id,
}
const VFPInjTable* VFPInjProperties::getTable(const int table_id) const {
return detail::getTable(m_tables, table_id);
}
} //Namespace Opm

73
opm/autodiff/VFPInjProperties.hpp
View File

@@ -33,13 +33,8 @@
namespace Opm {
template <class Scalar>
class AutoDiffBlock;
class VFPInjProperties {
public:
typedef AutoDiffBlock<double> ADB;
/**
* Empty constructor
*/
@@ -57,43 +52,8 @@ public:
* Takes *no* ownership of data.
* @param inj_tables A map of different VFPINJ tables.
*/
explicit VFPInjProperties(const std::map<int, std::shared_ptr<const VFPInjTable> >& inj_tables);
/**
* Linear interpolation of bhp as function of the input parameters.
* @param table_id Table number to use
* @param wells Wells structure with information about wells in qs
* @param qs Flow quantities
* @param thp Tubing head pressure
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
ADB bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp) const;
/**
* Linear interpolation of bhp as a function of the input parameters given as ADBs
* Each entry corresponds typically to one well.
* @param table_id Table number to use. A negative entry (e.g., -1)
* will indicate that no table is used, and the corresponding
* BHP will be calculated as a constant -1e100.
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param thp Tubing head pressure
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table, for each entry in the
* input ADB objects.
*/
ADB bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp) const;
using InjTable = std::map<int, std::shared_ptr<const VFPInjTable> >;
explicit VFPInjProperties(const InjTable& inj_tables);
/**
* Linear interpolation of bhp as a function of the input parameters given as
@@ -143,6 +103,19 @@ public:
return bhp;
}
/**
* Returns the table associated with the ID, or throws an exception if
* the table does not exist
*/
const VFPInjTable* getTable(const int table_id) const;
/**
* Returns true if no vfp tables are in the current map
*/
bool empty() const {
return m_tables.empty();
}
/**
* Linear interpolation of bhp as a function of the input parameters
* @param table_id Table number to use
@@ -160,7 +133,6 @@ public:
const double& vapour,
const double& thp) const;
/**
* Linear interpolation of thp as a function of the input parameters
* @param table_id Table number to use
@@ -178,20 +150,7 @@ public:
const double& vapour,
const double& bhp) const;
/**
* Returns the table associated with the ID, or throws an exception if
* the table does not exist
*/
const VFPInjTable* getTable(const int table_id) const;
/**
* Returns true if no vfp tables are in the current map
*/
bool empty() const {
return m_tables.empty();
}
private:
protected:
// Map which connects the table number with the table itself
std::map<int, const VFPInjTable*> m_tables;
};

130
opm/autodiff/VFPInjPropertiesLegacy.cpp Normal file
View File

@@ -0,0 +1,130 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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"
#include <opm/autodiff/VFPHelpersLegacy.hpp>
#include <opm/autodiff/VFPInjPropertiesLegacy.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
namespace Opm {
VFPInjPropertiesLegacy::ADB VFPInjPropertiesLegacy::bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp_val) const {
const int nw = wells.number_of_wells;
//Short-hands for water / oil / gas phases
//TODO enable support for two-phase.
assert(wells.number_of_phases == 3);
const ADB& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
return bhp(table_id, w, o, g, thp_val);
}
VFPInjPropertiesLegacy::ADB VFPInjPropertiesLegacy::bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp_arg) const {
const int nw = thp_arg.size();
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp_arg);
assert(static_cast<int>(table_id.size()) == nw);
assert(aqua.size() == nw);
assert(liquid.size() == nw);
assert(vapour.size() == nw);
assert(thp_arg.size() == nw);
//Allocate data for bhp's and partial derivatives
ADB::V value = ADB::V::Zero(nw);
ADB::V dthp = ADB::V::Zero(nw);
ADB::V dflo = ADB::V::Zero(nw);
//Get the table for each well
std::vector<const VFPInjTable*> well_tables(nw, nullptr);
for (int i=0; i<nw; ++i) {
if (table_id[i] > 0) {
well_tables[i] = detail::getTable(m_tables, table_id[i]);
}
}
//Get the right FLO variable for each well as a single ADB
const ADB flo = detail::combineADBVars<VFPInjTable::FLO_TYPE>(well_tables, aqua, liquid, vapour);
//Compute the BHP for each well independently
for (int i=0; i<nw; ++i) {
const VFPInjTable* table = well_tables[i];
if (table != nullptr) {
//First, find the values to interpolate between
auto flo_i = detail::findInterpData(flo.value()[i], table->getFloAxis());
auto thp_i = detail::findInterpData(thp_arg.value()[i], table->getTHPAxis());
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i);
value[i] = bhp_val.value;
dthp[i] = bhp_val.dthp;
dflo[i] = bhp_val.dflo;
}
else {
value[i] = -1e100; //Signal that this value has not been calculated properly, due to "missing" table
}
}
//Create diagonal matrices from ADB::Vs
ADB::M dthp_diag(dthp.matrix().asDiagonal());
ADB::M dflo_diag(dflo.matrix().asDiagonal());
//Calculate the Jacobians
const int num_blocks = block_pattern.size();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
//Could have used fastSparseProduct and temporary variables
//but may not save too much on that.
jacs[block] = ADB::M(nw, block_pattern[block]);
if (!thp_arg.derivative().empty()) {
jacs[block] += dthp_diag * thp_arg.derivative()[block];
}
if (!flo.derivative().empty()) {
jacs[block] += dflo_diag * flo.derivative()[block];
}
}
ADB retval = ADB::function(std::move(value), std::move(jacs));
return retval;
}
} //Namespace Opm

104
opm/autodiff/VFPInjPropertiesLegacy.hpp Normal file
View File

@@ -0,0 +1,104 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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_AUTODIFF_VFPINJPROPERTIESLEGACY_HPP_
#define OPM_AUTODIFF_VFPINJPROPERTIESLEGACY_HPP_
#include <opm/autodiff/VFPHelpersLegacy.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
#include <vector>
#include <map>
namespace Opm {
template <class Scalar>
class AutoDiffBlock;
class VFPInjPropertiesLegacy : public VFPInjProperties {
public:
typedef AutoDiffBlock<double> ADB;
/**
* Empty constructor
*/
VFPInjPropertiesLegacy() {}
/**
* Constructor
* Takes *no* ownership of data.
* @param inj_table A *single* VFPINJ table
*/
explicit VFPInjPropertiesLegacy(const VFPInjTable* inj_table) :
VFPInjProperties(inj_table) {}
/**
* Constructor
* Takes *no* ownership of data.
* @param inj_tables A map of different VFPINJ tables.
*/
explicit VFPInjPropertiesLegacy(const InjTable& inj_tables) :
VFPInjProperties(inj_tables) {}
using VFPInjProperties::bhp;
/**
* Linear interpolation of bhp as function of the input parameters.
* @param table_id Table number to use
* @param wells Wells structure with information about wells in qs
* @param qs Flow quantities
* @param thp Tubing head pressure
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
ADB bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp) const;
/**
* Linear interpolation of bhp as a function of the input parameters given as ADBs
* Each entry corresponds typically to one well.
* @param table_id Table number to use. A negative entry (e.g., -1)
* will indicate that no table is used, and the corresponding
* BHP will be calculated as a constant -1e100.
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param thp Tubing head pressure
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table, for each entry in the
* input ADB objects.
*/
ADB bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp) const;
};
} //namespace
#endif /* OPM_AUTODIFF_VFPINJPROPERTIES_HPP_ */

166
opm/autodiff/VFPProdProperties.cpp
View File

@@ -22,8 +22,6 @@
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/autodiff/VFPHelpers.hpp>
@@ -40,160 +38,24 @@ VFPProdProperties::VFPProdProperties() {
}
VFPProdProperties::VFPProdProperties(const VFPProdTable* table){
m_tables[table->getTableNum()] = table;
}
VFPProdProperties::VFPProdProperties(const std::map<int, std::shared_ptr<const VFPProdTable> >& tables) {
VFPProdProperties::VFPProdProperties(const VFPProdProperties::ProdTable& tables) {
for (const auto& table : tables) {
m_tables[table.first] = table.second.get();
}
}
VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp_arg,
const ADB& alq) const {
const int nw = wells.number_of_wells;
//Short-hands for water / oil / gas phases
//TODO enable support for two-phase.
assert(wells.number_of_phases == 3);
const ADB& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
return bhp(table_id, w, o, g, thp_arg, alq);
}
VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp_arg,
const ADB& alq) const {
const int nw = thp_arg.size();
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp_arg, alq);
assert(static_cast<int>(table_id.size()) == nw);
assert(aqua.size() == nw);
assert(liquid.size() == nw);
assert(vapour.size() == nw);
assert(thp_arg.size() == nw);
assert(alq.size() == nw);
//Allocate data for bhp's and partial derivatives
ADB::V value = ADB::V::Zero(nw);
ADB::V dthp = ADB::V::Zero(nw);
ADB::V dwfr = ADB::V::Zero(nw);
ADB::V dgfr = ADB::V::Zero(nw);
ADB::V dalq = ADB::V::Zero(nw);
ADB::V dflo = ADB::V::Zero(nw);
//Get the table for each well
std::vector<const VFPProdTable*> well_tables(nw, nullptr);
for (int i=0; i<nw; ++i) {
if (table_id[i] >= 0) {
well_tables[i] = detail::getTable(m_tables, table_id[i]);
}
}
//Get the right FLO/GFR/WFR variable for each well as a single ADB
const ADB flo = detail::combineADBVars<VFPProdTable::FLO_TYPE>(well_tables, aqua, liquid, vapour);
const ADB wfr = detail::combineADBVars<VFPProdTable::WFR_TYPE>(well_tables, aqua, liquid, vapour);
const ADB gfr = detail::combineADBVars<VFPProdTable::GFR_TYPE>(well_tables, aqua, liquid, vapour);
//Compute the BHP for each well independently
for (int i=0; i<nw; ++i) {
const VFPProdTable* table = well_tables[i];
if (table != nullptr) {
//First, find the values to interpolate between
//Value of FLO is negative in OPM for producers, but positive in VFP table
auto flo_i = detail::findInterpData(-flo.value()[i], table->getFloAxis());
auto thp_i = detail::findInterpData( thp_arg.value()[i], table->getTHPAxis());
auto wfr_i = detail::findInterpData( wfr.value()[i], table->getWFRAxis());
auto gfr_i = detail::findInterpData( gfr.value()[i], table->getGFRAxis());
auto alq_i = detail::findInterpData( alq.value()[i], table->getALQAxis());
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
value[i] = bhp_val.value;
dthp[i] = bhp_val.dthp;
dwfr[i] = bhp_val.dwfr;
dgfr[i] = bhp_val.dgfr;
dalq[i] = bhp_val.dalq;
dflo[i] = bhp_val.dflo;
}
else {
value[i] = -1e100; //Signal that this value has not been calculated properly, due to "missing" table
}
}
//Create diagonal matrices from ADB::Vs
ADB::M dthp_diag(dthp.matrix().asDiagonal());
ADB::M dwfr_diag(dwfr.matrix().asDiagonal());
ADB::M dgfr_diag(dgfr.matrix().asDiagonal());
ADB::M dalq_diag(dalq.matrix().asDiagonal());
ADB::M dflo_diag(dflo.matrix().asDiagonal());
//Calculate the Jacobians
const int num_blocks = block_pattern.size();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
//Could have used fastSparseProduct and temporary variables
//but may not save too much on that.
jacs[block] = ADB::M(nw, block_pattern[block]);
if (!thp_arg.derivative().empty()) {
jacs[block] += dthp_diag * thp_arg.derivative()[block];
}
if (!wfr.derivative().empty()) {
jacs[block] += dwfr_diag * wfr.derivative()[block];
}
if (!gfr.derivative().empty()) {
jacs[block] += dgfr_diag * gfr.derivative()[block];
}
if (!alq.derivative().empty()) {
jacs[block] += dalq_diag * alq.derivative()[block];
}
if (!flo.derivative().empty()) {
jacs[block] -= dflo_diag * flo.derivative()[block];
}
}
ADB retval = ADB::function(std::move(value), std::move(jacs));
return retval;
}
double VFPProdProperties::bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp_arg,
const double& alq) const {
const VFPProdTable* table = detail::getTable(m_tables, table_id);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq);
return retval.value;
}
double VFPProdProperties::thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp_arg,
const double& alq) const {
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp_arg,
const double& alq) const {
const VFPProdTable* table = detail::getTable(m_tables, table_id);
const VFPProdTable::array_type& data = table->getTable();
@@ -226,8 +88,17 @@ double VFPProdProperties::thp(int table_id,
}
double VFPProdProperties::bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp_arg,
const double& alq) const {
const VFPProdTable* table = detail::getTable(m_tables, table_id);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq);
return retval.value;
}
const VFPProdTable* VFPProdProperties::getTable(const int table_id) const {
@@ -235,9 +106,4 @@ const VFPProdTable* VFPProdProperties::getTable(const int table_id) const {
}
}

51
opm/autodiff/VFPProdProperties.hpp
View File

@@ -32,9 +32,6 @@
namespace Opm {
template <class Scalar>
class AutoDiffBlock;
/**
* Class which linearly interpolates BHP as a function of rate, tubing head pressure,
* water fraction, gas fraction, and artificial lift for production VFP tables, and similarly
@@ -42,8 +39,6 @@ class AutoDiffBlock;
*/
class VFPProdProperties {
public:
typedef AutoDiffBlock<double> ADB;
/**
* Empty constructor
*/
@@ -61,48 +56,8 @@ public:
* Takes *no* ownership of data.
* @param prod_tables A map of different VFPPROD tables.
*/
explicit VFPProdProperties(const std::map<int, std::shared_ptr<const VFPProdTable> >& prod_tables);
/**
* Linear interpolation of bhp as function of the input parameters.
* @param table_id Table number to use
* @param wells Wells structure with information about wells in qs
* @param qs Flow quantities
* @param thp Tubing head pressure
* @param alq Artificial lift or other parameter
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
ADB bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp,
const ADB& alq) const;
/**
* Linear interpolation of bhp as a function of the input parameters given as ADBs
* Each entry corresponds typically to one well.
* @param table_id Table number to use. A negative entry (e.g., -1)
* will indicate that no table is used, and the corresponding
* BHP will be calculated as a constant -1e100.
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param thp Tubing head pressure
* @param alq Artificial lift or other parameter
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table, for each entry in the
* input ADB objects.
*/
ADB bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp,
const ADB& alq) const;
using ProdTable = std::map<int, std::shared_ptr<const VFPProdTable> >;
explicit VFPProdProperties(const ProdTable& prod_tables);
/**
* Linear interpolation of bhp as a function of the input parameters given as
@@ -210,7 +165,7 @@ public:
return m_tables.empty();
}
private:
protected:
// Map which connects the table number with the table itself
std::map<int, const VFPProdTable*> m_tables;
};

154
opm/autodiff/VFPProdPropertiesLegacy.cpp Normal file
View File

@@ -0,0 +1,154 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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"
#include <opm/autodiff/VFPHelpersLegacy.hpp>
#include <opm/autodiff/VFPProdPropertiesLegacy.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
namespace Opm {
VFPProdPropertiesLegacy::ADB VFPProdPropertiesLegacy::bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp_arg,
const ADB& alq) const {
const int nw = wells.number_of_wells;
//Short-hands for water / oil / gas phases
//TODO enable support for two-phase.
assert(wells.number_of_phases == 3);
const ADB& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
return bhp(table_id, w, o, g, thp_arg, alq);
}
VFPProdPropertiesLegacy::ADB VFPProdPropertiesLegacy::bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp_arg,
const ADB& alq) const {
const int nw = thp_arg.size();
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp_arg, alq);
assert(static_cast<int>(table_id.size()) == nw);
assert(aqua.size() == nw);
assert(liquid.size() == nw);
assert(vapour.size() == nw);
assert(thp_arg.size() == nw);
assert(alq.size() == nw);
//Allocate data for bhp's and partial derivatives
ADB::V value = ADB::V::Zero(nw);
ADB::V dthp = ADB::V::Zero(nw);
ADB::V dwfr = ADB::V::Zero(nw);
ADB::V dgfr = ADB::V::Zero(nw);
ADB::V dalq = ADB::V::Zero(nw);
ADB::V dflo = ADB::V::Zero(nw);
//Get the table for each well
std::vector<const VFPProdTable*> well_tables(nw, nullptr);
for (int i=0; i<nw; ++i) {
if (table_id[i] >= 0) {
well_tables[i] = detail::getTable(m_tables, table_id[i]);
}
}
//Get the right FLO/GFR/WFR variable for each well as a single ADB
const ADB flo = detail::combineADBVars<VFPProdTable::FLO_TYPE>(well_tables, aqua, liquid, vapour);
const ADB wfr = detail::combineADBVars<VFPProdTable::WFR_TYPE>(well_tables, aqua, liquid, vapour);
const ADB gfr = detail::combineADBVars<VFPProdTable::GFR_TYPE>(well_tables, aqua, liquid, vapour);
//Compute the BHP for each well independently
for (int i=0; i<nw; ++i) {
const VFPProdTable* table = well_tables[i];
if (table != nullptr) {
//First, find the values to interpolate between
//Value of FLO is negative in OPM for producers, but positive in VFP table
auto flo_i = detail::findInterpData(-flo.value()[i], table->getFloAxis());
auto thp_i = detail::findInterpData( thp_arg.value()[i], table->getTHPAxis());
auto wfr_i = detail::findInterpData( wfr.value()[i], table->getWFRAxis());
auto gfr_i = detail::findInterpData( gfr.value()[i], table->getGFRAxis());
auto alq_i = detail::findInterpData( alq.value()[i], table->getALQAxis());
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
value[i] = bhp_val.value;
dthp[i] = bhp_val.dthp;
dwfr[i] = bhp_val.dwfr;
dgfr[i] = bhp_val.dgfr;
dalq[i] = bhp_val.dalq;
dflo[i] = bhp_val.dflo;
}
else {
value[i] = -1e100; //Signal that this value has not been calculated properly, due to "missing" table
}
}
//Create diagonal matrices from ADB::Vs
ADB::M dthp_diag(dthp.matrix().asDiagonal());
ADB::M dwfr_diag(dwfr.matrix().asDiagonal());
ADB::M dgfr_diag(dgfr.matrix().asDiagonal());
ADB::M dalq_diag(dalq.matrix().asDiagonal());
ADB::M dflo_diag(dflo.matrix().asDiagonal());
//Calculate the Jacobians
const int num_blocks = block_pattern.size();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
//Could have used fastSparseProduct and temporary variables
//but may not save too much on that.
jacs[block] = ADB::M(nw, block_pattern[block]);
if (!thp_arg.derivative().empty()) {
jacs[block] += dthp_diag * thp_arg.derivative()[block];
}
if (!wfr.derivative().empty()) {
jacs[block] += dwfr_diag * wfr.derivative()[block];
}
if (!gfr.derivative().empty()) {
jacs[block] += dgfr_diag * gfr.derivative()[block];
}
if (!alq.derivative().empty()) {
jacs[block] += dalq_diag * alq.derivative()[block];
}
if (!flo.derivative().empty()) {
jacs[block] -= dflo_diag * flo.derivative()[block];
}
}
ADB retval = ADB::function(std::move(value), std::move(jacs));
return retval;
}
}

110
opm/autodiff/VFPProdPropertiesLegacy.hpp Normal file
View File

@@ -0,0 +1,110 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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_AUTODIFF_VFPPRODPROPERTIESLEGACY_HPP_
#define OPM_AUTODIFF_VFPPRODPROPERTIESLEGACY_HPP_
#include <opm/autodiff/VFPHelpersLegacy.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <vector>
#include <map>
namespace Opm {
template <class Scalar>
class AutoDiffBlock;
/**
* Class which linearly interpolates BHP as a function of rate, tubing head pressure,
* water fraction, gas fraction, and artificial lift for production VFP tables, and similarly
* the BHP as a function of the rate and tubing head pressure.
*/
class VFPProdPropertiesLegacy : public VFPProdProperties {
public:
typedef AutoDiffBlock<double> ADB;
/**
* Empty constructor
*/
VFPProdPropertiesLegacy() {}
/**
* Constructor
* Takes *no* ownership of data.
* @param prod_table A *single* VFPPROD table
*/
explicit VFPProdPropertiesLegacy(const VFPProdTable* prod_table) :
VFPProdProperties(prod_table) {}
/**
* Constructor
* Takes *no* ownership of data.
* @param prod_tables A map of different VFPPROD tables.
*/
explicit VFPProdPropertiesLegacy(const ProdTable& prod_tables) :
VFPProdProperties(prod_tables) {}
using VFPProdProperties::bhp;
/**
* Linear interpolation of bhp as function of the input parameters.
* @param table_id Table number to use
* @param wells Wells structure with information about wells in qs
* @param qs Flow quantities
* @param thp Tubing head pressure
* @param alq Artificial lift or other parameter
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
ADB bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp,
const ADB& alq) const;
/**
* Linear interpolation of bhp as a function of the input parameters given as ADBs
* Each entry corresponds typically to one well.
* @param table_id Table number to use. A negative entry (e.g., -1)
* will indicate that no table is used, and the corresponding
* BHP will be calculated as a constant -1e100.
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param thp Tubing head pressure
* @param alq Artificial lift or other parameter
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table, for each entry in the
* input ADB objects.
*/
ADB bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp,
const ADB& alq) const;
};
} //namespace
#endif /* OPM_AUTODIFF_VFPPRODPROPERTIESLEGACY_HPP_ */

49
opm/autodiff/VFPProperties.cpp
View File

@@ -1,49 +0,0 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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"
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
namespace Opm {
VFPProperties::VFPProperties() {
}
VFPProperties::VFPProperties(const VFPInjTable* inj_table, const VFPProdTable* prod_table) {
if (inj_table != NULL) {
m_inj.reset(new VFPInjProperties(inj_table));
}
if (prod_table != NULL) {
m_prod.reset(new VFPProdProperties(prod_table));
}
}
VFPProperties::VFPProperties(const std::map<int, std::shared_ptr<const VFPInjTable> >& inj_tables,
const std::map<int, std::shared_ptr<const VFPProdTable> >& prod_tables) {
m_inj.reset(new VFPInjProperties(inj_tables));
m_prod.reset(new VFPProdProperties(prod_tables));
}
} //Namespace

27
opm/autodiff/VFPProperties.hpp
View File

@@ -27,16 +27,14 @@
namespace Opm {
class VFPProdProperties;
class VFPInjProperties;
/**
* A thin wrapper class that holds one VFPProdProperties and one
* VFPInjProperties object.
*/
template<typename VFPInjProp, typename VFPProdProp>
class VFPProperties {
public:
VFPProperties();
VFPProperties() {}
/**
* Constructor
@@ -44,7 +42,15 @@ public:
* @param inj_table A *single* VFPINJ table or NULL (no table)
* @param prod_table A *single* VFPPROD table or NULL (no table)
*/
explicit VFPProperties(const VFPInjTable* inj_table, const VFPProdTable* prod_table);
explicit VFPProperties(const VFPInjTable* inj_table, const VFPProdTable* prod_table)
{
if (inj_table != nullptr) {
m_inj.reset(new VFPInjProp(inj_table));
}
if (prod_table != nullptr) {
m_prod.reset(new VFPProdProp(prod_table));
}
}
/**
* Constructor
@@ -53,25 +59,26 @@ public:
* @param prod_tables A map of different VFPPROD tables.
*/
VFPProperties(const std::map<int, std::shared_ptr<const VFPInjTable> >& inj_tables,
const std::map<int, std::shared_ptr<const VFPProdTable> >& prod_tables);
const std::map<int, std::shared_ptr<const VFPProdTable> >& prod_tables)
: m_inj(new VFPInjProp(inj_tables)), m_prod(new VFPProdProp(prod_tables)) {}
/**
* Returns the VFP properties for injection wells
*/
const VFPInjProperties* getInj() const {
const VFPInjProp* getInj() const {
return m_inj.get();
}
/**
* Returns the VFP properties for production wells
*/
const VFPProdProperties* getProd() const {
const VFPProdProp* getProd() const {
return m_prod.get();
}
private:
std::shared_ptr<VFPInjProperties> m_inj;
std::shared_ptr<VFPProdProperties> m_prod;
std::shared_ptr<VFPInjProp> m_inj;
std::shared_ptr<VFPProdProp> m_prod;
};

6
opm/autodiff/WellInterface.hpp
View File

@@ -37,8 +37,6 @@
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/WellHelpers.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilModelParametersEbos.hpp>
@@ -130,7 +128,7 @@ namespace Opm
/// Well controls
WellControls* wellControls() const;
void setVFPProperties(const VFPProperties* vfp_properties_arg);
void setVFPProperties(const VFPProperties<VFPInjProperties,VFPProdProperties>* vfp_properties_arg);
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg,
@@ -287,7 +285,7 @@ namespace Opm
bool getAllowCrossFlow() const;
const VFPProperties* vfp_properties_;
const VFPProperties<VFPInjProperties,VFPProdProperties>* vfp_properties_;
double gravity_;

2
opm/autodiff/WellInterface_impl.hpp
View File

@@ -129,7 +129,7 @@ namespace Opm
template<typename TypeTag>
void
WellInterface<TypeTag>::
setVFPProperties(const VFPProperties* vfp_properties_arg)
setVFPProperties(const VFPProperties<VFPInjProperties,VFPProdProperties>* vfp_properties_arg)
{
vfp_properties_ = vfp_properties_arg;
}

448
tests/test_vfpproperties.cpp
View File

@@ -20,7 +20,7 @@
#include <config.h>
#define BOOST_TEST_MODULE AutoDiffBlockTest
#define BOOST_TEST_MODULE VFPTest
#include <algorithm>
#include <memory>
@@ -42,9 +42,8 @@
#include <opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/VFPHelpers.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
@@ -105,177 +104,6 @@ BOOST_AUTO_TEST_CASE(findInterpData)
}
BOOST_AUTO_TEST_SUITE_END() // HelperTests
struct ConversionFixture {
typedef Opm::VFPProdProperties::ADB ADB;
ConversionFixture() :
num_wells(5),
aqua(ADB::null()),
liquid(ADB::null()),
vapour(ADB::null())
{
ADB::V aqua_v(num_wells);
ADB::V liquid_v(num_wells);
ADB::V vapour_v(num_wells);
for (int i=0; i<num_wells; ++i) {
aqua_v[i] = 300+num_wells*15;
liquid_v[i] = 500+num_wells*15;
vapour_v[i] = 700+num_wells*15;
}
aqua = ADB::constant(aqua_v);
liquid = ADB::constant(liquid_v);
vapour = ADB::constant(vapour_v);
}
~ConversionFixture() {
}
int num_wells;
ADB aqua;
ADB liquid;
ADB vapour;
};
BOOST_FIXTURE_TEST_SUITE( ConversionTests, ConversionFixture )
BOOST_AUTO_TEST_CASE(getFlo)
{
//Compute reference solutions
std::vector<double> ref_flo_oil(num_wells);
std::vector<double> ref_flo_liq(num_wells);
std::vector<double> ref_flo_gas(num_wells);
for (int i=0; i<num_wells; ++i) {
ref_flo_oil[i] = liquid.value()[i];
ref_flo_liq[i] = aqua.value()[i] + liquid.value()[i];
ref_flo_gas[i] = vapour.value()[i];
}
{
ADB flo = Opm::detail::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_OIL);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_oil.begin(), ref_flo_oil.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_LIQ);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_liq.begin(), ref_flo_liq.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_GAS);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_gas.begin(), ref_flo_gas.end(), computed, computed+num_wells);
}
}
BOOST_AUTO_TEST_CASE(getWFR)
{
//Compute reference solutions
std::vector<double> ref_wfr_wor(num_wells);
std::vector<double> ref_wfr_wct(num_wells);
std::vector<double> ref_wfr_wgr(num_wells);
for (int i=0; i<num_wells; ++i) {
ref_wfr_wor[i] = aqua.value()[i] / liquid.value()[i];
ref_wfr_wct[i] = aqua.value()[i] / (aqua.value()[i] + liquid.value()[i]);
ref_wfr_wgr[i] = aqua.value()[i] / vapour.value()[i];
}
{
ADB flo = Opm::detail::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WOR);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_wfr_wor.begin(), ref_wfr_wor.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WCT);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_wfr_wct.begin(), ref_wfr_wct.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WGR);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_wfr_wgr.begin(), ref_wfr_wgr.end(), computed, computed+num_wells);
}
}
BOOST_AUTO_TEST_CASE(getGFR)
{
//Compute reference solutions
std::vector<double> ref_gfr_gor(num_wells);
std::vector<double> ref_gfr_glr(num_wells);
std::vector<double> ref_gfr_ogr(num_wells);
for (int i=0; i<num_wells; ++i) {
ref_gfr_gor[i] = vapour.value()[i] / liquid.value()[i];
ref_gfr_glr[i] = vapour.value()[i] / (liquid.value()[i] + aqua.value()[i]);
ref_gfr_ogr[i] = liquid.value()[i] / vapour.value()[i];
}
{
ADB flo = Opm::detail::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_GOR);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_gfr_gor.begin(), ref_gfr_gor.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_GLR);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_gfr_glr.begin(), ref_gfr_glr.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_OGR);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_gfr_ogr.begin(), ref_gfr_ogr.end(), computed, computed+num_wells);
}
}
BOOST_AUTO_TEST_SUITE_END() // unit tests
/**
@@ -284,7 +112,6 @@ BOOST_AUTO_TEST_SUITE_END() // unit tests
* values data is given at
*/
struct TrivialFixture {
typedef Opm::VFPProdProperties::ADB ADB;
typedef Opm::detail::VFPEvaluation VFPEvaluation;
TrivialFixture() : table_ids(1, 1),
@@ -391,15 +218,6 @@ struct TrivialFixture {
/**
* Helper function to simplify creating minimal ADB objects
*/
inline ADB createConstantScalarADB(double value) {
ADB::V v = ADB::V::Constant(1, value);
ADB adb = ADB::constant(std::move(v));
return adb;
}
std::shared_ptr<Opm::VFPProdProperties> properties;
std::shared_ptr<Opm::VFPProdTable> table;
std::vector<int> table_ids;
@@ -427,60 +245,6 @@ private:
BOOST_FIXTURE_TEST_SUITE( TrivialTests, TrivialFixture )
BOOST_AUTO_TEST_CASE(GetTable)
{
fillDataRandom();
initProperties();
//Create wells
const int nphases = 3;
const int nwells = 1;
const int nperfs = 1;
std::shared_ptr<Wells> wells(create_wells(nphases, nwells, nperfs),
destroy_wells);
const int cells[] = {5};
add_well(INJECTOR, 100, 1, NULL, cells, NULL, 0, NULL, true, wells.get());
//Create interpolation points
double aqua_d = -0.15;
double liquid_d = -0.25;
double vapour_d = -0.35;
double thp_d = 0.45;
double alq_d = 0.55;
ADB aqua_adb = createConstantScalarADB(aqua_d);
ADB liquid_adb = createConstantScalarADB(liquid_d);
ADB vapour_adb = createConstantScalarADB(vapour_d);
ADB thp_adb = createConstantScalarADB(thp_d);
ADB alq_adb = createConstantScalarADB(alq_d);
ADB::V qs_adb_v(3);
qs_adb_v << aqua_adb.value(), liquid_adb.value(), vapour_adb.value();
ADB qs_adb = ADB::constant(qs_adb_v);
//Check that our reference has not changed
Opm::detail::VFPEvaluation ref = Opm::detail::bhp(table.get(), aqua_d, liquid_d, vapour_d, thp_d, alq_d);
BOOST_CHECK_CLOSE(ref.value, 1.0923565702101556, max_d_tol);
BOOST_CHECK_CLOSE(ref.dthp, 0.13174065498177251, max_d_tol);
BOOST_CHECK_CLOSE(ref.dwfr, -1.2298177745501071, max_d_tol);
BOOST_CHECK_CLOSE(ref.dgfr, 0.82988935779290274, max_d_tol);
BOOST_CHECK_CLOSE(ref.dalq, 1.8148520254931713, max_d_tol);
BOOST_CHECK_CLOSE(ref.dflo, 9.0944843574181924, max_d_tol);
//Check that different versions of the prod_bph function work
ADB a = properties->bhp(table_ids, aqua_adb, liquid_adb, vapour_adb, thp_adb, alq_adb);
double b =properties->bhp(table_ids[0], aqua_d, liquid_d, vapour_d, thp_d, alq_d);
ADB c = properties->bhp(table_ids, *wells, qs_adb, thp_adb, alq_adb);
//Check that results are actually equal reference
BOOST_CHECK_EQUAL(a.value()[0], ref.value);
BOOST_CHECK_EQUAL(b, ref.value);
BOOST_CHECK_EQUAL(c.value()[0], ref.value);
//Table 2 does not exist.
std::vector<int> table_ids_wrong(1, 2);
BOOST_CHECK_THROW(properties->bhp(table_ids_wrong, *wells, qs_adb, thp_adb, alq_adb), std::invalid_argument);
}
/**
* Test that we can generate some dummy zero-data,
@@ -667,188 +431,6 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlane)
/**
* Test that we can generate some dummy data representing an ND plane,
* interpolate using ADBs as input, and compare against the analytic solution
*/
BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
{
fillDataPlane();
initProperties();
//Check linear extrapolation (i.e., using values of x, y, etc. outside our interpolant domain)
double sum = 0.0;
double reference_sum = 0.0;
double sad = 0.0; // Sum absolute difference
double max_d = 0.0; // Maximum difference
int n=1;
int o=5;
for (int i=0; i<=n+o; ++i) {
const double x = i / static_cast<double>(n);
for (int j=1; j<=n+o; ++j) {
const double aqua = -j / static_cast<double>(n);
for (int k=1; k<=n+o; ++k) {
const double vapour = -k / static_cast<double>(n);
for (int l=0; l<=n+o; ++l) {
const double u = l / static_cast<double>(n);
for (int m=1; m<=n+o; ++m) {
const double liquid = -m / static_cast<double>(n);
//Temporary variables used to represent independent wells
const int num_wells = 5;
ADB::V adb_v_x(num_wells);
ADB::V adb_v_aqua(num_wells);
ADB::V adb_v_vapour(num_wells);
ADB::V adb_v_u(num_wells);
ADB::V adb_v_liquid(num_wells);
table_ids.resize(num_wells);
for (unsigned int w=0; w<num_wells; ++w) {
table_ids[w] = 1;
adb_v_x[w] = x*(w+1);
adb_v_aqua[w] = aqua*(w+1);
adb_v_vapour[w] = vapour*(w+1);
adb_v_u[w] = u*(w+1);
adb_v_liquid[w] = liquid*(w+1);
}
ADB adb_x = ADB::constant(adb_v_x);
ADB adb_aqua = ADB::constant(adb_v_aqua);
ADB adb_vapour = ADB::constant(adb_v_vapour);
ADB adb_u = ADB::constant(adb_v_u);
ADB adb_liquid = ADB::constant(adb_v_liquid);
ADB bhp = properties->bhp(table_ids, adb_aqua, adb_liquid, adb_vapour, adb_x, adb_u);
ADB::V bhp_val = bhp.value();
double value = 0.0;
double reference = 0.0;
for (int w=0; w < num_wells; ++w) {
//Find values that should be in table
double v = Opm::detail::getFlo(aqua*(w+1), liquid*(w+1), vapour*(w+1), table->getFloType());
double y = Opm::detail::getWFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table->getWFRType());
double z = Opm::detail::getGFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table->getGFRType());
reference = x*(w+1) + 2*y + 3*z + 4*u*(w+1) - 5*v;
value = bhp_val[w];
sum += value;
reference_sum += reference;
double abs_diff = std::abs(value - reference);
sad += std::abs(abs_diff);
max_d = std::max(max_d, abs_diff);
}
}
}
}
}
}
BOOST_CHECK_CLOSE(sum, reference_sum, 0.0001);
BOOST_CHECK_SMALL(max_d, max_d_tol);
BOOST_CHECK_SMALL(sad, sad_tol);
}
/**
* Test that we can generate some dummy data representing an ND plane,
* interpolate using well flow rates and ADBs as input, and compare against the analytic solution
*/
BOOST_AUTO_TEST_CASE(InterpolateADBAndQs)
{
fillDataPlane();
initProperties();
//Create wells
const int nphases = 3;
const int nwells = 5;
const int nperfs = 1;
std::shared_ptr<Wells> wells(create_wells(nphases, nwells, nperfs),
destroy_wells);
int cells = 1;
for (int i=0; i<nwells; ++i) {
//Just give the cells a set of different indices
cells *= 2;
std::stringstream ss;
ss << "WELL_" << i;
const bool ok = add_well(INJECTOR, 0.0, 1, NULL, &cells,
NULL, 0, ss.str().c_str(), true, wells.get());
BOOST_REQUIRE(ok);
}
//Create some artificial flow values for our wells between 0 and 1
ADB::V qs_v(nphases*nwells);
for (int j=0; j<nphases; ++j) {
for (int i=0; i<nwells; ++i) {
qs_v[j*nwells+i] = -(j*nwells+i) / static_cast<double>(nwells*nphases-1.0);
}
}
ADB qs = ADB::constant(qs_v);
//Create the THP for each well
ADB::V thp_v(nwells);
for (int i=0; i<nwells; ++i) {
thp_v[i] = (i) / static_cast<double>(nwells-1.0);
}
ADB thp = ADB::constant(thp_v);
//Create the ALQ for each well
ADB::V alq_v(nwells);
for (int i=0; i<nwells; ++i) {
alq_v[i] = 0.0;
}
ADB alq = ADB::constant(alq_v);
//Set which VFP table to use for each well
table_ids.resize(nwells);
for (int i=0; i<nwells; ++i) {
table_ids[i] = 1;
}
//Call the bhp function
ADB::V bhp = properties->bhp(table_ids, *wells, qs, thp, alq).value();
//Calculate reference
//First, find the three phases
std::vector<double> water(nwells);
std::vector<double> oil(nwells);
std::vector<double> gas(nwells);
for (int i=0; i<nwells; ++i) {
water[i] = qs_v[i];
oil[i] = qs_v[nwells+i];
gas[i] = qs_v[2*nwells+i];
}
//Compute reference value
std::vector<double> reference(nwells);
for (int i=0; i<nwells; ++i) {
double flo = oil[i];
double wor = water[i]/oil[i];
double gor = gas[i]/oil[i];
reference[i] = thp_v[i] + 2*wor + 3*gor + 4*alq_v[i] - 5*flo;
}
//Check that interpolation matches
BOOST_REQUIRE_EQUAL(bhp.size(), nwells);
double sad = 0.0;
double max_d = 0.0;
for (int i=0; i<nwells; ++i) {
double value = bhp[i];
double ref = reference[i];
double abs_diff = std::abs(value-ref);
sad += abs_diff;
max_d = std::max(abs_diff, max_d);
}
BOOST_CHECK_SMALL(max_d, max_d_tol);
BOOST_CHECK_SMALL(sad, sad_tol);
}
/**
* Test that the partial derivatives are reasonable
@@ -976,32 +558,6 @@ BOOST_AUTO_TEST_CASE(THPToBHPAndBackNonTrivial)
BOOST_AUTO_TEST_SUITE_END() // Trivial tests
BOOST_AUTO_TEST_SUITE(IntegrationTests)
extern const double reference[];

605
tests/test_vfpproperties_legacy.cpp Normal file
View File

@@ -0,0 +1,605 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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>
#define BOOST_TEST_MODULE VFPTest
#include <algorithm>
#include <memory>
#include <map>
#include <sstream>
#include <limits>
#include <vector>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/test/unit_test.hpp>
#include <boost/filesystem.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/core/wells.h>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <opm/autodiff/VFPHelpersLegacy.hpp>
#include <opm/autodiff/VFPProdPropertiesLegacy.hpp>
const double max_d_tol = 1.0e-10;
const double sad_tol = 1.0e-8;
struct ConversionFixture {
typedef Opm::VFPProdPropertiesLegacy::ADB ADB;
ConversionFixture() :
num_wells(5),
aqua(ADB::null()),
liquid(ADB::null()),
vapour(ADB::null())
{
ADB::V aqua_v(num_wells);
ADB::V liquid_v(num_wells);
ADB::V vapour_v(num_wells);
for (int i=0; i<num_wells; ++i) {
aqua_v[i] = 300+num_wells*15;
liquid_v[i] = 500+num_wells*15;
vapour_v[i] = 700+num_wells*15;
}
aqua = ADB::constant(aqua_v);
liquid = ADB::constant(liquid_v);
vapour = ADB::constant(vapour_v);
}
~ConversionFixture() {
}
int num_wells;
ADB aqua;
ADB liquid;
ADB vapour;
};
BOOST_FIXTURE_TEST_SUITE( ConversionTests, ConversionFixture )
BOOST_AUTO_TEST_CASE(getFlo)
{
//Compute reference solutions
std::vector<double> ref_flo_oil(num_wells);
std::vector<double> ref_flo_liq(num_wells);
std::vector<double> ref_flo_gas(num_wells);
for (int i=0; i<num_wells; ++i) {
ref_flo_oil[i] = liquid.value()[i];
ref_flo_liq[i] = aqua.value()[i] + liquid.value()[i];
ref_flo_gas[i] = vapour.value()[i];
}
{
ADB flo = Opm::detail::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_OIL);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_oil.begin(), ref_flo_oil.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_LIQ);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_liq.begin(), ref_flo_liq.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_GAS);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_gas.begin(), ref_flo_gas.end(), computed, computed+num_wells);
}
}
BOOST_AUTO_TEST_CASE(getWFR)
{
//Compute reference solutions
std::vector<double> ref_wfr_wor(num_wells);
std::vector<double> ref_wfr_wct(num_wells);
std::vector<double> ref_wfr_wgr(num_wells);
for (int i=0; i<num_wells; ++i) {
ref_wfr_wor[i] = aqua.value()[i] / liquid.value()[i];
ref_wfr_wct[i] = aqua.value()[i] / (aqua.value()[i] + liquid.value()[i]);
ref_wfr_wgr[i] = aqua.value()[i] / vapour.value()[i];
}
{
ADB flo = Opm::detail::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WOR);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_wfr_wor.begin(), ref_wfr_wor.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WCT);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_wfr_wct.begin(), ref_wfr_wct.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WGR);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_wfr_wgr.begin(), ref_wfr_wgr.end(), computed, computed+num_wells);
}
}
BOOST_AUTO_TEST_CASE(getGFR)
{
//Compute reference solutions
std::vector<double> ref_gfr_gor(num_wells);
std::vector<double> ref_gfr_glr(num_wells);
std::vector<double> ref_gfr_ogr(num_wells);
for (int i=0; i<num_wells; ++i) {
ref_gfr_gor[i] = vapour.value()[i] / liquid.value()[i];
ref_gfr_glr[i] = vapour.value()[i] / (liquid.value()[i] + aqua.value()[i]);
ref_gfr_ogr[i] = liquid.value()[i] / vapour.value()[i];
}
{
ADB flo = Opm::detail::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_GOR);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_gfr_gor.begin(), ref_gfr_gor.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_GLR);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_gfr_glr.begin(), ref_gfr_glr.end(), computed, computed+num_wells);
}
{
ADB flo = Opm::detail::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_OGR);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_gfr_ogr.begin(), ref_gfr_ogr.end(), computed, computed+num_wells);
}
}
BOOST_AUTO_TEST_SUITE_END() // unit tests
/**
* Test fixture to set up axis etc.
* All of our axes go from 0 to 1, but with a varying number of
* values data is given at
*/
struct TrivialFixture {
typedef Opm::VFPProdPropertiesLegacy::ADB ADB;
typedef Opm::detail::VFPEvaluation VFPEvaluation;
TrivialFixture() : table_ids(1, 1),
thp_axis{0.0, 1.0},
wfr_axis{0.0, 0.5, 1.0},
gfr_axis{0.0, 0.25, 0.5, 0.75, 1},
alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1},
flo_axis{0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375,
0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1},
nx(thp_axis.size()),
ny(wfr_axis.size()),
nz(gfr_axis.size()),
nu(alq_axis.size()),
nv(flo_axis.size()),
size{{ nx, ny, nz, nu, nv }},
data(size)
{
}
~TrivialFixture() {
}
/**
* Fills our interpolation data with zeros
*/
inline void fillData(double value) {
for (int i=0; i<nx; ++i) {
for (int j=0; j<ny; ++j) {
for (int k=0; k<nz; ++k) {
for (int l=0; l<nu; ++l) {
for (int m=0; m<nv; ++m) {
data[i][j][k][l][m] = value;
}
}
}
}
}
}
/**
* Fills our interpolation data with an ND plane
*/
inline void fillDataPlane() {
for (int i=0; i<nx; ++i) {
double x = i / static_cast<double>(nx-1);
for (int j=0; j<ny; ++j) {
double y = j / static_cast<double>(ny-1);
for (int k=0; k<nz; ++k) {
double z = k / static_cast<double>(nz-1);
for (int l=0; l<nu; ++l) {
double u = l / static_cast<double>(nu-1);
for (int m=0; m<nv; ++m) {
double v = m / static_cast<double>(nv-1);
// table[thp_idx][wfr_idx][gfr_idx][alq_idx][flo_idx];
data[i][j][k][l][m] = x + 2*y + 3*z + 4*u + 5*v;
}
}
}
}
}
}
/**
* Fills our interpolation data with "random" values
*/
inline void fillDataRandom() {
unsigned long randx = 42;
static double max_val = static_cast<double>(std::numeric_limits<unsigned long>::max());
for (int i=0; i<nx; ++i) {
for (int j=0; j<ny; ++j) {
for (int k=0; k<nz; ++k) {
for (int l=0; l<nu; ++l) {
for (int m=0; m<nv; ++m) {
data[i][j][k][l][m] = randx / max_val;
randx = (randx*1103515245 + 12345);
}
}
}
}
}
}
inline void initProperties() {
table.reset(new Opm::VFPProdTable(1,
1000.0,
Opm::VFPProdTable::FLO_OIL,
Opm::VFPProdTable::WFR_WOR,
Opm::VFPProdTable::GFR_GOR,
Opm::VFPProdTable::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data));
properties.reset(new Opm::VFPProdPropertiesLegacy(table.get()));
}
/**
* Helper function to simplify creating minimal ADB objects
*/
inline ADB createConstantScalarADB(double value) {
ADB::V v = ADB::V::Constant(1, value);
ADB adb = ADB::constant(std::move(v));
return adb;
}
std::shared_ptr<Opm::VFPProdPropertiesLegacy> properties;
std::shared_ptr<Opm::VFPProdTable> table;
std::vector<int> table_ids;
private:
const std::vector<double> thp_axis;
const std::vector<double> wfr_axis;
const std::vector<double> gfr_axis;
const std::vector<double> alq_axis;
const std::vector<double> flo_axis;
int nx;
int ny;
int nz;
int nu;
int nv;
Opm::VFPProdTable::extents size;
Opm::VFPProdTable::array_type data;
};
//Set F to be our test suite fixture for our "trivial" tests
BOOST_FIXTURE_TEST_SUITE( TrivialTests, TrivialFixture )
BOOST_AUTO_TEST_CASE(GetTable)
{
fillDataRandom();
initProperties();
//Create wells
const int nphases = 3;
const int nwells = 1;
const int nperfs = 1;
std::shared_ptr<Wells> wells(create_wells(nphases, nwells, nperfs),
destroy_wells);
const int cells[] = {5};
add_well(INJECTOR, 100, 1, NULL, cells, NULL, 0, NULL, true, wells.get());
//Create interpolation points
double aqua_d = -0.15;
double liquid_d = -0.25;
double vapour_d = -0.35;
double thp_d = 0.45;
double alq_d = 0.55;
ADB aqua_adb = createConstantScalarADB(aqua_d);
ADB liquid_adb = createConstantScalarADB(liquid_d);
ADB vapour_adb = createConstantScalarADB(vapour_d);
ADB thp_adb = createConstantScalarADB(thp_d);
ADB alq_adb = createConstantScalarADB(alq_d);
ADB::V qs_adb_v(3);
qs_adb_v << aqua_adb.value(), liquid_adb.value(), vapour_adb.value();
ADB qs_adb = ADB::constant(qs_adb_v);
//Check that our reference has not changed
Opm::detail::VFPEvaluation ref = Opm::detail::bhp(table.get(), aqua_d, liquid_d, vapour_d, thp_d, alq_d);
BOOST_CHECK_CLOSE(ref.value, 1.0923565702101556, max_d_tol);
BOOST_CHECK_CLOSE(ref.dthp, 0.13174065498177251, max_d_tol);
BOOST_CHECK_CLOSE(ref.dwfr, -1.2298177745501071, max_d_tol);
BOOST_CHECK_CLOSE(ref.dgfr, 0.82988935779290274, max_d_tol);
BOOST_CHECK_CLOSE(ref.dalq, 1.8148520254931713, max_d_tol);
BOOST_CHECK_CLOSE(ref.dflo, 9.0944843574181924, max_d_tol);
//Check that different versions of the prod_bph function work
ADB a = properties->bhp(table_ids, aqua_adb, liquid_adb, vapour_adb, thp_adb, alq_adb);
double b =properties->bhp(table_ids[0], aqua_d, liquid_d, vapour_d, thp_d, alq_d);
ADB c = properties->bhp(table_ids, *wells, qs_adb, thp_adb, alq_adb);
//Check that results are actually equal reference
BOOST_CHECK_EQUAL(a.value()[0], ref.value);
BOOST_CHECK_EQUAL(b, ref.value);
BOOST_CHECK_EQUAL(c.value()[0], ref.value);
//Table 2 does not exist.
std::vector<int> table_ids_wrong(1, 2);
BOOST_CHECK_THROW(properties->bhp(table_ids_wrong, *wells, qs_adb, thp_adb, alq_adb), std::invalid_argument);
}
/**
* Test that we can generate some dummy data representing an ND plane,
* interpolate using ADBs as input, and compare against the analytic solution
*/
BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
{
fillDataPlane();
initProperties();
//Check linear extrapolation (i.e., using values of x, y, etc. outside our interpolant domain)
double sum = 0.0;
double reference_sum = 0.0;
double sad = 0.0; // Sum absolute difference
double max_d = 0.0; // Maximum difference
int n=1;
int o=5;
for (int i=0; i<=n+o; ++i) {
const double x = i / static_cast<double>(n);
for (int j=1; j<=n+o; ++j) {
const double aqua = -j / static_cast<double>(n);
for (int k=1; k<=n+o; ++k) {
const double vapour = -k / static_cast<double>(n);
for (int l=0; l<=n+o; ++l) {
const double u = l / static_cast<double>(n);
for (int m=1; m<=n+o; ++m) {
const double liquid = -m / static_cast<double>(n);
//Temporary variables used to represent independent wells
const int num_wells = 5;
ADB::V adb_v_x(num_wells);
ADB::V adb_v_aqua(num_wells);
ADB::V adb_v_vapour(num_wells);
ADB::V adb_v_u(num_wells);
ADB::V adb_v_liquid(num_wells);
table_ids.resize(num_wells);
for (unsigned int w=0; w<num_wells; ++w) {
table_ids[w] = 1;
adb_v_x[w] = x*(w+1);
adb_v_aqua[w] = aqua*(w+1);
adb_v_vapour[w] = vapour*(w+1);
adb_v_u[w] = u*(w+1);
adb_v_liquid[w] = liquid*(w+1);
}
ADB adb_x = ADB::constant(adb_v_x);
ADB adb_aqua = ADB::constant(adb_v_aqua);
ADB adb_vapour = ADB::constant(adb_v_vapour);
ADB adb_u = ADB::constant(adb_v_u);
ADB adb_liquid = ADB::constant(adb_v_liquid);
ADB bhp = properties->bhp(table_ids, adb_aqua, adb_liquid, adb_vapour, adb_x, adb_u);
ADB::V bhp_val = bhp.value();
double value = 0.0;
double reference = 0.0;
for (int w=0; w < num_wells; ++w) {
//Find values that should be in table
double v = Opm::detail::getFlo(aqua*(w+1), liquid*(w+1), vapour*(w+1), table->getFloType());
double y = Opm::detail::getWFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table->getWFRType());
double z = Opm::detail::getGFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table->getGFRType());
reference = x*(w+1) + 2*y + 3*z + 4*u*(w+1) - 5*v;
value = bhp_val[w];
sum += value;
reference_sum += reference;
double abs_diff = std::abs(value - reference);
sad += std::abs(abs_diff);
max_d = std::max(max_d, abs_diff);
}
}
}
}
}
}
BOOST_CHECK_CLOSE(sum, reference_sum, 0.0001);
BOOST_CHECK_SMALL(max_d, max_d_tol);
BOOST_CHECK_SMALL(sad, sad_tol);
}
/**
* Test that we can generate some dummy data representing an ND plane,
* interpolate using well flow rates and ADBs as input, and compare against the analytic solution
*/
BOOST_AUTO_TEST_CASE(InterpolateADBAndQs)
{
fillDataPlane();
initProperties();
//Create wells
const int nphases = 3;
const int nwells = 5;
const int nperfs = 1;
std::shared_ptr<Wells> wells(create_wells(nphases, nwells, nperfs),
destroy_wells);
int cells = 1;
for (int i=0; i<nwells; ++i) {
//Just give the cells a set of different indices
cells *= 2;
std::stringstream ss;
ss << "WELL_" << i;
const bool ok = add_well(INJECTOR, 0.0, 1, NULL, &cells,
NULL, 0, ss.str().c_str(), true, wells.get());
BOOST_REQUIRE(ok);
}
//Create some artificial flow values for our wells between 0 and 1
ADB::V qs_v(nphases*nwells);
for (int j=0; j<nphases; ++j) {
for (int i=0; i<nwells; ++i) {
qs_v[j*nwells+i] = -(j*nwells+i) / static_cast<double>(nwells*nphases-1.0);
}
}
ADB qs = ADB::constant(qs_v);
//Create the THP for each well
ADB::V thp_v(nwells);
for (int i=0; i<nwells; ++i) {
thp_v[i] = (i) / static_cast<double>(nwells-1.0);
}
ADB thp = ADB::constant(thp_v);
//Create the ALQ for each well
ADB::V alq_v(nwells);
for (int i=0; i<nwells; ++i) {
alq_v[i] = 0.0;
}
ADB alq = ADB::constant(alq_v);
//Set which VFP table to use for each well
table_ids.resize(nwells);
for (int i=0; i<nwells; ++i) {
table_ids[i] = 1;
}
//Call the bhp function
ADB::V bhp = properties->bhp(table_ids, *wells, qs, thp, alq).value();
//Calculate reference
//First, find the three phases
std::vector<double> water(nwells);
std::vector<double> oil(nwells);
std::vector<double> gas(nwells);
for (int i=0; i<nwells; ++i) {
water[i] = qs_v[i];
oil[i] = qs_v[nwells+i];
gas[i] = qs_v[2*nwells+i];
}
//Compute reference value
std::vector<double> reference(nwells);
for (int i=0; i<nwells; ++i) {
double flo = oil[i];
double wor = water[i]/oil[i];
double gor = gas[i]/oil[i];
reference[i] = thp_v[i] + 2*wor + 3*gor + 4*alq_v[i] - 5*flo;
}
//Check that interpolation matches
BOOST_REQUIRE_EQUAL(bhp.size(), nwells);
double sad = 0.0;
double max_d = 0.0;
for (int i=0; i<nwells; ++i) {
double value = bhp[i];
double ref = reference[i];
double abs_diff = std::abs(value-ref);
sad += abs_diff;
max_d = std::max(abs_diff, max_d);
}
BOOST_CHECK_SMALL(max_d, max_d_tol);
BOOST_CHECK_SMALL(sad, sad_tol);
}
BOOST_AUTO_TEST_SUITE_END() // Trivial tests