diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 87ee07591e0..2d800d4ad81 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -1360,6 +1360,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 @@ -1475,6 +1476,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/components/iapws/Common.hpp b/opm/material/components/iapws/Common.hpp index cc6444b61fd..367b214fb8f 100644 --- a/opm/material/components/iapws/Common.hpp +++ b/opm/material/components/iapws/Common.hpp @@ -103,8 +103,8 @@ class Common OPM_HOST_DEVICE static Evaluation viscosity(const Evaluation& temperature, const Evaluation& rho) { Evaluation rhoBar = rho/322.0; - Evaluation TBar = temperature/criticalTemperature; - // Evaluation TBar = temperature/Scalar(647.096); + // Evaluation TBar = temperature/criticalTemperature; + Evaluation TBar = temperature/Scalar(647.096); // muBar = muBar_1 const Scalar Hij[6][7] = { diff --git a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp index 47ac800ccad..348e737c1f2 100644 --- a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp @@ -38,6 +38,8 @@ #include +#include + #include #include 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/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/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 c5115f15d3a..f5e0d65b18d 100644 --- a/opm/material/fluidstates/BlackOilFluidState.hpp +++ b/opm/material/fluidstates/BlackOilFluidState.hpp @@ -139,6 +139,19 @@ template +friend class BlackOilFluidState; using FluidSystem = FluidSystemT; using ValueType = ValueT; @@ -160,6 +173,40 @@ 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) { @@ -167,11 +214,37 @@ class BlackOilFluidState } } + // Non-pointer version + template + requires (!std::is_pointer_v>) + auto withOtherFluidSystem(const OtherFluidSystemType& other) const + { + using FS = std::decay_t; + + auto bfstate = BlackOilFluidState< + ValueType, + FS, + storeTemperature, + storeEnthalpy, + enableDissolution, + enableVapwat, + enableBrine, + enableSaltPrecipitation, + enableDissolutionInWater, + enableSolvent, + NewNumStoragePhases + >(other); + + bfstate.assign(*this); + return bfstate; + } + // This is intended to be used when we are converting fluid // state from a version that uses the static fluidsystem to // a version that uses a dynamic fluid system. template - auto withOtherFluidSystem(const OtherFluidSystemType& other) const + auto withOtherFluidSystem(typename std::enable_if, + const OtherFluidSystemType&>::type other) const { auto bfstate = BlackOilFluidState(fs, pvtRegionIdx)); setRsSolw(BlackOil::getRsSolw_(fs, pvtRegionIdx)); } - for (unsigned storagePhaseIdx = 0; storagePhaseIdx < numStoragePhases; ++storagePhaseIdx) { - unsigned phaseIdx = storageToCanonicalPhaseIndex_(storagePhaseIdx, fluidSystem()); - setSaturation(phaseIdx, fs.saturation(phaseIdx)); - setPressure(phaseIdx, fs.pressure(phaseIdx)); - setDensity(phaseIdx, fs.density(phaseIdx)); - + // for (unsigned storagePhaseIdx = 0; storagePhaseIdx < numStoragePhases; ++storagePhaseIdx) { + // // unsigned phaseIdx = storageToCanonicalPhaseIndex_(storagePhaseIdx, fluidSystem()); + // // setSaturation(phaseIdx, fs.saturation(phaseIdx)); + // // setPressure(phaseIdx, fs.pressure(phaseIdx)); + // // setDensity(phaseIdx, fs.density(phaseIdx)); + + // // if constexpr (storeEnthalpy) + // // setEnthalpy(phaseIdx, fs.enthalpy(phaseIdx)); + + // // setInvB(phaseIdx, getInvB_(fs, phaseIdx, pvtRegionIdx, fluidSystem())); + // } + + for (unsigned int i = 0; i < numStoragePhases; ++i) { + pressure_[i] = fs.pressure_[i]; + saturation_[i] = fs.saturation_[i]; + density_[i] = fs.density_[i]; + invB_[i] = fs.invB_[i]; if constexpr (storeEnthalpy) - setEnthalpy(phaseIdx, fs.enthalpy(phaseIdx)); - - setInvB(phaseIdx, getInvB_(fs, phaseIdx, pvtRegionIdx, fluidSystem())); + (*enthalpy_)[i] = (*fs.enthalpy_)[i]; } } 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: 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