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

changed: split VFP classes in base (ebos) and legacy

Browse Source 
to get rid of eigen usage in ebos based classes
This commit is contained in:
Arne Morten Kvarving 2018-11-13 10:45:02 +01:00
parent 54cc52dcaf
commit c5ae3adbbf
22 changed files with 825 additions and 689 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
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

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;
}

13
tests/test_vfpproperties.cpp
View File

@ -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/VFPHelpersLegacy.hpp>
#include <opm/autodiff/VFPProdPropertiesLegacy.hpp>
@ -115,7 +114,7 @@ BOOST_AUTO_TEST_SUITE_END() // HelperTests
struct ConversionFixture {
typedef Opm::VFPProdProperties::ADB ADB;
typedef Opm::VFPProdPropertiesLegacy::ADB ADB;
ConversionFixture() :
num_wells(5),
@ -284,7 +283,7 @@ BOOST_AUTO_TEST_SUITE_END() // unit tests
* values data is given at
*/
struct TrivialFixture {
typedef Opm::VFPProdProperties::ADB ADB;
typedef Opm::VFPProdPropertiesLegacy::ADB ADB;
typedef Opm::detail::VFPEvaluation VFPEvaluation;
TrivialFixture() : table_ids(1, 1),
@ -386,7 +385,7 @@ struct TrivialFixture {
alq_axis,
data));
properties.reset(new Opm::VFPProdProperties(table.get()));
properties.reset(new Opm::VFPProdPropertiesLegacy(table.get()));
}
@ -400,7 +399,7 @@ struct TrivialFixture {
return adb;
}
std::shared_ptr<Opm::VFPProdProperties> properties;
std::shared_ptr<Opm::VFPProdPropertiesLegacy> properties;
std::shared_ptr<Opm::VFPProdTable> table;
std::vector<int> table_ids;