diff --git a/opm/common/utility/VectorWithDefaultAllocator.hpp b/opm/common/utility/VectorWithDefaultAllocator.hpp index 8f114a243bb..7e0e667b60e 100644 --- a/opm/common/utility/VectorWithDefaultAllocator.hpp +++ b/opm/common/utility/VectorWithDefaultAllocator.hpp @@ -18,6 +18,8 @@ #ifndef OPM_COMMON_VECTOR_WITH_DEFAULT_ALLOCATOR_HPP #define OPM_COMMON_VECTOR_WITH_DEFAULT_ALLOCATOR_HPP + +#include namespace Opm { // NVCC being weird about std::vector, so we need this workaround. diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp index 507911426f9..c6e8887834e 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp @@ -37,14 +37,29 @@ namespace Opm { -template -EclMaterialLawManager::EclMaterialLawManager() = default; - -template -EclMaterialLawManager::~EclMaterialLawManager() = default; - -template -void EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +EclMaterialLawManager::EclMaterialLawManager() = default; + +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +EclMaterialLawManager::~EclMaterialLawManager() = default; + +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +void EclMaterialLawManager:: initFromState(const EclipseState& eclState) { // get the number of saturation regions and the number of cells in the deck @@ -134,10 +149,15 @@ initFromState(const EclipseState& eclState) } } -template -void EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +void EclMaterialLawManager:: initParamsForElements(const EclipseState& eclState, size_t numCompressedElems, - const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner, + const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner, const std::function& lookupIdxOnLevelZeroAssigner) { InitParams initParams {*this, eclState, numCompressedElems}; @@ -146,8 +166,13 @@ initParamsForElements(const EclipseState& eclState, size_t numCompressedElems, // TODO: Better (proper?) handling of mixed wettability systems - see ecl kw OPTIONS switch 74 // Note: Without OPTIONS[74] the negative part of the Pcow curve is not scaled -template -std::pair EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +std::pair EclMaterialLawManager:: applySwatinit(unsigned elemIdx, Scalar pcow, Scalar Sw) @@ -220,9 +245,14 @@ applySwatinit(unsigned elemIdx, return {Sw, newSwatInit}; } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::applyRestartSwatInit(const unsigned elemIdx, +EclMaterialLawManager::applyRestartSwatInit(const unsigned elemIdx, const Scalar maxPcow) { // Maximum capillary pressure adjusted from SWATINIT data. @@ -238,9 +268,14 @@ EclMaterialLawManager::applyRestartSwatInit(const unsigned elemIdx, EclTwoPhaseSystemType::OilWater); } -template -const typename EclMaterialLawManager::MaterialLawParams& -EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +const typename EclMaterialLawManager::MaterialLawParams& +EclMaterialLawManager:: connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const { MaterialLawParams& mlp = const_cast(materialLawParams_[elemIdx]); @@ -328,12 +363,17 @@ connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const return mlp; } -template -int EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +int EclMaterialLawManager:: getKrnumSatIdx(unsigned elemIdx, FaceDir::DirEnum facedir) const { using Dir = FaceDir::DirEnum; - const std::vector* array = nullptr; + const Storage* array = nullptr; switch(facedir) { case Dir::XPlus: array = &krnumXArray_; @@ -355,8 +395,13 @@ getKrnumSatIdx(unsigned elemIdx, FaceDir::DirEnum facedir) const } } -template -void EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +void EclMaterialLawManager:: oilWaterHysteresisParams(Scalar& soMax, Scalar& swMax, Scalar& swMin, @@ -369,8 +414,13 @@ oilWaterHysteresisParams(Scalar& soMax, MaterialLaw::oilWaterHysteresisParams(soMax, swMax, swMin, params); } -template -void EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +void EclMaterialLawManager:: setOilWaterHysteresisParams(const Scalar& soMax, const Scalar& swMax, const Scalar& swMin, @@ -383,8 +433,13 @@ setOilWaterHysteresisParams(const Scalar& soMax, MaterialLaw::setOilWaterHysteresisParams(soMax, swMax, swMin, params); } -template -void EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +void EclMaterialLawManager:: gasOilHysteresisParams(Scalar& sgmax, Scalar& shmax, Scalar& somin, @@ -397,8 +452,13 @@ gasOilHysteresisParams(Scalar& sgmax, MaterialLaw::gasOilHysteresisParams(sgmax, shmax, somin, params); } -template -void EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +void EclMaterialLawManager:: setGasOilHysteresisParams(const Scalar& sgmax, const Scalar& shmax, const Scalar& somin, @@ -411,11 +471,17 @@ setGasOilHysteresisParams(const Scalar& sgmax, MaterialLaw::setGasOilHysteresisParams(sgmax, shmax, somin, params); } -template -EclEpsScalingPoints& -EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +EclEpsScalingPoints& +EclMaterialLawManager:: oilWaterScaledEpsPointsDrainage(unsigned elemIdx) { + #if !HAVE_CUDA auto& materialParams = materialLawParams_[elemIdx]; switch (materialParams.approach()) { case EclMultiplexerApproach::Stone1: { @@ -440,10 +506,18 @@ oilWaterScaledEpsPointsDrainage(unsigned elemIdx) default: throw std::logic_error("Enum value for material approach unknown!"); } + #else + OPM_THROW(NotImplementedError, "CUDA support is not available"); + #endif } -template -const typename EclMaterialLawManager::MaterialLawParams& EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +const typename EclMaterialLawManager::MaterialLawParams& EclMaterialLawManager:: materialLawParamsFunc_(unsigned elemIdx, FaceDir::DirEnum facedir) const { using Dir = FaceDir::DirEnum; @@ -467,8 +541,13 @@ materialLawParamsFunc_(unsigned elemIdx, FaceDir::DirEnum facedir) const } } -template -void EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +void EclMaterialLawManager:: readGlobalEpsOptions_(const EclipseState& eclState) { oilWaterEclEpsConfig_ = std::make_shared(); @@ -477,16 +556,26 @@ readGlobalEpsOptions_(const EclipseState& eclState) enableEndPointScaling_ = eclState.getTableManager().hasTables("ENKRVD"); } -template -void EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +void EclMaterialLawManager:: readGlobalHysteresisOptions_(const EclipseState& state) { hysteresisConfig_ = std::make_shared(); hysteresisConfig_->initFromState(state.runspec()); } -template -void EclMaterialLawManager:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +void EclMaterialLawManager:: readGlobalThreePhaseOptions_(const Runspec& runspec) { bool gasEnabled = runspec.phases().active(Phase::GAS); diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp index 184ffecf1a1..222db8339f2 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp @@ -31,6 +31,11 @@ #ifndef OPM_ECL_MATERIAL_LAW_MANAGER_HPP #define OPM_ECL_MATERIAL_LAW_MANAGER_HPP + +#include "EclTwoPhaseMaterial.hpp" // For GPU version we will directly access this material layer + +#include + #include #include @@ -69,7 +74,12 @@ class TableColumn; * \brief Provides an simple way to create and manage the material law objects * for a complete ECL deck. */ -template +template< + class TraitsT, + template class Storage = VectorWithDefaultAllocator, + template typename SharedPtr = std::shared_ptr, + template typename UniquePtr = std::unique_ptr +> class EclMaterialLawManager { private: @@ -84,6 +94,8 @@ class EclMaterialLawManager using OilWaterTraits = TwoPhaseMaterialTraits; using GasWaterTraits = TwoPhaseMaterialTraits; +#if !HAVE_CUDA + // the two-phase material law which is defined on effective (unscaled) saturations using GasOilEffectiveTwoPhaseLaw = SatCurveMultiplexer; using OilWaterEffectiveTwoPhaseLaw = SatCurveMultiplexer; @@ -109,84 +121,113 @@ class EclMaterialLawManager using OilWaterTwoPhaseHystParams = typename OilWaterTwoPhaseLaw::Params; using GasWaterTwoPhaseHystParams = typename GasWaterTwoPhaseLaw::Params; +#else // HAVE_CUDA + + using GasOilEffectiveTwoPhaseLaw = PiecewiseLinearTwoPhaseMaterial>>; + using OilWaterEffectiveTwoPhaseLaw = PiecewiseLinearTwoPhaseMaterial>>; + using GasWaterEffectiveTwoPhaseLaw = PiecewiseLinearTwoPhaseMaterial>>; + + using GasOilEffectiveTwoPhaseParams = typename GasOilEffectiveTwoPhaseLaw::Params; + using OilWaterEffectiveTwoPhaseParams = typename OilWaterEffectiveTwoPhaseLaw::Params; + using GasWaterEffectiveTwoPhaseParams = typename GasWaterEffectiveTwoPhaseLaw::Params; + +#endif // !HAVE_CUDA + public: // the three-phase material law used by the simulation + #if !HAVE_CUDA using MaterialLaw = EclMultiplexerMaterial; + #else // HAVE_CUDA + using MaterialLaw = EclTwoPhaseMaterial; + #endif using MaterialLawParams = typename MaterialLaw::Params; - using DirectionalMaterialLawParamsPtr = std::unique_ptr>; + using DirectionalMaterialLawParamsPtr = UniquePtr>; EclMaterialLawManager(); ~EclMaterialLawManager(); private: // internal typedefs - using GasOilEffectiveParamVector = std::vector>; - using OilWaterEffectiveParamVector = std::vector>; - using GasWaterEffectiveParamVector = std::vector>; - - using GasOilScalingPointsVector = std::vector>>; - using OilWaterScalingPointsVector = std::vector>>; - using GasWaterScalingPointsVector = std::vector>>; - using OilWaterScalingInfoVector = std::vector>; - using GasOilParamVector = std::vector>; - using OilWaterParamVector = std::vector>; - using GasWaterParamVector = std::vector>; - using MaterialLawParamsVector = std::vector>; + using GasOilEffectiveParamVector = Storage>; + using OilWaterEffectiveParamVector = Storage>; + using GasWaterEffectiveParamVector = Storage>; + + #if !HAVE_CUDA + using GasOilScalingPointsVector = Storage>>; + using OilWaterScalingPointsVector = Storage>>; + using GasWaterScalingPointsVector = Storage>>; + using OilWaterScalingInfoVector = Storage>; + using GasOilParamVector = Storage>; + using OilWaterParamVector = Storage>; + using GasWaterParamVector = Storage>; + using MaterialLawParamsVector = Storage>; + #endif // !HAVE_CUDA // helper classes // This class' implementation is defined in "EclMaterialLawManagerInitParams.cpp" class InitParams { public: - InitParams(EclMaterialLawManager& parent, const EclipseState& eclState, size_t numCompressedElems); + InitParams(EclMaterialLawManager& parent, const EclipseState& eclState, size_t numCompressedElems); // \brief Function argument 'fieldPropIntOnLeadAssigner' needed to lookup // field properties of cells on the leaf grid view for CpGrid with local grid refinement. // Function argument 'lookupIdxOnLevelZeroAssigner' is added to lookup, for each // leaf gridview cell with index 'elemIdx', its 'lookupIdx' (index of the parent/equivalent cell on level zero). - void run(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner, + void run(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner, const std::function& lookupIdxOnLevelZeroAssigner); private: + + #if !HAVE_CUDA class HystParams; + #endif // !HAVE_CUDA // \brief Function argument 'fieldPropIntOnLeadAssigner' needed to lookup // field properties of cells on the leaf grid view for CpGrid with local grid refinement. - void copySatnumArrays_(const std::function(const FieldPropsManager&, const std::string&, bool)>& + void copySatnumArrays_(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner); // \brief Function argument 'fieldPropIntOnLeadAssigner' needed to lookup // field properties of cells on the leaf grid view for CpGrid with local grid refinement. - void copyIntArray_(std::vector& dest, const std::string& keyword, - const std::function(const FieldPropsManager&, const std::string&, bool)>& + void copyIntArray_(Storage& dest, const std::string& keyword, + const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner); - unsigned imbRegion_(std::vector& array, unsigned elemIdx); + unsigned imbRegion_(Storage& array, unsigned elemIdx); void initArrays_( - std::vector*>& satnumArray, - std::vector*>& imbnumArray, - std::vector*>& mlpArray); + Storage*>& satnumArray, + Storage*>& imbnumArray, + Storage*>& mlpArray); void initMaterialLawParamVectors_(); void initOilWaterScaledEpsInfo_(); // \brief Function argument 'fieldProptOnLeadAssigner' needed to lookup // field properties of cells on the leaf grid view for CpGrid with local grid refinement. - void initSatnumRegionArray_(const std::function(const FieldPropsManager&, const std::string&, bool)>& + void initSatnumRegionArray_(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner); + #if !HAVE_CUDA void initThreePhaseParams_( HystParams &hystParams, MaterialLawParams& materialParams, unsigned satRegionIdx, unsigned elemIdx); + #else // HAVE_CUDA + void initThreePhaseParams_( + MaterialLawParams& materialParams, + unsigned satRegionIdx, + unsigned elemIdx); + #endif void readEffectiveParameters_(); void readUnscaledEpsPointsVectors_(); template - void readUnscaledEpsPoints_(Container& dest, std::shared_ptr config, EclTwoPhaseSystemType system_type); - unsigned satRegion_(std::vector& array, unsigned elemIdx); - unsigned satOrImbRegion_(std::vector& array, std::vector& default_vec, unsigned elemIdx); + void readUnscaledEpsPoints_(Container& dest, SharedPtr config, EclTwoPhaseSystemType system_type); + unsigned satRegion_(Storage& array, unsigned elemIdx); + unsigned satOrImbRegion_(Storage& array, Storage& default_vec, unsigned elemIdx); + #if !HAVE_CUDA // This class' implementation is defined in "EclMaterialLawManagerHystParams.cpp" class HystParams { public: - explicit HystParams(EclMaterialLawManager::InitParams& init_params); + explicit HystParams(EclMaterialLawManager::InitParams& init_params); void finalize(); - std::shared_ptr getGasOilParams(); - std::shared_ptr getOilWaterParams(); - std::shared_ptr getGasWaterParams(); + SharedPtr getGasOilParams(); + SharedPtr getOilWaterParams(); + SharedPtr getGasWaterParams(); void setConfig(unsigned satRegionIdx); // Function argument 'lookupIdxOnLevelZeroAssigner' is added to lookup, for each // leaf gridview cell with index 'elemIdx', its 'lookupIdx' (index of the parent/equivalent cell on level zero). @@ -219,21 +260,22 @@ class EclMaterialLawManager readScaledEpsPointsImbibition_(unsigned elemIdx, EclTwoPhaseSystemType type, const std::function& lookupIdxOnLevelZeroAssigner); - EclMaterialLawManager::InitParams& init_params_; - EclMaterialLawManager& parent_; + EclMaterialLawManager::InitParams& init_params_; + EclMaterialLawManager& parent_; const EclipseState& eclState_; - std::shared_ptr gasOilParams_; - std::shared_ptr oilWaterParams_; - std::shared_ptr gasWaterParams_; + SharedPtr gasOilParams_; + SharedPtr oilWaterParams_; + SharedPtr gasWaterParams_; }; + #endif // !HAVE_CUDA // This class' implementation is defined in "EclMaterialLawManagerReadEffectiveParams.cpp" class ReadEffectiveParams { public: - explicit ReadEffectiveParams(EclMaterialLawManager::InitParams& init_params); + explicit ReadEffectiveParams(EclMaterialLawManager::InitParams& init_params); void read(); private: - std::vector normalizeKrValues_(const double tolcrit, const TableColumn& krValues) const; + Storage normalizeKrValues_(const double tolcrit, const TableColumn& krValues) const; void readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx); template void readGasOilFamily2_( @@ -255,17 +297,17 @@ class EclMaterialLawManager void readGasWaterParameters_(GasWaterEffectiveParamVector& dest, unsigned satRegionIdx); void readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionIdx); - EclMaterialLawManager::InitParams& init_params_; - EclMaterialLawManager& parent_; + EclMaterialLawManager::InitParams& init_params_; + EclMaterialLawManager& parent_; const EclipseState& eclState_; }; // end of "class ReadEffectiveParams" - EclMaterialLawManager& parent_; + EclMaterialLawManager& parent_; const EclipseState& eclState_; size_t numCompressedElems_; - std::unique_ptr epsImbGridProperties_; //imbibition - std::unique_ptr epsGridProperties_; // drainage + UniquePtr epsImbGridProperties_; //imbibition + UniquePtr epsGridProperties_; // drainage }; // end of "class InitParams" @@ -277,7 +319,7 @@ class EclMaterialLawManager // Function argument 'lookupIdxOnLevelZeroAssigner' is added to lookup, for each // leaf gridview cell with index 'elemIdx', its 'lookupIdx' (index of the parent/equivalent cell on level zero). void initParamsForElements(const EclipseState& eclState, size_t numCompressedElems, - const std::function(const FieldPropsManager&, const std::string&, bool)>& + const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner, const std::function& lookupIdxOnLevelZeroAssigner); @@ -383,6 +425,7 @@ class EclMaterialLawManager OPM_TIMEFUNCTION_LOCAL(); if (!enableHysteresis()) return false; + #if !HAVE_CUDA bool changed = MaterialLaw::updateHysteresis(materialLawParams(elemIdx), fluidState); if (hasDirectionalRelperms() || hasDirectionalImbnum()) { using Dir = FaceDir::DirEnum; @@ -394,6 +437,7 @@ class EclMaterialLawManager } } return changed; + #endif // !HAVE_CUDA } void oilWaterHysteresisParams(Scalar& soMax, @@ -416,10 +460,18 @@ class EclMaterialLawManager const Scalar& somin, unsigned elemIdx); + #if !HAVE_CUDA EclEpsScalingPoints& oilWaterScaledEpsPointsDrainage(unsigned elemIdx); + #endif // !HAVE_CUDA const EclEpsScalingPointsInfo& oilWaterScaledEpsInfoDrainage(size_t elemIdx) const - { return oilWaterScaledEpsInfoDrainage_[elemIdx]; } + { + #if !HAVE_CUDA + return oilWaterScaledEpsInfoDrainage_[elemIdx]; + #else + OPM_THROW(std::runtime_error, "CUDA support is not available"); + #endif // !HAVE_CUDA + } template void serializeOp(Serializer& serializer) @@ -443,15 +495,15 @@ class EclMaterialLawManager void readGlobalThreePhaseOptions_(const Runspec& runspec); bool enableEndPointScaling_; - std::shared_ptr hysteresisConfig_; - std::vector> wagHystersisConfig_; + SharedPtr hysteresisConfig_; + Storage> wagHystersisConfig_; - - std::shared_ptr oilWaterEclEpsConfig_; - std::vector> unscaledEpsInfo_; + #if !HAVE_CUDA + SharedPtr oilWaterEclEpsConfig_; + Storage> unscaledEpsInfo_; OilWaterScalingInfoVector oilWaterScaledEpsInfoDrainage_; - std::shared_ptr gasWaterEclEpsConfig_; + SharedPtr gasWaterEclEpsConfig_; GasOilScalingPointsVector gasOilUnscaledPointsVector_; OilWaterScalingPointsVector oilWaterUnscaledPointsVector_; @@ -460,35 +512,35 @@ class EclMaterialLawManager GasOilEffectiveParamVector gasOilEffectiveParamVector_; OilWaterEffectiveParamVector oilWaterEffectiveParamVector_; GasWaterEffectiveParamVector gasWaterEffectiveParamVector_; - + #endif // !HAVE_CUDA EclMultiplexerApproach threePhaseApproach_ = EclMultiplexerApproach::Default; // this attribute only makes sense for twophase simulations! enum EclTwoPhaseApproach twoPhaseApproach_ = EclTwoPhaseApproach::GasOil; - std::vector materialLawParams_; + Storage materialLawParams_; DirectionalMaterialLawParamsPtr dirMaterialLawParams_; - std::vector satnumRegionArray_; - std::vector krnumXArray_; - std::vector krnumYArray_; - std::vector krnumZArray_; - std::vector imbnumXArray_; - std::vector imbnumYArray_; - std::vector imbnumZArray_; - std::vector imbnumRegionArray_; - std::vector stoneEtas_; + Storage satnumRegionArray_; + Storage krnumXArray_; + Storage krnumYArray_; + Storage krnumZArray_; + Storage imbnumXArray_; + Storage imbnumYArray_; + Storage imbnumZArray_; + Storage imbnumRegionArray_; + Storage stoneEtas_; bool enablePpcwmax_; - std::vector maxAllowPc_; - std::vector modifySwl_; + Storage maxAllowPc_; + Storage modifySwl_; bool hasGas; bool hasOil; bool hasWater; - std::shared_ptr gasOilConfig_; - std::shared_ptr oilWaterConfig_; - std::shared_ptr gasWaterConfig_; + SharedPtr gasOilConfig_; + SharedPtr oilWaterConfig_; + SharedPtr gasWaterConfig_; }; } // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerHystParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerHystParams.cpp index 921f78a12e9..9b1724d8192 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerHystParams.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerHystParams.cpp @@ -25,10 +25,17 @@ namespace Opm { +#if !HAVE_CUDA + /* constructors*/ -template -EclMaterialLawManager::InitParams::HystParams:: -HystParams(EclMaterialLawManager::InitParams& init_params) : +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +EclMaterialLawManager::InitParams::HystParams:: +HystParams(EclMaterialLawManager::InitParams& init_params) : init_params_{init_params}, parent_{init_params_.parent_}, eclState_{init_params_.eclState_} { @@ -39,9 +46,14 @@ HystParams(EclMaterialLawManager::InitParams& init_params) : /* public methods, alphabetically sorted */ -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: finalize() { if (hasGasOil_()) @@ -52,33 +64,53 @@ finalize() this->gasWaterParams_->finalize(); } -template -std::shared_ptr::GasOilTwoPhaseHystParams> -EclMaterialLawManager::InitParams::HystParams:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +SharedPtr::GasOilTwoPhaseHystParams> +EclMaterialLawManager::InitParams::HystParams:: getGasOilParams() { return gasOilParams_; } -template -std::shared_ptr::OilWaterTwoPhaseHystParams> -EclMaterialLawManager::InitParams::HystParams:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +SharedPtr::OilWaterTwoPhaseHystParams> +EclMaterialLawManager::InitParams::HystParams:: getOilWaterParams() { return oilWaterParams_; } -template -std::shared_ptr::GasWaterTwoPhaseHystParams> -EclMaterialLawManager::InitParams::HystParams:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +SharedPtr::GasWaterTwoPhaseHystParams> +EclMaterialLawManager::InitParams::HystParams:: getGasWaterParams() { return gasWaterParams_; } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: setConfig(unsigned satRegionIdx) { this->gasOilParams_->setConfig(this->parent_.hysteresisConfig_); @@ -93,9 +125,14 @@ setConfig(unsigned satRegionIdx) } // namespace Opm -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: setDrainageParamsGasWater(unsigned elemIdx, unsigned satRegionIdx, const std::function& lookupIdxOnLevelZeroAssigner) { @@ -112,9 +149,14 @@ setDrainageParamsGasWater(unsigned elemIdx, unsigned satRegionIdx, } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: setDrainageParamsOilGas(unsigned elemIdx, unsigned satRegionIdx, const std::function& lookupIdxOnLevelZeroAssigner) { @@ -131,9 +173,14 @@ setDrainageParamsOilGas(unsigned elemIdx, unsigned satRegionIdx, } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: setDrainageParamsOilWater(unsigned elemIdx, unsigned satRegionIdx, const std::function& lookupIdxOnLevelZeroAssigner) { @@ -159,9 +206,14 @@ setDrainageParamsOilWater(unsigned elemIdx, unsigned satRegionIdx, } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: setImbibitionParamsGasWater(unsigned elemIdx, unsigned imbRegionIdx, const std::function& lookupIdxOnLevelZeroAssigner) { @@ -181,9 +233,14 @@ setImbibitionParamsGasWater(unsigned elemIdx, unsigned imbRegionIdx, } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: setImbibitionParamsOilGas(unsigned elemIdx, unsigned imbRegionIdx, const std::function& lookupIdxOnLevelZeroAssigner) { @@ -203,9 +260,14 @@ setImbibitionParamsOilGas(unsigned elemIdx, unsigned imbRegionIdx, } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: setImbibitionParamsOilWater(unsigned elemIdx, unsigned imbRegionIdx, const std::function& lookupIdxOnLevelZeroAssigner) { @@ -227,36 +289,56 @@ setImbibitionParamsOilWater(unsigned elemIdx, unsigned imbRegionIdx, /* private methods, alphabetically sorted */ -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> bool -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: hasGasOil_() { return this->parent_.hasGas && this->parent_.hasOil; } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> bool -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: hasGasWater_() { return this->parent_.hasGas && this->parent_.hasWater && !this->parent_.hasOil; } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> bool -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: hasOilWater_() { return this->parent_.hasOil && this->parent_.hasWater; } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> std::tuple< - EclEpsScalingPointsInfo::Scalar>, - EclEpsScalingPoints::Scalar> + EclEpsScalingPointsInfo::Scalar>, + EclEpsScalingPoints::Scalar> > -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: readScaledEpsPoints_(const EclEpsGridProperties& epsGridProperties, unsigned elemIdx, EclTwoPhaseSystemType type, const std::function& fieldPropIdxOnLevelZero) { @@ -277,12 +359,17 @@ readScaledEpsPoints_(const EclEpsGridProperties& epsGridProperties, unsigned ele return {destInfo, destPoint}; } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> std::tuple< - EclEpsScalingPointsInfo::Scalar>, - EclEpsScalingPoints::Scalar> + EclEpsScalingPointsInfo::Scalar>, + EclEpsScalingPoints::Scalar> > -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: readScaledEpsPointsDrainage_(unsigned elemIdx, EclTwoPhaseSystemType type, const std::function& fieldPropIdxOnLevelZero) { @@ -290,12 +377,17 @@ readScaledEpsPointsDrainage_(unsigned elemIdx, EclTwoPhaseSystemType type, return readScaledEpsPoints_(*epsGridProperties, elemIdx, type, fieldPropIdxOnLevelZero); } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> std::tuple< - EclEpsScalingPointsInfo::Scalar>, - EclEpsScalingPoints::Scalar> + EclEpsScalingPointsInfo::Scalar>, + EclEpsScalingPoints::Scalar> > -EclMaterialLawManager::InitParams::HystParams:: +EclMaterialLawManager::InitParams::HystParams:: readScaledEpsPointsImbibition_(unsigned elemIdx, EclTwoPhaseSystemType type, const std::function& fieldPropIdxOnLevelZero) { @@ -309,5 +401,7 @@ template class EclMaterialLawManager>::Ini template class EclMaterialLawManager>::InitParams::HystParams; template class EclMaterialLawManager>::InitParams::HystParams; +#endif // !HAVE_CUDA + } // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerInitParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerInitParams.cpp index e0c4eb5d9fa..d72f612b065 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerInitParams.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerInitParams.cpp @@ -31,9 +31,14 @@ namespace Opm { /* constructors*/ -template -EclMaterialLawManager::InitParams:: -InitParams(EclMaterialLawManager& parent, const EclipseState& eclState, size_t numCompressedElems) : +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +EclMaterialLawManager::InitParams:: +InitParams(EclMaterialLawManager& parent, const EclipseState& eclState, size_t numCompressedElems) : parent_{parent}, eclState_{eclState}, numCompressedElems_{numCompressedElems} @@ -50,10 +55,15 @@ InitParams(EclMaterialLawManager& parent, const EclipseState& eclState, /* public methods */ -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams:: -run(const std::function(const FieldPropsManager&, const std::string&, bool)>& +EclMaterialLawManager::InitParams:: +run(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner, const std::function& lookupIdxOnLevelZeroAssigner) { @@ -63,9 +73,9 @@ run(const std::function(const FieldPropsManager&, const std::st copySatnumArrays_(fieldPropIntOnLeafAssigner); initOilWaterScaledEpsInfo_(); initMaterialLawParamVectors_(); - std::vector*> satnumArray; - std::vector*> imbnumArray; - std::vector*> mlpArray; + Storage*> satnumArray; + Storage*> imbnumArray; + Storage*> mlpArray; initArrays_(satnumArray, imbnumArray, mlpArray); const auto num_arrays = mlpArray.size(); for (unsigned i = 0; i < num_arrays; i++) { @@ -94,10 +104,15 @@ run(const std::function(const FieldPropsManager&, const std::st /* private methods alphabetically sorted*/ -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams:: -copySatnumArrays_(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner) +EclMaterialLawManager::InitParams:: +copySatnumArrays_(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner) { copyIntArray_(this->parent_.krnumXArray_, "KRNUMX", fieldPropIntOnLeafAssigner); copyIntArray_(this->parent_.krnumYArray_, "KRNUMY", fieldPropIntOnLeafAssigner); @@ -113,33 +128,48 @@ copySatnumArrays_(const std::function(const FieldPropsManager&, assert(!this->parent_.enableHysteresis() || this->numCompressedElems_ == this->parent_.imbnumRegionArray_.size()); } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams:: -copyIntArray_(std::vector& dest, const std::string& keyword, - const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner) +EclMaterialLawManager::InitParams:: +copyIntArray_(Storage& dest, const std::string& keyword, + const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner) { if (this->eclState_.fieldProps().has_int(keyword)) { dest = fieldPropIntOnLeafAssigner(this->eclState_.fieldProps(), keyword, /*needsTranslation*/true); } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> unsigned -EclMaterialLawManager::InitParams:: -imbRegion_(std::vector& array, unsigned elemIdx) +EclMaterialLawManager::InitParams:: +imbRegion_(Storage& array, unsigned elemIdx) { - std::vector& default_vec = this->parent_.imbnumRegionArray_; + Storage& default_vec = this->parent_.imbnumRegionArray_; return satOrImbRegion_(array, default_vec, elemIdx); } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams:: +EclMaterialLawManager::InitParams:: initArrays_( - std::vector*>& satnumArray, - std::vector*>& imbnumArray, - std::vector*>& mlpArray) + Storage*>& satnumArray, + Storage*>& imbnumArray, + Storage*>& mlpArray) { satnumArray.push_back(&this->parent_.satnumRegionArray_); imbnumArray.push_back(&this->parent_.imbnumRegionArray_); @@ -161,9 +191,14 @@ initArrays_( } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams:: +EclMaterialLawManager::InitParams:: initMaterialLawParamVectors_() { this->parent_.materialLawParams_.resize(this->numCompressedElems_); @@ -173,19 +208,29 @@ initMaterialLawParamVectors_() } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams:: +EclMaterialLawManager::InitParams:: initOilWaterScaledEpsInfo_() { // This vector will be updated in the hystParams.setDrainageOilWater() in the run() method this->parent_.oilWaterScaledEpsInfoDrainage_.resize(this->numCompressedElems_); } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams:: -initSatnumRegionArray_(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner) +EclMaterialLawManager::InitParams:: +initSatnumRegionArray_(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner) { // copy the SATNUM grid property. in some cases this is not necessary, but it // should not require much memory anyway... @@ -199,9 +244,16 @@ initSatnumRegionArray_(const std::function(const FieldPropsMana } } -template +#if !HAVE_CUDA + +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams:: +EclMaterialLawManager::InitParams:: initThreePhaseParams_(HystParams &hystParams, MaterialLawParams& materialParams, unsigned satRegionIdx, @@ -263,10 +315,82 @@ initThreePhaseParams_(HystParams &hystParams, } } // end switch() } +#else +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +void +EclMaterialLawManager::InitParams:: +initThreePhaseParams_(MaterialLawParams& materialParams, + unsigned satRegionIdx, + unsigned elemIdx) +{ + const auto& epsInfo = this->parent_.oilWaterScaledEpsInfoDrainage_[elemIdx]; -template + materialParams.setApproach(this->parent_.threePhaseApproach_); + switch (materialParams.approach()) { + case EclMultiplexerApproach::Stone1: { + auto& realParams = materialParams.template getRealParams(); + realParams.setGasOilParams(gasOilParams); + realParams.setOilWaterParams(oilWaterParams); + realParams.setSwl(epsInfo.Swl); + + if (!this->parent_.stoneEtas_.empty()) { + realParams.setEta(this->parent_.stoneEtas_[satRegionIdx]); + } + else + realParams.setEta(1.0); + realParams.finalize(); + break; + } + + case EclMultiplexerApproach::Stone2: { + auto& realParams = materialParams.template getRealParams(); + realParams.setGasOilParams(gasOilParams); + realParams.setOilWaterParams(oilWaterParams); + realParams.setSwl(epsInfo.Swl); + realParams.finalize(); + break; + } + + case EclMultiplexerApproach::Default: { + auto& realParams = materialParams.template getRealParams(); + realParams.setGasOilParams(gasOilParams); + realParams.setOilWaterParams(oilWaterParams); + realParams.setSwl(epsInfo.Swl); + realParams.finalize(); + break; + } + + case EclMultiplexerApproach::TwoPhase: { + auto& realParams = materialParams.template getRealParams(); + realParams.setGasOilParams(gasOilParams); + realParams.setOilWaterParams(oilWaterParams); + realParams.setGasWaterParams(gasWaterParams); + realParams.setApproach(this->parent_.twoPhaseApproach_); + realParams.finalize(); + break; + } + + case EclMultiplexerApproach::OnePhase: { + // Nothing to do, no parameters. + break; + } + } // end switch() +} +#endif // !HAVE_CUDA + +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams:: +EclMaterialLawManager::InitParams:: readEffectiveParameters_() { ReadEffectiveParams effectiveReader {*this}; @@ -274,9 +398,14 @@ readEffectiveParameters_() effectiveReader.read(); } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams:: +EclMaterialLawManager::InitParams:: readUnscaledEpsPointsVectors_() { if (this->parent_.hasGas && this->parent_.hasOil) { @@ -302,33 +431,48 @@ readUnscaledEpsPointsVectors_() } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> template void -EclMaterialLawManager::InitParams:: -readUnscaledEpsPoints_(Container& dest, std::shared_ptr config, EclTwoPhaseSystemType system_type) +EclMaterialLawManager::InitParams:: +readUnscaledEpsPoints_(Container& dest, SharedPtr config, EclTwoPhaseSystemType system_type) { const size_t numSatRegions = this->eclState_.runspec().tabdims().getNumSatTables(); dest.resize(numSatRegions); for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { - dest[satRegionIdx] = std::make_shared >(); + dest[satRegionIdx] = SharedPtr >(); dest[satRegionIdx]->init(this->parent_.unscaledEpsInfo_[satRegionIdx], *config, system_type); } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> unsigned -EclMaterialLawManager::InitParams:: -satRegion_(std::vector& array, unsigned elemIdx) +EclMaterialLawManager::InitParams:: +satRegion_(Storage& array, unsigned elemIdx) { - std::vector& default_vec = this->parent_.satnumRegionArray_; + Storage& default_vec = this->parent_.satnumRegionArray_; return satOrImbRegion_(array, default_vec, elemIdx); } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> unsigned -EclMaterialLawManager::InitParams:: -satOrImbRegion_(std::vector& array, std::vector& default_vec, unsigned elemIdx) +EclMaterialLawManager::InitParams:: +satOrImbRegion_(Storage& array, Storage& default_vec, unsigned elemIdx) { int value; if (array.size() > 0) { diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerReadEffectiveParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerReadEffectiveParams.cpp index 6b9418f4a77..a2e4f1f1e21 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerReadEffectiveParams.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerReadEffectiveParams.cpp @@ -37,18 +37,28 @@ namespace Opm { /* constructors*/ -template -EclMaterialLawManager::InitParams::ReadEffectiveParams:: -ReadEffectiveParams(EclMaterialLawManager::InitParams& init_params) : +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +EclMaterialLawManager::InitParams::ReadEffectiveParams:: +ReadEffectiveParams(EclMaterialLawManager::InitParams& init_params) : init_params_{init_params}, parent_{init_params_.parent_}, eclState_{init_params_.eclState_} { } /* public methods */ -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::ReadEffectiveParams:: +EclMaterialLawManager::InitParams::ReadEffectiveParams:: read() { auto& gasOilVector = this->parent_.gasOilEffectiveParamVector_; auto& oilWaterVector = this->parent_.oilWaterEffectiveParamVector_; @@ -68,9 +78,14 @@ read() { /* private methods, alphabetically sorted*/ // Relative permeability values not strictly greater than 'tolcrit' treated as zero. -template -std::vector -EclMaterialLawManager::InitParams::ReadEffectiveParams:: +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> +Storage +EclMaterialLawManager::InitParams::ReadEffectiveParams:: normalizeKrValues_(const double tolcrit, const TableColumn& krValues) const { auto kr = krValues.vectorCopy(); @@ -83,9 +98,14 @@ normalizeKrValues_(const double tolcrit, const TableColumn& krValues) const return kr; } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::ReadEffectiveParams:: +EclMaterialLawManager::InitParams::ReadEffectiveParams:: readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx) { if (!this->parent_.hasGas || !this->parent_.hasOil) @@ -115,7 +135,7 @@ readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx) readGasOilSlgof_(effParams, Swco, tolcrit, slgofTables.template getTable(satRegionIdx)); else if ( !tableManager.getSgofletTable().empty() ) { const auto& letSgofTab = tableManager.getSgofletTable()[satRegionIdx]; - const std::vector dum; // dummy arg to comform with existing interface + const Storage dum; // dummy arg to comform with existing interface effParams.setApproach(SatCurveMultiplexerApproach::LET); auto& realParams = effParams.template getRealParams(); @@ -123,7 +143,7 @@ readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx) // S=(So-Sogcr)/(1-Sogcr-Sgcr-Swco), krog = Krt*S^L/[S^L+E*(1.0-S)^T] const Scalar s_min_w = letSgofTab.s2_critical; const Scalar s_max_w = 1.0-letSgofTab.s1_critical-Swco; - const std::vector& letCoeffsOil = {s_min_w, s_max_w, + const Storage& letCoeffsOil = {s_min_w, s_max_w, static_cast(letSgofTab.l2_relperm), static_cast(letSgofTab.e2_relperm), static_cast(letSgofTab.t2_relperm), @@ -133,7 +153,7 @@ readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx) // S=(1-So-Sgcr-Swco)/(1-Sogcr-Sgcr-Swco), krg = Krt*S^L/[S^L+E*(1.0-S)^T] const Scalar s_min_nw = letSgofTab.s1_critical+Swco; const Scalar s_max_nw = 1.0-letSgofTab.s2_critical; - const std::vector& letCoeffsGas = {s_min_nw, s_max_nw, + const Storage& letCoeffsGas = {s_min_nw, s_max_nw, static_cast(letSgofTab.l1_relperm), static_cast(letSgofTab.e1_relperm), static_cast(letSgofTab.t1_relperm), @@ -141,7 +161,7 @@ readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx) realParams.setKrnSamples(letCoeffsGas, dum); // S=(So-Sorg)/(1-Sorg-Sgl-Swco), Pc = Pct + (pcir_pc-Pct)*(1-S)^L/[(1-S)^L+E*S^T] - const std::vector& letCoeffsPc = {static_cast(letSgofTab.s2_residual), + const Storage& letCoeffsPc = {static_cast(letSgofTab.s2_residual), static_cast(letSgofTab.s1_residual+Swco), static_cast(letSgofTab.l_pc), static_cast(letSgofTab.e_pc), @@ -180,10 +200,15 @@ readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx) } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> template void -EclMaterialLawManager::InitParams::ReadEffectiveParams:: +EclMaterialLawManager::InitParams::ReadEffectiveParams:: readGasOilFamily2_(GasOilEffectiveTwoPhaseParams& effParams, const Scalar Swco, const double tolcrit, @@ -192,8 +217,8 @@ readGasOilFamily2_(GasOilEffectiveTwoPhaseParams& effParams, const std::string& columnName) { // convert the saturations of the SGFN keyword from gas to oil saturations - std::vector SoSamples(sgfnTable.numRows()); - std::vector SoColumn = sofTable.getColumn("SO").vectorCopy(); + Storage SoSamples(sgfnTable.numRows()); + Storage SoColumn = sofTable.getColumn("SO").vectorCopy(); for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) { SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx); } @@ -207,16 +232,21 @@ readGasOilFamily2_(GasOilEffectiveTwoPhaseParams& effParams, realParams.finalize(); } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::ReadEffectiveParams:: +EclMaterialLawManager::InitParams::ReadEffectiveParams:: readGasOilSgof_(GasOilEffectiveTwoPhaseParams& effParams, const Scalar Swco, const double tolcrit, const SgofTable& sgofTable) { // convert the saturations of the SGOF keyword from gas to oil saturations - std::vector SoSamples(sgofTable.numRows()); + Storage SoSamples(sgofTable.numRows()); for (size_t sampleIdx = 0; sampleIdx < sgofTable.numRows(); ++ sampleIdx) { SoSamples[sampleIdx] = (1.0 - Swco) - sgofTable.get("SG", sampleIdx); } @@ -230,16 +260,21 @@ readGasOilSgof_(GasOilEffectiveTwoPhaseParams& effParams, realParams.finalize(); } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::ReadEffectiveParams:: +EclMaterialLawManager::InitParams::ReadEffectiveParams:: readGasOilSlgof_(GasOilEffectiveTwoPhaseParams& effParams, const Scalar Swco, const double tolcrit, const SlgofTable& slgofTable) { // convert the saturations of the SLGOF keyword from "liquid" to oil saturations - std::vector SoSamples(slgofTable.numRows()); + Storage SoSamples(slgofTable.numRows()); for (size_t sampleIdx = 0; sampleIdx < slgofTable.numRows(); ++ sampleIdx) { SoSamples[sampleIdx] = slgofTable.get("SL", sampleIdx) - Swco; } @@ -253,9 +288,14 @@ readGasOilSlgof_(GasOilEffectiveTwoPhaseParams& effParams, realParams.finalize(); } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::ReadEffectiveParams:: +EclMaterialLawManager::InitParams::ReadEffectiveParams:: readGasWaterParameters_(GasWaterEffectiveParamVector& dest, unsigned satRegionIdx) { if (!this->parent_.hasGas || !this->parent_.hasWater || this->parent_.hasOil) @@ -284,7 +324,7 @@ readGasWaterParameters_(GasWaterEffectiveParamVector& dest, unsigned satRegionId auto& realParams = effParams.template getRealParams(); if (!sgwfnTables.empty()){ const SgwfnTable& sgwfnTable = tableManager.getSgwfnTables().template getTable( satRegionIdx ); - std::vector SwSamples(sgwfnTable.numRows()); + Storage SwSamples(sgwfnTable.numRows()); for (size_t sampleIdx = 0; sampleIdx < sgwfnTable.numRows(); ++ sampleIdx) SwSamples[sampleIdx] = 1 - sgwfnTable.get("SG", sampleIdx); realParams.setKrwSamples(SwSamples, normalizeKrValues_(tolcrit, sgwfnTable.getColumn("KRGW"))); @@ -295,10 +335,10 @@ readGasWaterParameters_(GasWaterEffectiveParamVector& dest, unsigned satRegionId const SgfnTable& sgfnTable = tableManager.getSgfnTables().template getTable( satRegionIdx ); const SwfnTable& swfnTable = tableManager.getSwfnTables().template getTable( satRegionIdx ); - std::vector SwColumn = swfnTable.getColumn("SW").vectorCopy(); + Storage SwColumn = swfnTable.getColumn("SW").vectorCopy(); realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW"))); - std::vector SwSamples(sgfnTable.numRows()); + Storage SwSamples(sgfnTable.numRows()); for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) SwSamples[sampleIdx] = 1 - sgfnTable.get("SG", sampleIdx); realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG"))); @@ -318,10 +358,10 @@ readGasWaterParameters_(GasWaterEffectiveParamVector& dest, unsigned satRegionId effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); auto& realParams = effParams.template getRealParams(); - std::vector SwColumn = wsfTable.getColumn("SW").vectorCopy(); + Storage SwColumn = wsfTable.getColumn("SW").vectorCopy(); realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, wsfTable.getColumn("KRW"))); - std::vector SwSamples(gsfTable.numRows()); + Storage SwSamples(gsfTable.numRows()); for (size_t sampleIdx = 0; sampleIdx < gsfTable.numRows(); ++ sampleIdx) SwSamples[sampleIdx] = 1 - gsfTable.get("SG", sampleIdx); realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, gsfTable.getColumn("KRG"))); @@ -336,9 +376,14 @@ readGasWaterParameters_(GasWaterEffectiveParamVector& dest, unsigned satRegionId } } -template +template< + class Traits, + template class Storage, + template typename SharedPtr, + template typename UniquePtr +> void -EclMaterialLawManager::InitParams::ReadEffectiveParams:: +EclMaterialLawManager::InitParams::ReadEffectiveParams:: readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionIdx) { if (!this->parent_.hasOil || !this->parent_.hasWater) @@ -358,7 +403,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId { if (tableManager.hasTables("SWOF")) { const auto& swofTable = tableManager.getSwofTables().template getTable(satRegionIdx); - const std::vector SwColumn = swofTable.getColumn("SW").vectorCopy(); + const Storage SwColumn = swofTable.getColumn("SW").vectorCopy(); effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); auto& realParams = effParams.template getRealParams(); @@ -370,7 +415,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId } else if ( !tableManager.getSwofletTable().empty() ) { const auto& letTab = tableManager.getSwofletTable()[satRegionIdx]; - const std::vector dum; // dummy arg to conform with existing interface + const Storage dum; // dummy arg to conform with existing interface effParams.setApproach(SatCurveMultiplexerApproach::LET); auto& realParams = effParams.template getRealParams(); @@ -378,7 +423,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId // S=(Sw-Swcr)/(1-Sowcr-Swcr), krw = Krt*S^L/[S^L+E*(1.0-S)^T] const Scalar s_min_w = letTab.s1_critical; const Scalar s_max_w = 1.0-letTab.s2_critical; - const std::vector& letCoeffsWat = {s_min_w, s_max_w, + const Storage& letCoeffsWat = {s_min_w, s_max_w, static_cast(letTab.l1_relperm), static_cast(letTab.e1_relperm), static_cast(letTab.t1_relperm), @@ -388,7 +433,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId // S=(So-Sowcr)/(1-Sowcr-Swcr), krow = Krt*S^L/[S^L+E*(1.0-S)^T] const Scalar s_min_nw = letTab.s2_critical; const Scalar s_max_nw = 1.0-letTab.s1_critical; - const std::vector& letCoeffsOil = {s_min_nw, s_max_nw, + const Storage& letCoeffsOil = {s_min_nw, s_max_nw, static_cast(letTab.l2_relperm), static_cast(letTab.e2_relperm), static_cast(letTab.t2_relperm), @@ -396,7 +441,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId realParams.setKrnSamples(letCoeffsOil, dum); // S=(Sw-Swco)/(1-Swco-Sorw), Pc = Pct + (Pcir-Pct)*(1-S)^L/[(1-S)^L+E*S^T] - const std::vector& letCoeffsPc = {static_cast(letTab.s1_residual), + const Storage& letCoeffsPc = {static_cast(letTab.s1_residual), static_cast(letTab.s2_residual), static_cast(letTab.l_pc), static_cast(letTab.e_pc), @@ -413,7 +458,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId case SatFuncControls::KeywordFamily::Family_II: { const auto& swfnTable = tableManager.getSwfnTables().template getTable(satRegionIdx); - const std::vector SwColumn = swfnTable.getColumn("SW").vectorCopy(); + const Storage SwColumn = swfnTable.getColumn("SW").vectorCopy(); effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); auto& realParams = effParams.template getRealParams(); @@ -424,7 +469,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId if (!this->parent_.hasGas) { const auto& sof2Table = tableManager.getSof2Tables().template getTable(satRegionIdx); // convert the saturations of the SOF2 keyword from oil to water saturations - std::vector SwSamples(sof2Table.numRows()); + Storage SwSamples(sof2Table.numRows()); for (size_t sampleIdx = 0; sampleIdx < sof2Table.numRows(); ++ sampleIdx) SwSamples[sampleIdx] = 1 - sof2Table.get("SO", sampleIdx); @@ -432,7 +477,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId } else { const auto& sof3Table = tableManager.getSof3Tables().template getTable(satRegionIdx); // convert the saturations of the SOF3 keyword from oil to water saturations - std::vector SwSamples(sof3Table.numRows()); + Storage SwSamples(sof3Table.numRows()); for (size_t sampleIdx = 0; sampleIdx < sof3Table.numRows(); ++ sampleIdx) SwSamples[sampleIdx] = 1 - sof3Table.get("SO", sampleIdx); diff --git a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp index 1042fcf957e..2ecc1aaaf57 100644 --- a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp @@ -251,35 +251,52 @@ class EclTwoPhaseMaterial : public TraitsT } static Scalar trappedGasSaturation(const Params& params, bool maximumTrapping){ + #if !HAVE_CUDA + // For now we have not implemented hysteresis for CUDA, and SnTrapped is only implemented by EclHysteresisTwoPhaseLawParams if(params.approach() == EclTwoPhaseApproach::GasOil) return params.gasOilParams().SnTrapped(maximumTrapping); if(params.approach() == EclTwoPhaseApproach::GasWater) return params.gasWaterParams().SnTrapped(maximumTrapping); return 0.0; // oil-water case + #else + return 0.0; // oil-water case + #endif } static Scalar strandedGasSaturation(const Params& params, Scalar Sg, Scalar Kg){ + #if !HAVE_CUDA if(params.approach() == EclTwoPhaseApproach::GasOil) return params.gasOilParams().SnStranded(Sg, Kg); if(params.approach() == EclTwoPhaseApproach::GasWater) return params.gasWaterParams().SnStranded(Sg, Kg); return 0.0; // oil-water case + #else + return 0.0; // oil-water case + #endif } static Scalar trappedOilSaturation(const Params& params, bool maximumTrapping){ + #if !HAVE_CUDA if(params.approach() == EclTwoPhaseApproach::GasOil) return params.gasOilParams().SwTrapped(); if(params.approach() == EclTwoPhaseApproach::OilWater) return params.oilWaterParams().SnTrapped(maximumTrapping); return 0.0; // gas-water case + #else + return 0.0; // gas-water case + #endif } static Scalar trappedWaterSaturation(const Params& params){ + #if !HAVE_CUDA if(params.approach() == EclTwoPhaseApproach::GasWater) return params.gasWaterParams().SwTrapped(); if(params.approach() == EclTwoPhaseApproach::OilWater) return params.oilWaterParams().SwTrapped(); return 0.0; // gas-oil case + #else + return 0.0; // gas-oil case + #endif }