From 698f6195d32e7bd77c9d324a38090fa77af7d98b Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Fri, 17 Apr 2026 11:01:37 +0200 Subject: [PATCH 1/6] First version --- opm/material/densead/DynamicEvaluation.hpp | 1 + .../EclDefaultMaterial.hpp | 38 +++ .../EclEpsTwoPhaseLaw.hpp | 30 +++ .../EclMultiplexerMaterial.hpp | 24 ++ .../EclTwoPhaseMaterialParams.hpp | 135 +++++++--- .../SatCurveMultiplexer.hpp | 3 + .../fluidstates/BlackOilFluidState.hpp | 40 ++- .../BlackOilFluidSystem_macrotemplate.hpp | 14 +- opm/material/thermal/EclSpecrockLaw.hpp | 6 +- opm/material/thermal/EclSpecrockLawParams.hpp | 234 +++++++++++++++--- opm/material/thermal/EclThconrLaw.hpp | 30 ++- opm/material/thermal/EclThconrLawParams.hpp | 9 +- 12 files changed, 464 insertions(+), 100 deletions(-) diff --git a/opm/material/densead/DynamicEvaluation.hpp b/opm/material/densead/DynamicEvaluation.hpp index 7280ef9d8ad..19a6f669c81 100644 --- a/opm/material/densead/DynamicEvaluation.hpp +++ b/opm/material/densead/DynamicEvaluation.hpp @@ -36,6 +36,7 @@ #include #endif #include +#include #include #include diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp index f27084213f9..9cb182be53a 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp @@ -136,8 +136,13 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static void capillaryPressures(ContainerT& values, +<<<<<<< HEAD const Params& params, const FluidState& state) +======= + const Params& params, + const FluidState& state) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); using Evaluation = typename std::remove_reference::type; @@ -263,7 +268,11 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation pcgn(const Params& params, +<<<<<<< HEAD const FluidState& fs) +======= + const FluidState& fs) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); // Maximum attainable oil saturation is 1-SWL. @@ -282,7 +291,11 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation pcnw(const Params& params, +<<<<<<< HEAD const FluidState& fs) +======= + const FluidState& fs) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); const auto Sw = decay(fs.saturation(waterPhaseIdx)); @@ -347,8 +360,13 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static void relativePermeabilities(ContainerT& values, +<<<<<<< HEAD const Params& params, const FluidState& fluidState) +======= + const Params& params, + const FluidState& fluidState) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); using Evaluation = typename std::remove_reference::type; @@ -363,7 +381,11 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation krg(const Params& params, +<<<<<<< HEAD const FluidState& fluidState) +======= + const FluidState& fluidState) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); // Maximum attainable oil saturation is 1-SWL. @@ -376,7 +398,11 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation krw(const Params& params, +<<<<<<< HEAD const FluidState& fluidState) +======= + const FluidState& fluidState) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); const Evaluation sw = decay(fluidState.saturation(waterPhaseIdx)); @@ -388,7 +414,11 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation krn(const Params& params, +<<<<<<< HEAD const FluidState& fluidState) +======= + const FluidState& fluidState) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); const Scalar Swco = params.Swl(); @@ -427,7 +457,11 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation relpermOilInOilGasSystem(const Params& params, +<<<<<<< HEAD const FluidState& fluidState) +======= + const FluidState& fluidState) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); const Evaluation sw = @@ -445,7 +479,11 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation relpermOilInOilWaterSystem(const Params& params, +<<<<<<< HEAD const FluidState& fluidState) +======= + const FluidState& fluidState) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); const Evaluation sw = diff --git a/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp b/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp index d9c64c2dbe5..bd2ef93000f 100644 --- a/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp +++ b/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp @@ -374,8 +374,13 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits private: template OPM_HOST_DEVICE static Evaluation scaledToUnscaledSatTwoPoint_(const Evaluation& scaledSat, +<<<<<<< HEAD const PointsContainer& unscaledSats, const PointsContainer& scaledSats) +======= + const PointsContainer& unscaledSats, + const PointsContainer& scaledSats) +>>>>>>> 2ae970437 (First version) { return unscaledSats[0] @@ -385,8 +390,13 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits template OPM_HOST_DEVICE static Evaluation unscaledToScaledSatTwoPoint_(const Evaluation& unscaledSat, +<<<<<<< HEAD const PointsContainer& unscaledSats, const PointsContainer& scaledSats) +======= + const PointsContainer& unscaledSats, + const PointsContainer& scaledSats) +>>>>>>> 2ae970437 (First version) { return scaledSats[0] @@ -396,8 +406,13 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits template OPM_HOST_DEVICE static Evaluation scaledToUnscaledSatThreePoint_(const Evaluation& scaledSat, +<<<<<<< HEAD const PointsContainer& unscaledSats, const PointsContainer& scaledSats) +======= + const PointsContainer& unscaledSats, + const PointsContainer& scaledSats) +>>>>>>> 2ae970437 (First version) { using UnscaledSat = std::remove_cv_t>; @@ -436,8 +451,13 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits template OPM_HOST_DEVICE static Evaluation unscaledToScaledSatThreePoint_(const Evaluation& unscaledSat, +<<<<<<< HEAD const PointsContainer& unscaledSats, const PointsContainer& scaledSats) +======= + const PointsContainer& unscaledSats, + const PointsContainer& scaledSats) +>>>>>>> 2ae970437 (First version) { using ScaledSat = std::remove_cv_t>; @@ -525,8 +545,13 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits */ template OPM_HOST_DEVICE static Evaluation unscaledToScaledKrw_(const Evaluation& SwScaled, +<<<<<<< HEAD const Params& params, const Evaluation& unscaledKrw) +======= + const Params& params, + const Evaluation& unscaledKrw) +>>>>>>> 2ae970437 (First version) { const auto& cfg = params.config(); @@ -598,8 +623,13 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits */ template OPM_HOST_DEVICE static Evaluation unscaledToScaledKrn_(const Evaluation& SwScaled, +<<<<<<< HEAD const Params& params, const Evaluation& unscaledKrn) +======= + const Params& params, + const Evaluation& unscaledKrn) +>>>>>>> 2ae970437 (First version) { const auto& cfg = params.config(); diff --git a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp index 47ac800ccad..62b4a6fda34 100644 --- a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp @@ -38,6 +38,8 @@ #include +#include + #include #include @@ -45,7 +47,11 @@ namespace Opm { #if OPM_IS_INSIDE_DEVICE_FUNCTION #define OPM_ECL_MULTIPLEXER_MATERIAL_CALL(codeToCall, onePhaseCode) \ { \ +<<<<<<< HEAD constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default; \ +======= + [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default; \ +>>>>>>> 2ae970437 (First version) auto& realParams = params.template getRealParams(); \ using ActualLaw = DefaultMaterial; \ codeToCall; \ @@ -99,13 +105,21 @@ namespace Opm { #if OPM_IS_INSIDE_DEVICE_FUNCTION #define OPM_ECL_MULTIPLEXER_MATERIAL_CALL_COMPILETIME(codeToCall, onePhaseCode) \ if constexpr (Head::approach == EclMultiplexerApproach::Default) { \ +<<<<<<< HEAD constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default; \ +======= + [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default; \ +>>>>>>> 2ae970437 (First version) auto& realParams = params.template getRealParams(); \ using ActualLaw = DefaultMaterial; \ codeToCall; \ } +<<<<<<< HEAD #else +======= +#else +>>>>>>> 2ae970437 (First version) #define OPM_ECL_MULTIPLEXER_MATERIAL_CALL_COMPILETIME(codeToCall, onePhaseCode) \ if constexpr (Head::approach == EclMultiplexerApproach::Stone1) { \ [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Stone1; \ @@ -230,8 +244,13 @@ class EclMultiplexerMaterial : public TraitsT */ template OPM_HOST_DEVICE static void capillaryPressures(ContainerT& values, +<<<<<<< HEAD const Params& params, const FluidState& fluidState) +======= + const Params& params, + const FluidState& fluidState) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); if constexpr (FrontIsEclMultiplexerDispatchV) { @@ -445,8 +464,13 @@ class EclMultiplexerMaterial : public TraitsT */ template OPM_HOST_DEVICE static void relativePermeabilities(ContainerT& values, +<<<<<<< HEAD const Params& params, const FluidState& fluidState) +======= + const Params& params, + const FluidState& fluidState) +>>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); if constexpr (FrontIsEclMultiplexerDispatchV) { diff --git a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterialParams.hpp b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterialParams.hpp index 5c0255826b8..e8157a6d847 100644 --- a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterialParams.hpp +++ b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterialParams.hpp @@ -29,15 +29,29 @@ #include +#include #include -namespace Opm { +namespace Opm +{ -enum class EclTwoPhaseApproach { - GasOil, - OilWater, - GasWater -}; +enum class EclTwoPhaseApproach { GasOil, OilWater, GasWater }; + +namespace EclTwoPhaseMaterialParamsDetail +{ + + /*! + * \brief Single-argument alias for std::shared_ptr. + * + * Used as the default storage policy for EclTwoPhaseMaterialParams so the + * class can be templated on a `template class` parameter while + * retaining full backward compatibility with existing call sites that pass + * std::shared_ptr instances. + */ + template + using SharedPointer = std::shared_ptr; + +} // namespace EclTwoPhaseMaterialParamsDetail /*! * \brief Implementation for the parameters required by the material law for two-phase @@ -45,87 +59,128 @@ enum class EclTwoPhaseApproach { * * Essentially, this class just stores the two parameter objects for * the twophase capillary pressure laws. + * + * The \c StoragePolicy template parameter selects how the per-sub-law + * parameter objects are stored. By default they are heap-allocated and + * referenced through std::shared_ptr; instantiations targeting the GPU may + * pass an inline-storage policy such as \c Opm::gpuistl::ValueAsPointer so + * that the resulting object is trivially copyable to the device. */ -template +template class StoragePolicy = EclTwoPhaseMaterialParamsDetail::SharedPointer> class EclTwoPhaseMaterialParams : public EnsureFinalized { using Scalar = typename Traits::Scalar; enum { numPhases = 3 }; + public: - using EnsureFinalized :: finalize; + using EnsureFinalized::finalize; using GasOilParams = GasOilParamsT; using OilWaterParams = OilWaterParamsT; using GasWaterParams = GasWaterParamsT; + using GasOilParamsStorage = StoragePolicy; + using OilWaterParamsStorage = StoragePolicy; + using GasWaterParamsStorage = StoragePolicy; + /*! * \brief The default constructor. */ - EclTwoPhaseMaterialParams() + OPM_HOST_DEVICE EclTwoPhaseMaterialParams() = default; + + OPM_HOST_DEVICE void setApproach(EclTwoPhaseApproach newApproach) { + approach_ = newApproach; } - void setApproach(EclTwoPhaseApproach newApproach) - { approach_ = newApproach; } - - EclTwoPhaseApproach approach() const - { return approach_; } + OPM_HOST_DEVICE EclTwoPhaseApproach approach() const + { + return approach_; + } /*! * \brief The parameter object for the gas-oil twophase law. */ - const GasOilParams& gasOilParams() const - { EnsureFinalized::check(); return *gasOilParams_; } + OPM_HOST_DEVICE const GasOilParams& gasOilParams() const + { + EnsureFinalized::check(); + return *gasOilParams_; + } /*! * \brief The parameter object for the gas-oil twophase law. */ - GasOilParams& gasOilParams() - { EnsureFinalized::check(); return *gasOilParams_; } + OPM_HOST_DEVICE GasOilParams& gasOilParams() + { + EnsureFinalized::check(); + return *gasOilParams_; + } /*! * \brief Set the parameter object for the gas-oil twophase law. */ - void setGasOilParams(std::shared_ptr val) - { gasOilParams_ = val; } + OPM_HOST_DEVICE void setGasOilParams(GasOilParamsStorage val) + { + gasOilParams_ = std::move(val); + } /*! * \brief The parameter object for the oil-water twophase law. */ - const OilWaterParams& oilWaterParams() const - { EnsureFinalized::check(); return *oilWaterParams_; } + OPM_HOST_DEVICE const OilWaterParams& oilWaterParams() const + { + EnsureFinalized::check(); + return *oilWaterParams_; + } /*! * \brief The parameter object for the oil-water twophase law. */ - OilWaterParams& oilWaterParams() - { EnsureFinalized::check(); return *oilWaterParams_; } + OPM_HOST_DEVICE OilWaterParams& oilWaterParams() + { + EnsureFinalized::check(); + return *oilWaterParams_; + } /*! * \brief Set the parameter object for the oil-water twophase law. */ - void setOilWaterParams(std::shared_ptr val) - { oilWaterParams_ = val; } + OPM_HOST_DEVICE void setOilWaterParams(OilWaterParamsStorage val) + { + oilWaterParams_ = std::move(val); + } - /*! + /*! * \brief The parameter object for the gas-water twophase law. */ - const GasWaterParams& gasWaterParams() const - { EnsureFinalized::check(); return *gasWaterParams_; } + OPM_HOST_DEVICE const GasWaterParams& gasWaterParams() const + { + EnsureFinalized::check(); + return *gasWaterParams_; + } /*! * \brief The parameter object for the gas-water twophase law. */ - GasWaterParams& gasWaterParams() - { EnsureFinalized::check(); return *gasWaterParams_; } + OPM_HOST_DEVICE GasWaterParams& gasWaterParams() + { + EnsureFinalized::check(); + return *gasWaterParams_; + } /*! * \brief Set the parameter object for the gas-water twophase law. */ - void setGasWaterParams(std::shared_ptr val) - { gasWaterParams_ = val; } + OPM_HOST_DEVICE void setGasWaterParams(GasWaterParamsStorage val) + { + gasWaterParams_ = std::move(val); + } - template + template void serializeOp(Serializer& serializer) { // This is for restart serialization. @@ -135,14 +190,16 @@ class EclTwoPhaseMaterialParams : public EnsureFinalized serializer(*gasWaterParams_); } - void setSwl(Scalar) {} + void setSwl(Scalar) + { + } private: - EclTwoPhaseApproach approach_{EclTwoPhaseApproach::GasOil}; + EclTwoPhaseApproach approach_ {EclTwoPhaseApproach::GasOil}; - std::shared_ptr gasOilParams_; - std::shared_ptr oilWaterParams_; - std::shared_ptr gasWaterParams_; + GasOilParamsStorage gasOilParams_ {}; + OilWaterParamsStorage oilWaterParams_ {}; + GasWaterParamsStorage gasWaterParams_ {}; }; } // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/SatCurveMultiplexer.hpp b/opm/material/fluidmatrixinteractions/SatCurveMultiplexer.hpp index ff2b0955dc0..b48212d531b 100644 --- a/opm/material/fluidmatrixinteractions/SatCurveMultiplexer.hpp +++ b/opm/material/fluidmatrixinteractions/SatCurveMultiplexer.hpp @@ -29,7 +29,10 @@ #include "SatCurveMultiplexerParams.hpp" +<<<<<<< HEAD #include +======= +>>>>>>> 2ae970437 (First version) #include #include diff --git a/opm/material/fluidstates/BlackOilFluidState.hpp b/opm/material/fluidstates/BlackOilFluidState.hpp index 5f14f14bc30..cba54526ba7 100644 --- a/opm/material/fluidstates/BlackOilFluidState.hpp +++ b/opm/material/fluidstates/BlackOilFluidState.hpp @@ -152,7 +152,10 @@ template friend class BlackOilFluidState; +<<<<<<< HEAD +======= +>>>>>>> 2cee7df46 (First version) using FluidSystem = FluidSystemT; using ValueType = ValueT; @@ -174,6 +177,40 @@ friend class BlackOilFluidState; * * \param fluidSystem The fluid system which is used to compute various quantities */ + explicit OPM_HOST_DEVICE BlackOilFluidState(const FluidSystem* fluidSystem) + { + if constexpr (!fluidSystemIsStatic) { + fluidSystemPtr_ = fluidSystem; + } + } + + // Pointer version + template + requires std::is_pointer_v> + auto withOtherFluidSystem(OtherFluidSystemType other) const + { + using Raw = std::decay_t; + using FS = std::remove_pointer_t; + + auto bfstate = BlackOilFluidState< + ValueType, + FS, + storeTemperature, + storeEnthalpy, + enableDissolution, + enableVapwat, + enableBrine, + enableSaltPrecipitation, + enableDissolutionInWater, + enableSolvent, + NewNumStoragePhases + >(other); + + bfstate.assign(*this); + return bfstate; + } + + explicit OPM_HOST_DEVICE BlackOilFluidState(const FluidSystem& fluidSystem) { if constexpr (!fluidSystemIsStatic) { @@ -188,7 +225,8 @@ friend class BlackOilFluidState; * \return A new BlackOilFluidState object with the specified fluid system. */ template - auto withOtherFluidSystem(const OtherFluidSystemType& other) const + auto withOtherFluidSystem(typename std::enable_if, + const OtherFluidSystemType&>::type other) const { using FS = std::decay_t; diff --git a/opm/material/fluidsystems/BlackOilFluidSystem_macrotemplate.hpp b/opm/material/fluidsystems/BlackOilFluidSystem_macrotemplate.hpp index 61f10acd66a..452bb970461 100644 --- a/opm/material/fluidsystems/BlackOilFluidSystem_macrotemplate.hpp +++ b/opm/material/fluidsystems/BlackOilFluidSystem_macrotemplate.hpp @@ -179,6 +179,8 @@ class FLUIDSYSTEM_CLASSNAME : public BaseFluidSystem>&& _referenceDensity_, - Storage>&& _molarMass_, - Storage>&& _diffusionCoefficients_, + const Storage>& _referenceDensity_, + const Storage>& _molarMass_, + const Storage>& _diffusionCoefficients_, const PhaseUsageInfo& _phaseUsageInfo_, bool _isInitialized_, bool _useSaturatedTables_, @@ -210,9 +212,9 @@ class FLUIDSYSTEM_CLASSNAME : public BaseFluidSystem>(_referenceDensity_)) + , molarMass_(Storage>(_molarMass_)) + , diffusionCoefficients_(Storage>(_diffusionCoefficients_)) , phaseUsageInfo_(_phaseUsageInfo_) , isInitialized_(_isInitialized_) , useSaturatedTables_(_useSaturatedTables_) diff --git a/opm/material/thermal/EclSpecrockLaw.hpp b/opm/material/thermal/EclSpecrockLaw.hpp index dc850d68077..0a16acb02e2 100644 --- a/opm/material/thermal/EclSpecrockLaw.hpp +++ b/opm/material/thermal/EclSpecrockLaw.hpp @@ -29,6 +29,8 @@ #include "EclSpecrockLawParams.hpp" +#include + namespace Opm { @@ -51,10 +53,10 @@ class EclSpecrockLaw * \brief Given a fluid state, compute the volumetric internal energy of the rock [W/m^3]. */ template - static Evaluation solidInternalEnergy(const Params& params, const FluidState& fluidState) + OPM_HOST_DEVICE static Evaluation solidInternalEnergy(const Params& params, const FluidState& fluidState) { const auto& T = fluidState.temperature(/*phaseIdx=*/0); - return params.internalEnergyFunction().eval(T, /*extrapolate=*/true); + return params.template eval(T); } }; diff --git a/opm/material/thermal/EclSpecrockLawParams.hpp b/opm/material/thermal/EclSpecrockLawParams.hpp index 7e659a170e6..6e2f1418ff3 100644 --- a/opm/material/thermal/EclSpecrockLawParams.hpp +++ b/opm/material/thermal/EclSpecrockLawParams.hpp @@ -27,78 +27,234 @@ #ifndef OPM_ECL_SPECROCK_LAW_PARAMS_HPP #define OPM_ECL_SPECROCK_LAW_PARAMS_HPP +#include +#include +#include +#include + #include -#include #include +#include +#include +#include +#include + +namespace Opm { + +template class Storage = ::Opm::VectorWithDefaultAllocator> +class EclSpecrockLawParams; + +} // namespace Opm + +#if HAVE_CUDA +namespace Opm::gpuistl { + +template +::Opm::EclSpecrockLawParams +copy_to_gpu(const ::Opm::EclSpecrockLawParams& cpu); + +template class ContainerT> +::Opm::EclSpecrockLawParams +make_view(::Opm::EclSpecrockLawParams& gpuBuffers); + +} // namespace Opm::gpuistl +#endif // HAVE_CUDA namespace Opm { /*! * \brief The default implementation of a parameter object for the * ECL thermal law based on SPECROCK. + * + * Stores the temperature-vs-volumetric-internal-energy table in a + * templated \c Storage container so the same class can be instantiated + * as a CPU object (\c VectorWithDefaultAllocator), an owning GPU object + * (\c GpuBuffer) and a non-owning GPU view (\c GpuView) usable from a + * kernel. */ -template +template class Storage> class EclSpecrockLawParams : public EnsureFinalized { - using InternalEnergyFunction = Tabulated1DFunction; - public: using Scalar = ScalarT; + using ValueVector = Storage; - EclSpecrockLawParams(const EclSpecrockLawParams&) = default; + OPM_HOST_DEVICE EclSpecrockLawParams() = default; - EclSpecrockLawParams() - { } + OPM_HOST_DEVICE EclSpecrockLawParams(ValueVector temperatureSamples, + ValueVector internalEnergySamples) + : temperatureSamples_(std::move(temperatureSamples)) + , internalEnergySamples_(std::move(internalEnergySamples)) + { + EnsureFinalized::finalize(); + } /*! * \brief Specify the volumetric internal energy of rock via heat capacities. + * + * Available only on the CPU instantiation since GPU storage types are + * not constructible from arbitrary host containers. Integrates the + * piecewise-linear heat capacity to obtain the volumetric internal + * energy at the same temperature samples. */ - template - void setHeatCapacities(const Container& temperature, - const Container& heatCapacity) + template , + std::enable_if_t>, + int> = 0> + void setHeatCapacities(const ContainerT& temperature, + const ContainerT& heatCapacity) { - assert(temperature.size() == heatCapacity.size()); - - // integrate the heat capacity to compute the internal energy - Scalar curU = temperature[0]*heatCapacity[0]; - unsigned n = temperature.size(); - std::vector T(n); - std::vector u(n); - for (unsigned i = 0; i < temperature.size(); ++ i) { - T[i] = temperature[i]; - u[i] = curU; - - if (i >= temperature.size() - 1) + if (temperature.size() != heatCapacity.size()) { + OPM_THROW(std::invalid_argument, + "EclSpecrockLawParams: temperature and heat-capacity arrays must have " + "matching sizes"); + } + + const std::size_t n = temperature.size(); + temperatureSamples_.resize(n); + internalEnergySamples_.resize(n); + + // Integrate the heat capacity to compute the internal energy. + Scalar curU = static_cast(temperature[0]) * static_cast(heatCapacity[0]); + for (std::size_t i = 0; i < n; ++i) { + temperatureSamples_[i] = static_cast(temperature[i]); + internalEnergySamples_[i] = curU; + + if (i + 1 >= n) { break; + } - // integrate to the heat capacity from the current sampling point to the next - // one. this leads to a quadratic polynomial. - Scalar c_v0 = heatCapacity[i]; - Scalar c_v1 = heatCapacity[i + 1]; - Scalar T0 = temperature[i]; - Scalar T1 = temperature[i + 1]; - curU += 0.5*(c_v0 + c_v1)*(T1 - T0); + // Trapezoidal integration of the heat capacity from the + // current sample to the next one. + const Scalar c_v0 = static_cast(heatCapacity[i]); + const Scalar c_v1 = static_cast(heatCapacity[i + 1]); + const Scalar T0 = static_cast(temperature[i]); + const Scalar T1 = static_cast(temperature[i + 1]); + curU += Scalar(0.5) * (c_v0 + c_v1) * (T1 - T0); } - - internalEnergyFunction_.setXYContainers(T, u); } /*! - * \brief Return the function which maps temparature to the rock's volumetric - * internal energy + * \brief Set the sample tables directly. Marks the object as finalized. * - * Currently we assume this function to be piecewise linear. (Assuming piecewise - * linear heat capacity, the real function is quadratic, but the difference should be - * negligible.) + * Available only on the CPU instantiation. */ - const InternalEnergyFunction& internalEnergyFunction() const - { EnsureFinalized::check(); return internalEnergyFunction_; } + template , + std::enable_if_t>, + int> = 0> + void setSamples(const ContainerT& temperature, const ContainerT& internalEnergy) + { + if (temperature.size() != internalEnergy.size()) { + OPM_THROW(std::invalid_argument, + "EclSpecrockLawParams: temperature and internal-energy arrays must have " + "matching sizes"); + } + const std::size_t n = temperature.size(); + temperatureSamples_.resize(n); + internalEnergySamples_.resize(n); + for (std::size_t i = 0; i < n; ++i) { + temperatureSamples_[i] = static_cast(temperature[i]); + internalEnergySamples_[i] = static_cast(internalEnergy[i]); + } + EnsureFinalized::finalize(); + } + + OPM_HOST_DEVICE std::size_t numSamples() const + { + return temperatureSamples_.size(); + } + + /*! + * \brief Linearly interpolate the volumetric internal energy at a + * given temperature. The sample table is assumed sorted in + * ascending order; values outside the range are extrapolated + * linearly using the first/last segment (matching the + * previous Tabulated1DFunction::eval(T, true) behaviour. + */ + template + OPM_HOST_DEVICE Evaluation eval(const Evaluation& x) const + { + EnsureFinalized::check(); + const std::size_t n = temperatureSamples_.size(); + // n >= 2 by construction (SPECROCK tables always have >= 2 rows). + std::size_t segIdx = 0; + if (x <= temperatureSamples_[1]) { + segIdx = 0; + } else if (x >= temperatureSamples_[n - 2]) { + segIdx = n - 2; + } else { + std::size_t lo = 1; + std::size_t hi = n - 2; + while (lo + 1 < hi) { + const std::size_t mid = (lo + hi) / 2; + if (x < temperatureSamples_[mid]) { + hi = mid; + } else { + lo = mid; + } + } + segIdx = lo; + } + const Scalar x0 = temperatureSamples_[segIdx]; + const Scalar x1 = temperatureSamples_[segIdx + 1]; + const Scalar y0 = internalEnergySamples_[segIdx]; + const Scalar y1 = internalEnergySamples_[segIdx + 1]; + return y0 + (y1 - y0) * (x - x0) / (x1 - x0); + } + + OPM_HOST_DEVICE const ValueVector& temperatureSamples() const + { + EnsureFinalized::check(); + return temperatureSamples_; + } + + OPM_HOST_DEVICE const ValueVector& internalEnergySamples() const + { + EnsureFinalized::check(); + return internalEnergySamples_; + } + + ValueVector& temperatureSamplesMutable() + { + return temperatureSamples_; + } + + ValueVector& internalEnergySamplesMutable() + { + return internalEnergySamples_; + } private: - InternalEnergyFunction internalEnergyFunction_; + ValueVector temperatureSamples_ {}; + ValueVector internalEnergySamples_ {}; }; } // namespace Opm +#if HAVE_CUDA +namespace Opm::gpuistl { + +template +::Opm::EclSpecrockLawParams +copy_to_gpu(const ::Opm::EclSpecrockLawParams& cpu) +{ + return ::Opm::EclSpecrockLawParams( + GpuBuffer(cpu.temperatureSamples()), + GpuBuffer(cpu.internalEnergySamples())); +} + +template class ContainerT> +::Opm::EclSpecrockLawParams +make_view(::Opm::EclSpecrockLawParams& gpuBuffers) +{ + auto tView = make_view(gpuBuffers.temperatureSamplesMutable()); + auto eView = make_view(gpuBuffers.internalEnergySamplesMutable()); + return ::Opm::EclSpecrockLawParams(tView, eView); +} + +} // namespace Opm::gpuistl +#endif // HAVE_CUDA + #endif diff --git a/opm/material/thermal/EclThconrLaw.hpp b/opm/material/thermal/EclThconrLaw.hpp index 845e99691a1..e3802dce54d 100644 --- a/opm/material/thermal/EclThconrLaw.hpp +++ b/opm/material/thermal/EclThconrLaw.hpp @@ -29,6 +29,7 @@ #include "EclThconrLawParams.hpp" +#include #include namespace Opm @@ -52,19 +53,30 @@ class EclThconrLaw * medium. */ template - static Evaluation thermalConductivity(const Params& params, - const FluidState& fluidState) + OPM_HOST_DEVICE static Evaluation thermalConductivity(const Params& params, + const FluidState& fluidState) { // THCONR + THCONSF approach. - Scalar lambdaRef = params.referenceTotalThermalConductivity(); - static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx; - if (FluidSystem::phaseIsActive(gasPhaseIdx)) { - Scalar alpha = params.dTotalThermalConductivity_dSg(); - const Evaluation& Sg = decay(fluidState.saturation(gasPhaseIdx)); - return lambdaRef*(1.0 - alpha*Sg); + const Scalar lambdaRef = params.referenceTotalThermalConductivity(); + constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx; + + // Some fluid systems (e.g. BlackOilFluidSystemNonStatic used on the + // GPU) only expose phaseIsActive as a non-static member. Fall back + // to a fluid-state-provided fluidSystem() instance accessor in that + // case. + bool gasActive = false; + if constexpr (requires { FluidSystem::phaseIsActive(gasPhaseIdx); }) { + gasActive = FluidSystem::phaseIsActive(gasPhaseIdx); } else { - return lambdaRef; + gasActive = fluidState.fluidSystem().phaseIsActive(gasPhaseIdx); + } + + if (gasActive) { + const Scalar alpha = params.dTotalThermalConductivity_dSg(); + const Evaluation& Sg = decay(fluidState.saturation(gasPhaseIdx)); + return lambdaRef * (Scalar(1) - alpha * Sg); } + return Evaluation(lambdaRef); } }; diff --git a/opm/material/thermal/EclThconrLawParams.hpp b/opm/material/thermal/EclThconrLawParams.hpp index 93e66a1c962..a1ef483301a 100644 --- a/opm/material/thermal/EclThconrLawParams.hpp +++ b/opm/material/thermal/EclThconrLawParams.hpp @@ -27,6 +27,7 @@ #ifndef OPM_ECL_THCONR_LAW_PARAMS_HPP #define OPM_ECL_THCONR_LAW_PARAMS_HPP +#include #include namespace Opm { @@ -41,9 +42,9 @@ class EclThconrLawParams : public EnsureFinalized public: using Scalar = ScalarT; - EclThconrLawParams(const EclThconrLawParams&) = default; + OPM_HOST_DEVICE EclThconrLawParams(const EclThconrLawParams&) = default; - EclThconrLawParams() + OPM_HOST_DEVICE EclThconrLawParams() { } /*! @@ -55,7 +56,7 @@ class EclThconrLawParams : public EnsureFinalized /*! * \brief The total thermal conductivity [J/m^2 / (K/m)] of at Sg = 0 */ - Scalar referenceTotalThermalConductivity() const + OPM_HOST_DEVICE Scalar referenceTotalThermalConductivity() const { EnsureFinalized::check(); return referenceTotalThermalConductivity_; } /*! @@ -67,7 +68,7 @@ class EclThconrLawParams : public EnsureFinalized /*! * \brief The gas saturation dependence of thermal conductivity [-] */ - Scalar dTotalThermalConductivity_dSg() const + OPM_HOST_DEVICE Scalar dTotalThermalConductivity_dSg() const { EnsureFinalized::check(); return dTotalThermalConductivity_dSg_; } private: From c1f682ad457dd7f651e0e3cf2eb5e334c4d87f5d Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Wed, 6 May 2026 08:56:55 +0200 Subject: [PATCH 2/6] GPU portability: add OPM_HOST_DEVICE/OPM_THROW and refresh BlackOilFluidSystem singleton - Add OPM_HOST_DEVICE annotations and convert raw `throw` sites to OPM_THROW in Evaluation, DynamicEvaluation, BlackOilFluidState, BlackOilFluidSystem and the BrineCo2/Co2Gas/NullOil PVT classes. - Refresh BlackOilFluidSystem_macrotemplate's getNonStaticInstance() so the non-static singleton mirrors the current static state after every getInstance() call (fixes stale PVT/density tables seen by the GPU intensive-quantities path). - Add a small `mixingEnergy()` constexpr probe and an `inverseFormationVolumeFactorAndViscosity` helper used by GPU code. - Update bin/genEvalSpecializations.py so generated Evaluation*.hpp pick up ErrorMacros.hpp. CopyablePtr: GPU-friendly storage and host/device decorations Extends `CopyablePtr` so it can be used both as a plain CPU smart pointer and as a GPU-side handle (host/device-decorated accessors, allocator-aware construction, explicit raw-pointer access for kernels). The behaviour on the host is unchanged. Add missing decorator Add polymorphism check to enforce commented warning remove unneeded header remove diff improve robustness and formatting remove checks per cell add todo finish rebasing add more rebasing stuff GPU assembly support on AMD and CUDA minor fix --- CMakeLists_files.cmake | 2 + opm/material/densead/DynamicEvaluation.hpp | 1 - .../EclDefaultMaterial.hpp | 38 -- .../EclEpsTwoPhaseLaw.hpp | 30 -- .../EclMultiplexerMaterial.hpp | 22 - .../GpuEclMaterialLawManager.hpp | 461 ++++++++++++++++++ .../SatCurveMultiplexer.hpp | 3 - .../SatCurveMultiplexerParams.hpp | 1 + .../fluidstates/BlackOilFluidState.hpp | 2 +- opm/material/thermal/EclThermalLawManager.hpp | 27 + .../thermal/GpuEclThermalLawManager.hpp | 365 ++++++++++++++ 11 files changed, 857 insertions(+), 95 deletions(-) create mode 100644 opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp create mode 100644 opm/material/thermal/GpuEclThermalLawManager.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index a9cd5f4dc42..958c3c873cb 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -1355,6 +1355,7 @@ list(APPEND PUBLIC_HEADER_FILES opm/material/fluidmatrixinteractions/EclMaterialLawHystParams.hpp opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.hpp opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp + opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp opm/material/fluidmatrixinteractions/EclMaterialLawReadEffectiveParams.hpp opm/material/fluidmatrixinteractions/EclMaterialLawTwoPhaseTypes.hpp opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp @@ -1470,6 +1471,7 @@ list(APPEND PUBLIC_HEADER_FILES opm/material/thermal/EclThermalConductionLawMultiplexer.hpp opm/material/thermal/EclThermalConductionLawMultiplexerParams.hpp opm/material/thermal/EclThermalLawManager.hpp + opm/material/thermal/GpuEclThermalLawManager.hpp opm/material/thermal/EnergyModuleType.hpp opm/material/thermal/FluidThermalConductionLaw.hpp opm/material/thermal/FluidThermalConductionLawParams.hpp diff --git a/opm/material/densead/DynamicEvaluation.hpp b/opm/material/densead/DynamicEvaluation.hpp index 19a6f669c81..7280ef9d8ad 100644 --- a/opm/material/densead/DynamicEvaluation.hpp +++ b/opm/material/densead/DynamicEvaluation.hpp @@ -36,7 +36,6 @@ #include #endif #include -#include #include #include diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp index 9cb182be53a..f27084213f9 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp @@ -136,13 +136,8 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static void capillaryPressures(ContainerT& values, -<<<<<<< HEAD const Params& params, const FluidState& state) -======= - const Params& params, - const FluidState& state) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); using Evaluation = typename std::remove_reference::type; @@ -268,11 +263,7 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation pcgn(const Params& params, -<<<<<<< HEAD const FluidState& fs) -======= - const FluidState& fs) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); // Maximum attainable oil saturation is 1-SWL. @@ -291,11 +282,7 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation pcnw(const Params& params, -<<<<<<< HEAD const FluidState& fs) -======= - const FluidState& fs) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); const auto Sw = decay(fs.saturation(waterPhaseIdx)); @@ -360,13 +347,8 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static void relativePermeabilities(ContainerT& values, -<<<<<<< HEAD const Params& params, const FluidState& fluidState) -======= - const Params& params, - const FluidState& fluidState) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); using Evaluation = typename std::remove_reference::type; @@ -381,11 +363,7 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation krg(const Params& params, -<<<<<<< HEAD const FluidState& fluidState) -======= - const FluidState& fluidState) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); // Maximum attainable oil saturation is 1-SWL. @@ -398,11 +376,7 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation krw(const Params& params, -<<<<<<< HEAD const FluidState& fluidState) -======= - const FluidState& fluidState) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); const Evaluation sw = decay(fluidState.saturation(waterPhaseIdx)); @@ -414,11 +388,7 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation krn(const Params& params, -<<<<<<< HEAD const FluidState& fluidState) -======= - const FluidState& fluidState) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); const Scalar Swco = params.Swl(); @@ -457,11 +427,7 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation relpermOilInOilGasSystem(const Params& params, -<<<<<<< HEAD const FluidState& fluidState) -======= - const FluidState& fluidState) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); const Evaluation sw = @@ -479,11 +445,7 @@ class EclDefaultMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation relpermOilInOilWaterSystem(const Params& params, -<<<<<<< HEAD const FluidState& fluidState) -======= - const FluidState& fluidState) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); const Evaluation sw = diff --git a/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp b/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp index bd2ef93000f..d9c64c2dbe5 100644 --- a/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp +++ b/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp @@ -374,13 +374,8 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits private: template OPM_HOST_DEVICE static Evaluation scaledToUnscaledSatTwoPoint_(const Evaluation& scaledSat, -<<<<<<< HEAD const PointsContainer& unscaledSats, const PointsContainer& scaledSats) -======= - const PointsContainer& unscaledSats, - const PointsContainer& scaledSats) ->>>>>>> 2ae970437 (First version) { return unscaledSats[0] @@ -390,13 +385,8 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits template OPM_HOST_DEVICE static Evaluation unscaledToScaledSatTwoPoint_(const Evaluation& unscaledSat, -<<<<<<< HEAD const PointsContainer& unscaledSats, const PointsContainer& scaledSats) -======= - const PointsContainer& unscaledSats, - const PointsContainer& scaledSats) ->>>>>>> 2ae970437 (First version) { return scaledSats[0] @@ -406,13 +396,8 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits template OPM_HOST_DEVICE static Evaluation scaledToUnscaledSatThreePoint_(const Evaluation& scaledSat, -<<<<<<< HEAD const PointsContainer& unscaledSats, const PointsContainer& scaledSats) -======= - const PointsContainer& unscaledSats, - const PointsContainer& scaledSats) ->>>>>>> 2ae970437 (First version) { using UnscaledSat = std::remove_cv_t>; @@ -451,13 +436,8 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits template OPM_HOST_DEVICE static Evaluation unscaledToScaledSatThreePoint_(const Evaluation& unscaledSat, -<<<<<<< HEAD const PointsContainer& unscaledSats, const PointsContainer& scaledSats) -======= - const PointsContainer& unscaledSats, - const PointsContainer& scaledSats) ->>>>>>> 2ae970437 (First version) { using ScaledSat = std::remove_cv_t>; @@ -545,13 +525,8 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits */ template OPM_HOST_DEVICE static Evaluation unscaledToScaledKrw_(const Evaluation& SwScaled, -<<<<<<< HEAD const Params& params, const Evaluation& unscaledKrw) -======= - const Params& params, - const Evaluation& unscaledKrw) ->>>>>>> 2ae970437 (First version) { const auto& cfg = params.config(); @@ -623,13 +598,8 @@ class EclEpsTwoPhaseLaw : public EffLawT::Traits */ template OPM_HOST_DEVICE static Evaluation unscaledToScaledKrn_(const Evaluation& SwScaled, -<<<<<<< HEAD const Params& params, const Evaluation& unscaledKrn) -======= - const Params& params, - const Evaluation& unscaledKrn) ->>>>>>> 2ae970437 (First version) { const auto& cfg = params.config(); diff --git a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp index 62b4a6fda34..348e737c1f2 100644 --- a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp @@ -47,11 +47,7 @@ namespace Opm { #if OPM_IS_INSIDE_DEVICE_FUNCTION #define OPM_ECL_MULTIPLEXER_MATERIAL_CALL(codeToCall, onePhaseCode) \ { \ -<<<<<<< HEAD constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default; \ -======= - [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default; \ ->>>>>>> 2ae970437 (First version) auto& realParams = params.template getRealParams(); \ using ActualLaw = DefaultMaterial; \ codeToCall; \ @@ -105,21 +101,13 @@ namespace Opm { #if OPM_IS_INSIDE_DEVICE_FUNCTION #define OPM_ECL_MULTIPLEXER_MATERIAL_CALL_COMPILETIME(codeToCall, onePhaseCode) \ if constexpr (Head::approach == EclMultiplexerApproach::Default) { \ -<<<<<<< HEAD constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default; \ -======= - [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default; \ ->>>>>>> 2ae970437 (First version) auto& realParams = params.template getRealParams(); \ using ActualLaw = DefaultMaterial; \ codeToCall; \ } -<<<<<<< HEAD #else -======= -#else ->>>>>>> 2ae970437 (First version) #define OPM_ECL_MULTIPLEXER_MATERIAL_CALL_COMPILETIME(codeToCall, onePhaseCode) \ if constexpr (Head::approach == EclMultiplexerApproach::Stone1) { \ [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Stone1; \ @@ -244,13 +232,8 @@ class EclMultiplexerMaterial : public TraitsT */ template OPM_HOST_DEVICE static void capillaryPressures(ContainerT& values, -<<<<<<< HEAD const Params& params, const FluidState& fluidState) -======= - const Params& params, - const FluidState& fluidState) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); if constexpr (FrontIsEclMultiplexerDispatchV) { @@ -464,13 +447,8 @@ class EclMultiplexerMaterial : public TraitsT */ template OPM_HOST_DEVICE static void relativePermeabilities(ContainerT& values, -<<<<<<< HEAD const Params& params, const FluidState& fluidState) -======= - const Params& params, - const FluidState& fluidState) ->>>>>>> 2ae970437 (First version) { OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); if constexpr (FrontIsEclMultiplexerDispatchV) { diff --git a/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp new file mode 100644 index 00000000000..46611f4c688 --- /dev/null +++ b/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp @@ -0,0 +1,461 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2026 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +/*! + * \file + * \brief GPU-compatible, simplified version of Opm::EclMaterialLaw::Manager. + * + * This class only supports the actions needed by FlowProblem, + * FlowProblemBlackoil and the BlackOilIntensiveQuantities. It re-uses the + * generic two-phase material types from opm-common + * (\c Opm::EclTwoPhaseMaterial and \c Opm::EclTwoPhaseMaterialParams) with + * a value-storage policy (\c Opm::gpuistl::ValueAsPointer) so the per-cell + * material law parameters are trivially copyable to the device. + * + * The chain of multiplexers used by the CPU manager + * (\c EclMultiplexerMaterial \f$\to\f$ \c EclEpsTwoPhaseLaw \f$\to\f$ + * \c EclHysteresisTwoPhaseLaw \f$\to\f$ \c SatCurveMultiplexer) is bypassed: + * the GPU manager unwraps the CPU multiplexer parameters down to the + * underlying \c PiecewiseLinearTwoPhaseMaterialParams and uploads only the + * piecewise-linear sample tables. + */ +#ifndef OPM_GPU_ECL_MATERIAL_LAW_MANAGER_HPP +#define OPM_GPU_ECL_MATERIAL_LAW_MANAGER_HPP + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm::EclMaterialLaw +{ + +namespace detail +{ + + /*! + * \brief Walk down a CPU material-law parameter object until the enclosed + * \c PiecewiseLinearTwoPhaseMaterialParams is reached. + * + * Supports the standard layering used by the CPU \c EclMaterialLaw::Manager: + * \c EclHysteresisTwoPhaseLawParams \f$\to\f$ \c EclEpsTwoPhaseLawParams + * \f$\to\f$ \c SatCurveMultiplexerParams + * \f$\to\f$ \c PiecewiseLinearTwoPhaseMaterialParams. + */ + template + const auto& extractCpuPlParams(const CpuParams& params) + { + if constexpr (requires { params.drainageParams(); }) { + return extractCpuPlParams(params.drainageParams()); + } else if constexpr (requires { params.effectiveLawParams(); }) { + return extractCpuPlParams(params.effectiveLawParams()); + } else if constexpr (requires { + params.template getRealParams< + ::Opm::SatCurveMultiplexerApproach::PiecewiseLinear>(); + }) { + return params + .template getRealParams<::Opm::SatCurveMultiplexerApproach::PiecewiseLinear>(); + } else { + return params; + } + } + + /*! + * \brief Optional host-side base used by \c GpuManager when its storage is + * owning device memory: keeps the per-cell sample buffers alive so + * that the \c GpuView pointers stored inside the cells' + * \c PiecewiseLinearTwoPhaseMaterialParams remain valid for as long + * as the manager exists. Empty for any non-owning storage so that + * \c GpuView based managers stay device-trivially-copyable. + */ + template + struct GpuPiecewiseLinearSampleHolder { + }; + + template + struct GpuPiecewiseLinearSampleHolder { + std::vector<::Opm::gpuistl::GpuBuffer> sampleBuffers {}; + }; + +} // namespace detail + +/*! + * \brief A minimal, GPU-compatible material-law manager. + * + * Only the EclDefaultMaterial three-phase law and the EclTwoPhaseMaterial + * two-phase law are supported. The two-phase gas/oil and oil/water sub-law + * types are template parameters so the caller can choose any GPU compatible + * implementation (e.g. PiecewiseLinear). + * + * \tparam TraitsT Three-phase material traits. + * \tparam GasOilLawT Two-phase gas/oil law type (must implement the + * saturation-only API). + * \tparam OilWaterLawT Two-phase oil/water law type (must implement the + * saturation-only API). + * \tparam Storage Storage container template; defaults to a CPU vector. + * Use \c Opm::gpuistl::GpuBuffer for owning GPU storage + * and \c Opm::gpuistl::GpuView for non-owning GPU + * storage. + * \tparam MaterialLawT The three-phase material-law type to use; defaults to + * \c Opm::EclDefaultMaterial. + */ +template class Storage = ::Opm::VectorWithDefaultAllocator, + class MaterialLawT = ::Opm::EclDefaultMaterial> +class GpuManager : private detail::GpuPiecewiseLinearSampleHolder< + typename TraitsT::Scalar, + std::is_same_v, ::Opm::gpuistl::GpuBuffer>> +{ +public: + using Traits = TraitsT; + using Scalar = typename Traits::Scalar; + using GasOilLaw = GasOilLawT; + using OilWaterLaw = OilWaterLawT; + using GasOilParams = typename GasOilLaw::Params; + using OilWaterParams = typename OilWaterLaw::Params; + + /*! + * \brief The actual material-law type used by the IQ. + * + * Defaults to the three-phase \c EclDefaultMaterial. Pass an + * \c Opm::EclTwoPhaseMaterial<...> instantiation (with a value-storage + * \c EclTwoPhaseMaterialParams, e.g. one using + * \c Opm::gpuistl::ValueAsPointer) to instead use the two-phase law. + * Only the GasWater sub-approach is currently supported by the + * from-CPU constructors below; this matches the CO2STORE setup. + */ + using MaterialLaw = MaterialLawT; + using MaterialLawParams = typename MaterialLaw::Params; + + static constexpr int waterPhaseIdx = Traits::wettingPhaseIdx; + static constexpr int oilPhaseIdx = Traits::nonWettingPhaseIdx; + static constexpr int gasPhaseIdx = Traits::gasPhaseIdx; + static constexpr int numPhases = Traits::numPhases; + +private: + /*! \brief Trait: true if the material-law parameters are an + * \c EclTwoPhaseMaterialParams (i.e. expose a nested + * \c GasWaterParams type). */ + template + struct IsTwoPhaseMaterialLawParams : std::false_type { + }; + + template + struct IsTwoPhaseMaterialLawParams> + : std::true_type { + }; + + static constexpr bool isTwoPhase = IsTwoPhaseMaterialLawParams::value; + +public: + OPM_HOST_DEVICE GpuManager() = default; + + OPM_HOST_DEVICE GpuManager(Storage materialLawParams, + Storage satnumRegionArray) + : materialLawParams_(std::move(materialLawParams)) + , satnumRegionArray_(std::move(satnumRegionArray)) + { + } + + /*! + * \brief Construct from a CPU \c Opm::EclMaterialLaw::Manager directly + * into device-resident \c GpuBuffer storage. + * + * For each cell, the per-cell piecewise-linear sample arrays are + * uploaded to the GPU as individual \c GpuBuffer instances kept + * alive by this manager (via its private + * \c detail::GpuPiecewiseLinearSampleHolder base). The cell's + * \c MaterialLawParams is populated with \c GpuView views + * referencing those device buffers and is then bulk-copied to the + * device alongside the satnum array. + * + * Only enabled when this manager itself uses \c GpuBuffer storage and + * the two-phase laws use \c GpuView sample storage. + */ + template , + class GasOilParamsArg = GasOilParams, + class OilWaterParamsArg = OilWaterParams, + class IntegerStorage = Storage, + class MaterialLawParamsStorageT = Storage, + std::enable_if_t> + && std::is_same_v> + && std::is_same_v> + && std::is_same_v>, + int> + = 0> + explicit GpuManager(const CpuManager& cpu, std::size_t numElements) + : detail::GpuPiecewiseLinearSampleHolder {} + , materialLawParams_( + buildHostMaterialLawParamsForGpu(cpu, numElements, this->sampleBuffers)) + , satnumRegionArray_(buildHostSatnumRegionArray(cpu, numElements)) + { + } + + /*! \brief Material-law parameters of an active cell. */ + OPM_HOST_DEVICE MaterialLawParams materialLawParams(unsigned elemIdx) const + { + // Return by value: GpuView::operator[] const already returns by + // value, so binding a reference to materialLawParams_[idx] inside + // this function would dangle. Returning by value lets the caller + // (which typically does `const auto& mp = ...`) safely extend the + // temporary's lifetime. + return materialLawParams_[elemIdx]; + } + + OPM_HOST_DEVICE MaterialLawParams& materialLawParams(unsigned elemIdx) + { + return materialLawParams_[elemIdx]; + } + + /*! \brief Saturation-region index of an active cell. */ + OPM_HOST_DEVICE int satnumRegionIdx(unsigned elemIdx) const + { + return satnumRegionArray_[elemIdx]; + } + + /*! \brief Direct access to the underlying storages + * (used by copy_to_gpu / make_view). */ + const Storage& materialLawParamsStorage() const + { + return materialLawParams_; + } + Storage& materialLawParamsStorage() + { + return materialLawParams_; + } + + const Storage& satnumRegionArrayStorage() const + { + return satnumRegionArray_; + } + Storage& satnumRegionArrayStorage() + { + return satnumRegionArray_; + } + +private: + /*! + * \brief Build the per-cell host-side \c MaterialLawParams vector for + * the GpuBuffer-storage from-CPU constructor. Allocates one + * \c GpuBuffer per per-cell sample array and wraps each + * in a \c GpuView stored inside the cell's MLP. + * The owning device buffers are pushed into \p sampleBuffers, + * which is the manager's holder member. + */ + template + static std::vector + buildHostMaterialLawParamsForGpu(const CpuManager& cpu, + std::size_t numElements, + std::vector<::Opm::gpuistl::GpuBuffer>& sampleBuffers) + { + // The CPU manager produces one piecewise-linear sample table per + // sub-law per cell. For a three-phase EclDefaultMaterial that is up + // to twelve buffers per cell (six per gas/oil and oil/water sub-law); + // for the two-phase EclTwoPhaseMaterial it is six (a single + // gas/water sub-law). Reserving the upper bound avoids reallocation. + sampleBuffers.reserve(numElements * 12); + std::vector materialLawParams(numElements); + + auto pushSampleBuffer = [&](const auto& sampleVector) { + std::vector hostCopy(sampleVector.begin(), sampleVector.end()); + sampleBuffers.emplace_back(hostCopy); + const auto& deviceBuffer = sampleBuffers.back(); + return ::Opm::gpuistl::GpuView(deviceBuffer.data(), deviceBuffer.size()); + }; + + for (std::size_t i = 0; i < numElements; ++i) { + const auto& cpuMaterialParams = cpu.materialLawParams(static_cast(i)); + + if constexpr (isTwoPhase) { + buildTwoPhaseCellParams(cpuMaterialParams, materialLawParams[i], pushSampleBuffer); + } else { + buildThreePhaseCellParams( + cpuMaterialParams, materialLawParams[i], pushSampleBuffer); + } + } + return materialLawParams; + } + + /*! + * \brief Build a single cell's two-phase \c EclTwoPhaseMaterialParams + * instance with GPU-resident piecewise-linear sample tables. + * + * Only the \c GasWater sub-approach is currently supported, since this + * is the only configuration produced by CO2STORE-style decks. + */ + template + static void buildTwoPhaseCellParams(const CpuMaterialParams& cpuMaterialParams, + MaterialLawParams& cellParams, + PushBuffer& pushSampleBuffer) + { + if (cpuMaterialParams.approach() != ::Opm::EclMultiplexerApproach::TwoPhase) { + throw std::logic_error("Opm::EclMaterialLaw::GpuManager configured for the two-phase " + "material law requires a CPU material-law parameter set that " + "uses EclMultiplexerApproach::TwoPhase."); + } + const auto& cpuTwoPhaseParams + = cpuMaterialParams.template getRealParams<::Opm::EclMultiplexerApproach::TwoPhase>(); + if (cpuTwoPhaseParams.approach() != ::Opm::EclTwoPhaseApproach::GasWater) { + throw std::logic_error("Opm::EclMaterialLaw::GpuManager only supports the GasWater " + "two-phase sub-approach."); + } + + using GasWaterParams = typename MaterialLawParams::GasWaterParams; + const auto& cpuGasWaterPiecewiseLinear + = detail::extractCpuPlParams(cpuTwoPhaseParams.gasWaterParams()); + GasWaterParams gasWaterParams(pushSampleBuffer(cpuGasWaterPiecewiseLinear.SwPcwnSamples()), + pushSampleBuffer(cpuGasWaterPiecewiseLinear.pcwnSamples()), + pushSampleBuffer(cpuGasWaterPiecewiseLinear.SwKrwSamples()), + pushSampleBuffer(cpuGasWaterPiecewiseLinear.krwSamples()), + pushSampleBuffer(cpuGasWaterPiecewiseLinear.SwKrnSamples()), + pushSampleBuffer(cpuGasWaterPiecewiseLinear.krnSamples())); + + cellParams.setGasWaterParams( + typename MaterialLawParams::GasWaterParamsStorage(std::move(gasWaterParams))); + cellParams.setApproach(::Opm::EclTwoPhaseApproach::GasWater); + cellParams.finalize(); + } + + /*! + * \brief Build a single cell's three-phase \c EclDefaultMaterialParams + * instance with GPU-resident piecewise-linear sample tables. + */ + template + static void buildThreePhaseCellParams(const CpuMaterialParams& cpuMaterialParams, + MaterialLawParams& cellParams, + PushBuffer& pushSampleBuffer) + { + if (cpuMaterialParams.approach() != ::Opm::EclMultiplexerApproach::Default) { + throw std::logic_error("Opm::EclMaterialLaw::GpuManager only supports the Default " + "three-phase material law approach."); + } + const auto& cpuDefaultParams + = cpuMaterialParams.template getRealParams<::Opm::EclMultiplexerApproach::Default>(); + const auto& cpuGasOilPiecewiseLinear + = detail::extractCpuPlParams(cpuDefaultParams.gasOilParams()); + const auto& cpuOilWaterPiecewiseLinear + = detail::extractCpuPlParams(cpuDefaultParams.oilWaterParams()); + + GasOilParams gasOilParams(pushSampleBuffer(cpuGasOilPiecewiseLinear.SwPcwnSamples()), + pushSampleBuffer(cpuGasOilPiecewiseLinear.pcwnSamples()), + pushSampleBuffer(cpuGasOilPiecewiseLinear.SwKrwSamples()), + pushSampleBuffer(cpuGasOilPiecewiseLinear.krwSamples()), + pushSampleBuffer(cpuGasOilPiecewiseLinear.SwKrnSamples()), + pushSampleBuffer(cpuGasOilPiecewiseLinear.krnSamples())); + OilWaterParams oilWaterParams(pushSampleBuffer(cpuOilWaterPiecewiseLinear.SwPcwnSamples()), + pushSampleBuffer(cpuOilWaterPiecewiseLinear.pcwnSamples()), + pushSampleBuffer(cpuOilWaterPiecewiseLinear.SwKrwSamples()), + pushSampleBuffer(cpuOilWaterPiecewiseLinear.krwSamples()), + pushSampleBuffer(cpuOilWaterPiecewiseLinear.SwKrnSamples()), + pushSampleBuffer(cpuOilWaterPiecewiseLinear.krnSamples())); + + cellParams.setGasOilParams(std::make_shared(std::move(gasOilParams))); + cellParams.setOilWaterParams(std::make_shared(std::move(oilWaterParams))); + cellParams.setSwl(cpuDefaultParams.Swl()); + cellParams.finalize(); + } + + template + static std::vector buildHostSatnumRegionArray(const CpuManager& cpu, + std::size_t numElements) + { + std::vector satnumRegionArray(numElements); + for (std::size_t i = 0; i < numElements; ++i) { + satnumRegionArray[i] = static_cast(cpu.satnumRegionIdx(static_cast(i))); + } + return satnumRegionArray; + } + + Storage materialLawParams_ {}; + Storage satnumRegionArray_ {}; +}; + +} // namespace Opm::EclMaterialLaw + +namespace Opm::gpuistl +{ + +/*! + * \brief Copy a CPU GpuManager to GPU-resident GpuBuffer storage. + * + * The MaterialLawParams element type is assumed to be the same on the CPU + * and the GPU, i.e. the caller is responsible for making the GasOilLaw and + * OilWaterLaw GPU-compatible (typically by templating their parameter type + * on a GPU storage). + */ +template +::Opm::EclMaterialLaw::GpuManager +copy_to_gpu(const ::Opm::EclMaterialLaw::GpuManager& cpu) +{ + using GpuManagerBuffer = ::Opm::EclMaterialLaw:: + GpuManager; + using MaterialLawParams = typename GpuManagerBuffer::MaterialLawParams; + return GpuManagerBuffer(GpuBuffer(cpu.materialLawParamsStorage()), + GpuBuffer(cpu.satnumRegionArrayStorage())); +} + +/*! + * \brief Make a non-owning GpuView based GpuManager from an owning GpuBuffer + * based GpuManager. + */ +template +::Opm::EclMaterialLaw::GpuManager +make_view( + ::Opm::EclMaterialLaw::GpuManager& + buf) +{ + using GpuManagerView = ::Opm::EclMaterialLaw:: + GpuManager; + using MaterialLawParams = typename GpuManagerView::MaterialLawParams; + return GpuManagerView( + GpuView(buf.materialLawParamsStorage().data(), + buf.materialLawParamsStorage().size()), + GpuView(buf.satnumRegionArrayStorage().data(), buf.satnumRegionArrayStorage().size())); +} + +} // namespace Opm::gpuistl + +#endif // OPM_GPU_ECL_MATERIAL_LAW_MANAGER_HPP diff --git a/opm/material/fluidmatrixinteractions/SatCurveMultiplexer.hpp b/opm/material/fluidmatrixinteractions/SatCurveMultiplexer.hpp index b48212d531b..ff2b0955dc0 100644 --- a/opm/material/fluidmatrixinteractions/SatCurveMultiplexer.hpp +++ b/opm/material/fluidmatrixinteractions/SatCurveMultiplexer.hpp @@ -29,10 +29,7 @@ #include "SatCurveMultiplexerParams.hpp" -<<<<<<< HEAD #include -======= ->>>>>>> 2ae970437 (First version) #include #include diff --git a/opm/material/fluidmatrixinteractions/SatCurveMultiplexerParams.hpp b/opm/material/fluidmatrixinteractions/SatCurveMultiplexerParams.hpp index 001109a6881..c48e6da24b3 100644 --- a/opm/material/fluidmatrixinteractions/SatCurveMultiplexerParams.hpp +++ b/opm/material/fluidmatrixinteractions/SatCurveMultiplexerParams.hpp @@ -120,6 +120,7 @@ class SatCurveMultiplexerParams : public EnsureFinalized void setApproach(SatCurveMultiplexerApproach newApproach) { + // TODO: have some logic here to ensure we are choosing a approach available on GPU if we are on GPU. assert(realParams_ == 0); approach_ = newApproach; diff --git a/opm/material/fluidstates/BlackOilFluidState.hpp b/opm/material/fluidstates/BlackOilFluidState.hpp index cba54526ba7..c75236d53ce 100644 --- a/opm/material/fluidstates/BlackOilFluidState.hpp +++ b/opm/material/fluidstates/BlackOilFluidState.hpp @@ -190,7 +190,7 @@ friend class BlackOilFluidState; auto withOtherFluidSystem(OtherFluidSystemType other) const { using Raw = std::decay_t; - using FS = std::remove_pointer_t; + using FS = std::remove_cv_t>; auto bfstate = BlackOilFluidState< ValueType, diff --git a/opm/material/thermal/EclThermalLawManager.hpp b/opm/material/thermal/EclThermalLawManager.hpp index 639b46ea26e..c95b4bed261 100644 --- a/opm/material/thermal/EclThermalLawManager.hpp +++ b/opm/material/thermal/EclThermalLawManager.hpp @@ -69,6 +69,33 @@ class EclThermalLawManager const ThermalConductionLawParams& thermalConductionLawParams(unsigned elemIdx) const; + /*! \brief Return the approach used for solid energy storage. */ + EclSolidEnergyApproach solidEnergyApproach() const + { return solidEnergyApproach_; } + + /*! \brief Return the approach used for thermal conduction. */ + EclThermalConductionApproach thermalConductionApproach() const + { return thermalConductivityApproach_; } + + /*! + * \brief Return the element-index → SATNUM-region-index mapping. + * + * Only populated (non-empty) when the SPECROCK approach is used. + * Each entry is a 0-based satnum region index. + */ + const std::vector& elemToSatnumIdx() const + { return elemToSatnumIdx_; } + + /*! + * \brief Return the per-region solid-energy law parameter vector. + * + * For the SPECROCK approach this is indexed by satnum region (use + * \c elemToSatnumIdx() to map an element to its region). For the + * HEATCR approach it is indexed directly by element index. + */ + const std::vector& solidEnergyLawParamsVector() const + { return solidEnergyLawParams_; } + private: /*! * \brief Initialize the parameters for the solid energy law using using HEATCR and friends. diff --git a/opm/material/thermal/GpuEclThermalLawManager.hpp b/opm/material/thermal/GpuEclThermalLawManager.hpp new file mode 100644 index 00000000000..15700fd974a --- /dev/null +++ b/opm/material/thermal/GpuEclThermalLawManager.hpp @@ -0,0 +1,365 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright TODO ADD YEAR AND NAME OF AUTHOR + + 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 . +*/ +/*! + * \file + * + * GPU-portable, simplified version of \c Opm::EclThermalLawManager. + * + * Mirrors the design of \c Opm::EclMaterialLaw::GpuManager: + * + * - Templated on the outer storage container (per-cell / per-region + * bulk arrays of params): \c VectorWithDefaultAllocator on the CPU, + * \c GpuBuffer for owning device memory, \c GpuView for non-owning + * device storage usable from a kernel. + * - Templated on the inner storage used by the per-region SPECROCK + * sample tables (one tiny array per region). + * - The CPU \c Opm::EclThermalLawManager is treated as a *builder*: a + * dedicated constructor extracts the per-region SPECROCK sample + * tables and the per-cell THCONR coefficients and uploads them to + * device memory. + * + * Currently only the SPECROCK solid-energy approach and the THCONR + * thermal-conduction approach are supported by the builder, since these + * are the only ones used by the CO2STORE+THERMAL setup that the GPU + * dispatcher targets. The builder throws via \c OPM_THROW for any other + * approach. + */ +#ifndef OPM_GPU_ECL_THERMAL_LAW_MANAGER_HPP +#define OPM_GPU_ECL_THERMAL_LAW_MANAGER_HPP + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm::EclThermalLaw { + +namespace detail { + + /*! + * \brief Holder for the per-region SPECROCK sample buffers when the + * manager owns its device memory (i.e. the outer storage is + * \c GpuBuffer). Each per-region + * \c EclSpecrockLawParams stored in the + * bulk array references one of the buffers in this holder, so + * the holder must outlive every kernel invocation. + * + * For non-owning outer storage (host or GpuView), the holder + * is empty. + */ + template + struct SpecrockSampleHolder { + }; + + template + struct SpecrockSampleHolder { + std::vector<::Opm::gpuistl::GpuBuffer> sampleBuffers {}; + }; + +} // namespace detail + +/*! + * \brief Minimal, GPU-portable thermal-law manager. + * + * \tparam ScalarT Floating point type used for the per-cell and + * per-region sample tables. + * \tparam FluidSystemT Fluid system associated with this manager. + * \tparam OuterStorage Container used for the per-cell / per-region + * bulk arrays of params and the SATNUM index map. + * \tparam SampleStorage Container used for the per-region SPECROCK + * internal-energy / temperature sample tables. + * + * Typical instantiations: + * - CPU-only: \c GpuManager + * - GPU owning (built from CPU): \c GpuManager + * - GPU non-owning (kernel arg): \c GpuManager + */ +template class OuterStorage = ::Opm::VectorWithDefaultAllocator, + template class SampleStorage = ::Opm::VectorWithDefaultAllocator> +class GpuManager + : private detail::SpecrockSampleHolder< + ScalarT, + std::is_same_v, ::Opm::gpuistl::GpuBuffer>> +{ +public: + using Scalar = ScalarT; + using FluidSystem = FluidSystemT; + using SolidEnergyLawParams = ::Opm::EclSpecrockLawParams; + using ThermalConductionLawParams = ::Opm::EclThconrLawParams; + + static constexpr bool isOwningGpu + = std::is_same_v, ::Opm::gpuistl::GpuBuffer>; + + OPM_HOST_DEVICE GpuManager() = default; + + OPM_HOST_DEVICE + GpuManager(OuterStorage solidEnergyParams, + OuterStorage elementToSolidRegionIdx, + OuterStorage thermalConductionParams) + : solidEnergyParams_(std::move(solidEnergyParams)) + , elementToSolidRegionIdx_(std::move(elementToSolidRegionIdx)) + , thermalConductionParams_(std::move(thermalConductionParams)) + { + } + + /*! + * \brief Owning constructor that also takes ownership of the + * per-region SPECROCK sample buffers (only available for the + * owning GpuBuffer outer-storage instantiation). + */ + template = 0> + GpuManager(std::vector<::Opm::gpuistl::GpuBuffer> sampleBuffers, + OuterStorage solidEnergyParams, + OuterStorage elementToSolidRegionIdx, + OuterStorage thermalConductionParams) + : detail::SpecrockSampleHolder{std::move(sampleBuffers)} + , solidEnergyParams_(std::move(solidEnergyParams)) + , elementToSolidRegionIdx_(std::move(elementToSolidRegionIdx)) + , thermalConductionParams_(std::move(thermalConductionParams)) + { + } + + /*! \brief Solid-energy law parameters for an active cell. */ + OPM_HOST_DEVICE SolidEnergyLawParams solidEnergyLawParams(unsigned elemIdx) const + { + const int regionIdx = elementToSolidRegionIdx_[elemIdx]; + return solidEnergyParams_[regionIdx]; + } + + /*! \brief Thermal-conduction law parameters for an active cell. */ + OPM_HOST_DEVICE ThermalConductionLawParams + thermalConductionLawParams(unsigned elemIdx) const + { + return thermalConductionParams_[elemIdx]; + } + + /*! \brief Direct storage accessors used by copy_to_gpu / make_view. */ + const OuterStorage& solidEnergyParamsStorage() const + { + return solidEnergyParams_; + } + OuterStorage& solidEnergyParamsStorage() + { + return solidEnergyParams_; + } + + const OuterStorage& elementToSolidRegionIdxStorage() const + { + return elementToSolidRegionIdx_; + } + OuterStorage& elementToSolidRegionIdxStorage() + { + return elementToSolidRegionIdx_; + } + + const OuterStorage& thermalConductionParamsStorage() const + { + return thermalConductionParams_; + } + OuterStorage& thermalConductionParamsStorage() + { + return thermalConductionParams_; + } + + /*! \brief Mutable accessor for the per-region SPECROCK sample + * buffer holder; used by \c copy_to_gpu to populate the + * owning device buffers. Only available when this manager + * owns its device memory. + */ + template = 0> + std::vector<::Opm::gpuistl::GpuBuffer>& specrockSampleBuffers() + { + return this->sampleBuffers; + } + +private: + OuterStorage solidEnergyParams_ {}; + OuterStorage elementToSolidRegionIdx_ {}; + OuterStorage thermalConductionParams_ {}; +}; + +/*! + * \brief Build a CPU \c GpuManager from a CPU \c Opm::FlowProblem. + * + * Reads the solid-energy and thermal-conduction law data directly from + * the CPU problem's \c EclThermalLawManager (via \c thermalLawManager()). + * + * The per-region SPECROCK sample tables are extracted from the manager's + * per-region vector (one entry per SATNUM region), and the element-to-region + * index map is copied directly from \c EclThermalLawManager::elemToSatnumIdx(). + * The per-element THCONR thermal-conduction coefficients are then populated + * in a single, focused loop. + * + * Throws via \c OPM_THROW for any solid-energy approach other than + * SPECROCK or any thermal-conduction approach other than THCONR. + * + * \note \c CpuFlowProblemT must expose a \c thermalLawManager() accessor + * returning a (possibly const) pointer or reference to an + * \c EclThermalLawManager. + */ +template +GpuManager +buildCpuManagerFromFlowProblem(const CpuFlowProblemT& cpu, std::size_t numElements) +{ + using ManagerCpu = GpuManager; + using SolidEnergyLawParams = typename ManagerCpu::SolidEnergyLawParams; + using ThermalConductionLawParams = typename ManagerCpu::ThermalConductionLawParams; + + const auto& thermalMgr = *cpu.thermalLawManager(); + + // Validate both approaches upfront, not once per element. + if (thermalMgr.solidEnergyApproach() != ::Opm::EclSolidEnergyApproach::Specrock) { + OPM_THROW(std::logic_error, + "Opm::EclThermalLaw::GpuManager only supports the SPECROCK " + "solid-energy approach."); + } + if (thermalMgr.thermalConductionApproach() != ::Opm::EclThermalConductionApproach::Thconr) { + OPM_THROW(std::logic_error, + "Opm::EclThermalLaw::GpuManager only supports the THCONR " + "thermal-conduction approach."); + } + + // Build per-region solid-energy params by iterating over the region table + // (one entry per SATNUM region, not one per element). + const auto& cpuRegionParams = thermalMgr.solidEnergyLawParamsVector(); + std::vector hostSolidEnergyParams; + hostSolidEnergyParams.reserve(cpuRegionParams.size()); + for (const auto& cpuRegion : cpuRegionParams) { + const auto& cpuSpecrock = cpuRegion.template getRealParams< + ::Opm::EclSolidEnergyApproach::Specrock>(); + SolidEnergyLawParams params; + params.setSamples(cpuSpecrock.temperatureSamples(), + cpuSpecrock.internalEnergySamples()); + hostSolidEnergyParams.emplace_back(std::move(params)); + } + + // Copy the element→region index map directly from the CPU manager. + const auto& srcMap = thermalMgr.elemToSatnumIdx(); + std::vector hostElementToSolidRegionIdx(srcMap.begin(), srcMap.end()); + + // Build per-element THCONR thermal-conduction params. + std::vector hostThermalConductionParams(numElements); + for (std::size_t i = 0; i < numElements; ++i) { + const auto& cpuThconr = thermalMgr + .thermalConductionLawParams(static_cast(i)) + .template getRealParams<::Opm::EclThermalConductionApproach::Thconr>(); + hostThermalConductionParams[i].setReferenceTotalThermalConductivity( + cpuThconr.referenceTotalThermalConductivity()); + hostThermalConductionParams[i].setDTotalThermalConductivity_dSg( + cpuThconr.dTotalThermalConductivity_dSg()); + hostThermalConductionParams[i].finalize(); + } + + return ManagerCpu(std::move(hostSolidEnergyParams), + std::move(hostElementToSolidRegionIdx), + std::move(hostThermalConductionParams)); +} + +} // namespace Opm::EclThermalLaw + +namespace Opm::gpuistl { + +/*! + * \brief Copy a CPU \c GpuManager (plain host storage) to GPU-resident + * \c GpuBuffer storage. + * + * Each per-region SPECROCK sample table is uploaded to its own + * \c GpuBuffer; those buffers are owned by the returned + * manager's \c SpecrockSampleHolder base. The corresponding per-region + * \c EclSpecrockLawParams objects are then assembled + * on the host and uploaded as a single bulk + * \c GpuBuffer. + */ +template +::Opm::EclThermalLaw::GpuManager +copy_to_gpu(const ::Opm::EclThermalLaw::GpuManager& cpu) +{ + using ManagerBuf + = ::Opm::EclThermalLaw::GpuManager; + using SolidEnergyLawParamsView = typename ManagerBuf::SolidEnergyLawParams; + using ThermalConductionLawParams = typename ManagerBuf::ThermalConductionLawParams; + + // multiplied by two to account for both temperature and internal energy samples + std::vector> sampleBuffers; + sampleBuffers.reserve(2u * cpu.solidEnergyParamsStorage().size()); + + std::vector hostStaging; + hostStaging.reserve(cpu.solidEnergyParamsStorage().size()); + + for (const auto& cpuRegion : cpu.solidEnergyParamsStorage()) { + sampleBuffers.emplace_back(GpuBuffer(cpuRegion.temperatureSamples())); + auto& tBuf = sampleBuffers.back(); + sampleBuffers.emplace_back(GpuBuffer(cpuRegion.internalEnergySamples())); + auto& eBuf = sampleBuffers.back(); + SolidEnergyLawParamsView params(GpuView(tBuf.data(), tBuf.size()), + GpuView(eBuf.data(), eBuf.size())); + hostStaging.emplace_back(std::move(params)); + } + + return ManagerBuf(std::move(sampleBuffers), + GpuBuffer(hostStaging), + GpuBuffer(cpu.elementToSolidRegionIdxStorage()), + GpuBuffer(cpu.thermalConductionParamsStorage())); +} + +/*! + * \brief Make a non-owning \c GpuView based \c GpuManager from an owning + * \c GpuBuffer based \c GpuManager. The per-region solid-energy + * params element type is unchanged (\c GpuView sample storage in + * both cases). + */ +template +::Opm::EclThermalLaw::GpuManager +make_view(::Opm::EclThermalLaw::GpuManager& buf) +{ + using ManagerView + = ::Opm::EclThermalLaw::GpuManager; + using SolidEnergyLawParamsView = typename ManagerView::SolidEnergyLawParams; + using ThermalConductionLawParams = typename ManagerView::ThermalConductionLawParams; + return ManagerView(GpuView(buf.solidEnergyParamsStorage().data(), + buf.solidEnergyParamsStorage().size()), + GpuView(buf.elementToSolidRegionIdxStorage().data(), + buf.elementToSolidRegionIdxStorage().size()), + GpuView( + buf.thermalConductionParamsStorage().data(), + buf.thermalConductionParamsStorage().size())); +} + +} // namespace Opm::gpuistl + +#endif // OPM_GPU_ECL_THERMAL_LAW_MANAGER_HPP From 677f6fd10e621c601840219b88187f2d2f24e72d Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Thu, 28 May 2026 10:57:27 +0200 Subject: [PATCH 3/6] fix rebasing issues --- .../fluidstates/BlackOilFluidState.hpp | 45 +------------------ 1 file changed, 1 insertion(+), 44 deletions(-) diff --git a/opm/material/fluidstates/BlackOilFluidState.hpp b/opm/material/fluidstates/BlackOilFluidState.hpp index c75236d53ce..b7d53b5710e 100644 --- a/opm/material/fluidstates/BlackOilFluidState.hpp +++ b/opm/material/fluidstates/BlackOilFluidState.hpp @@ -152,10 +152,7 @@ template friend class BlackOilFluidState; -<<<<<<< HEAD -======= ->>>>>>> 2cee7df46 (First version) using FluidSystem = FluidSystemT; using ValueType = ValueT; @@ -172,45 +169,6 @@ friend class BlackOilFluidState; static constexpr bool fluidSystemIsStatic = std::is_empty_v; - /** - * \brief Construct a fluid state object. - * - * \param fluidSystem The fluid system which is used to compute various quantities - */ - explicit OPM_HOST_DEVICE BlackOilFluidState(const FluidSystem* fluidSystem) - { - if constexpr (!fluidSystemIsStatic) { - fluidSystemPtr_ = fluidSystem; - } - } - - // Pointer version - template - requires std::is_pointer_v> - auto withOtherFluidSystem(OtherFluidSystemType other) const - { - using Raw = std::decay_t; - using FS = std::remove_cv_t>; - - auto bfstate = BlackOilFluidState< - ValueType, - FS, - storeTemperature, - storeEnthalpy, - enableDissolution, - enableVapwat, - enableBrine, - enableSaltPrecipitation, - enableDissolutionInWater, - enableSolvent, - NewNumStoragePhases - >(other); - - bfstate.assign(*this); - return bfstate; - } - - explicit OPM_HOST_DEVICE BlackOilFluidState(const FluidSystem& fluidSystem) { if constexpr (!fluidSystemIsStatic) { @@ -225,8 +183,7 @@ friend class BlackOilFluidState; * \return A new BlackOilFluidState object with the specified fluid system. */ template - auto withOtherFluidSystem(typename std::enable_if, - const OtherFluidSystemType&>::type other) const + auto withOtherFluidSystem(const OtherFluidSystemType& other) const { using FS = std::decay_t; From ac82ba8aa418af33aea620bed3de5a141afd95c8 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Fri, 29 May 2026 14:06:16 +0200 Subject: [PATCH 4/6] remove duplicate include --- opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp index 348e737c1f2..47ac800ccad 100644 --- a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp @@ -38,8 +38,6 @@ #include -#include - #include #include From 2514a35a89bbbc99aec6e98b91d8a8099a133035 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Fri, 29 May 2026 14:06:24 +0200 Subject: [PATCH 5/6] remove unneeded templated struct --- .../EclTwoPhaseMaterialParams.hpp | 18 +----------------- 1 file changed, 1 insertion(+), 17 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterialParams.hpp b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterialParams.hpp index e8157a6d847..ec163c8d497 100644 --- a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterialParams.hpp +++ b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterialParams.hpp @@ -37,22 +37,6 @@ namespace Opm enum class EclTwoPhaseApproach { GasOil, OilWater, GasWater }; -namespace EclTwoPhaseMaterialParamsDetail -{ - - /*! - * \brief Single-argument alias for std::shared_ptr. - * - * Used as the default storage policy for EclTwoPhaseMaterialParams so the - * class can be templated on a `template class` parameter while - * retaining full backward compatibility with existing call sites that pass - * std::shared_ptr instances. - */ - template - using SharedPointer = std::shared_ptr; - -} // namespace EclTwoPhaseMaterialParamsDetail - /*! * \brief Implementation for the parameters required by the material law for two-phase * simulations. @@ -70,7 +54,7 @@ template class StoragePolicy = EclTwoPhaseMaterialParamsDetail::SharedPointer> + template class StoragePolicy = std::shared_ptr> class EclTwoPhaseMaterialParams : public EnsureFinalized { using Scalar = typename Traits::Scalar; From c78111084cb9e07d062664e2b25049476cc16241 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Fri, 29 May 2026 14:10:33 +0200 Subject: [PATCH 6/6] re-add removed comment --- opm/material/fluidstates/BlackOilFluidState.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/opm/material/fluidstates/BlackOilFluidState.hpp b/opm/material/fluidstates/BlackOilFluidState.hpp index b7d53b5710e..5f14f14bc30 100644 --- a/opm/material/fluidstates/BlackOilFluidState.hpp +++ b/opm/material/fluidstates/BlackOilFluidState.hpp @@ -169,6 +169,11 @@ friend class BlackOilFluidState; static constexpr bool fluidSystemIsStatic = std::is_empty_v; + /** + * \brief Construct a fluid state object. + * + * \param fluidSystem The fluid system which is used to compute various quantities + */ explicit OPM_HOST_DEVICE BlackOilFluidState(const FluidSystem& fluidSystem) { if constexpr (!fluidSystemIsStatic) {