From 6ec5b2d676763a4e0c758ad5dfa0ee41449ab16e Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Fri, 17 Apr 2026 11:01:37 +0200 Subject: [PATCH 01/10] First version 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 fix rebasing issues remove duplicate include remove unneeded templated struct re-add removed comment make the gpueclmateriallawmanager store tables per region avoid simulators includes remove unneeded gpu decorators --- CMakeLists_files.cmake | 2 + .../GpuEclMaterialLawManager.hpp | 563 ++++++++++++++++++ .../SatCurveMultiplexerParams.hpp | 1 + .../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/EclThermalLawManager.hpp | 27 + .../thermal/GpuEclThermalLawManager.hpp | 381 ++++++++++++ 9 files changed, 1202 insertions(+), 56 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 e1c084785c8..85eda02bfe9 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -1372,6 +1372,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 @@ -1487,6 +1488,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/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp new file mode 100644 index 00000000000..4469c7a0527 --- /dev/null +++ b/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp @@ -0,0 +1,563 @@ +// -*- 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::gpuistl +{ +template +class GpuBuffer; + +template +class GpuView; +} // namespace Opm::gpuistl + +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 { + }; + +#if HAVE_CUDA + template + struct GpuPiecewiseLinearSampleHolder { + std::vector<::Opm::gpuistl::GpuBuffer> sampleBuffers {}; + }; +#endif + +} // 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>> +{ +private: + #if HAVE_CUDA + struct GpuData { + std::vector<::Opm::gpuistl::GpuBuffer> sampleBuffers; + std::vector materialLawParamsByRegion; + std::vector satnumRegionArray; + }; + #endif + +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 bool isOwningGpu + = std::is_same_v, ::Opm::gpuistl::GpuBuffer>; + + 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: + GpuManager() = default; + + GpuManager(Storage materialLawParamsByRegion, + Storage satnumRegionArray) + : materialLawParamsByRegion_(std::move(materialLawParamsByRegion)) + , satnumRegionArray_(std::move(satnumRegionArray)) + { + } + + /*! + * \brief Construct from a CPU \c Opm::EclMaterialLaw::Manager directly + * into device-resident \c GpuBuffer storage. + * + * For each SATNUM region, the 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. + */ + /*! \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 materialLawParamsByRegion_[satnumRegionArray_[elemIdx]]; + } + + OPM_HOST_DEVICE MaterialLawParams& materialLawParams(unsigned elemIdx) + { + return materialLawParamsByRegion_[satnumRegionArray_[elemIdx]]; + } + + /*! \brief Saturation-region index of an active cell. */ + OPM_HOST_DEVICE int satnumRegionIdx(unsigned elemIdx) const + { + return satnumRegionArray_[elemIdx]; + } + + std::size_t materialLawParamsRegionCount() const + { + return materialLawParamsByRegion_.size(); + } + + /*! \brief Direct access to the underlying storages + * (used by copy_to_gpu / make_view). */ + const Storage& materialLawParamsStorage() const + { + return materialLawParamsByRegion_; + } + Storage& materialLawParamsStorage() + { + return materialLawParamsByRegion_; + } + + const Storage& satnumRegionArrayStorage() const + { + return satnumRegionArray_; + } + Storage& satnumRegionArrayStorage() + { + return satnumRegionArray_; + } + +#if HAVE_CUDA + template = 0> + GpuManager(GpuData data) + : detail::GpuPiecewiseLinearSampleHolder{std::move(data.sampleBuffers)} + , materialLawParamsByRegion_(::Opm::gpuistl::GpuBuffer( + data.materialLawParamsByRegion)) + , satnumRegionArray_(::Opm::gpuistl::GpuBuffer(data.satnumRegionArray)) + { + } +#endif + + template + static bool samePiecewiseLinearParams(const CpuParams& lhs, const CpuParams& rhs) + { + return lhs.SwPcwnSamples() == rhs.SwPcwnSamples() + && lhs.pcwnSamples() == rhs.pcwnSamples() + && lhs.SwKrwSamples() == rhs.SwKrwSamples() + && lhs.krwSamples() == rhs.krwSamples() + && lhs.SwKrnSamples() == rhs.SwKrnSamples() + && lhs.krnSamples() == rhs.krnSamples(); + } + + template + static bool sameSupportedMaterialLawParams(const CpuMaterialParams& lhs, + const CpuMaterialParams& rhs) + { + if (lhs.approach() != rhs.approach()) { + return false; + } + + if constexpr (isTwoPhase) { + if (lhs.approach() != ::Opm::EclMultiplexerApproach::TwoPhase) { + return false; + } + const auto& lhsTwoPhase + = lhs.template getRealParams<::Opm::EclMultiplexerApproach::TwoPhase>(); + const auto& rhsTwoPhase + = rhs.template getRealParams<::Opm::EclMultiplexerApproach::TwoPhase>(); + if (lhsTwoPhase.approach() != rhsTwoPhase.approach() + || lhsTwoPhase.approach() != ::Opm::EclTwoPhaseApproach::GasWater) { + return false; + } + return samePiecewiseLinearParams( + detail::extractCpuPlParams(lhsTwoPhase.gasWaterParams()), + detail::extractCpuPlParams(rhsTwoPhase.gasWaterParams())); + } + + if (lhs.approach() != ::Opm::EclMultiplexerApproach::Default) { + return false; + } + const auto& lhsDefault + = lhs.template getRealParams<::Opm::EclMultiplexerApproach::Default>(); + const auto& rhsDefault + = rhs.template getRealParams<::Opm::EclMultiplexerApproach::Default>(); + return lhsDefault.Swl() == rhsDefault.Swl() + && samePiecewiseLinearParams( + detail::extractCpuPlParams(lhsDefault.gasOilParams()), + detail::extractCpuPlParams(rhsDefault.gasOilParams())) + && samePiecewiseLinearParams( + detail::extractCpuPlParams(lhsDefault.oilWaterParams()), + detail::extractCpuPlParams(rhsDefault.oilWaterParams())); + } + +#if HAVE_CUDA + template + static GpuData buildGpuData(const CpuManager& cpu, std::size_t numElements) + { + GpuData data; + data.satnumRegionArray = buildHostSatnumRegionArray(cpu, numElements); + + const auto maxRegion = data.satnumRegionArray.empty() + ? -1 + : *std::max_element(data.satnumRegionArray.begin(), data.satnumRegionArray.end()); + if (maxRegion < 0) { + return data; + } + data.materialLawParamsByRegion.resize(static_cast(maxRegion) + 1u); + std::vector initialized(data.materialLawParamsByRegion.size(), false); + std::vector firstElement(data.materialLawParamsByRegion.size(), 0u); + data.sampleBuffers.reserve(data.materialLawParamsByRegion.size() * 12u); + + auto pushSampleBuffer = [&](const auto& sampleVector) { + std::vector hostCopy(sampleVector.begin(), sampleVector.end()); + data.sampleBuffers.emplace_back(hostCopy); + const auto& deviceBuffer = data.sampleBuffers.back(); + return ::Opm::gpuistl::GpuView(deviceBuffer.data(), deviceBuffer.size()); + }; + + for (std::size_t elemIdx = 0; elemIdx < numElements; ++elemIdx) { + const int region = data.satnumRegionArray[elemIdx]; + if (region < 0 || static_cast(region) >= initialized.size()) { + OPM_THROW(std::logic_error, "Invalid SATNUM region index in GPU material manager"); + } + + const auto& cpuMaterialParams = cpu.materialLawParams(static_cast(elemIdx)); + if (!initialized[region]) { + if constexpr (isTwoPhase) { + buildTwoPhaseCellParams(cpuMaterialParams, + data.materialLawParamsByRegion[region], + pushSampleBuffer); + } else { + buildThreePhaseCellParams(cpuMaterialParams, + data.materialLawParamsByRegion[region], + pushSampleBuffer); + } + initialized[region] = true; + firstElement[region] = elemIdx; + } else { + const auto& firstParams + = cpu.materialLawParams(static_cast(firstElement[region])); + if (!sameSupportedMaterialLawParams(cpuMaterialParams, firstParams)) { + OPM_THROW(std::logic_error, + "GPU material-law tables vary within one SATNUM region"); + } + } + } + return data; + } + + 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) + : GpuManager(buildGpuData(cpu, numElements)) + { + } +#endif + + /*! + * \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) { + OPM_THROW(std::logic_error, + "GPU material manager requires EclMultiplexerApproach::TwoPhase"); + } + const auto& cpuTwoPhaseParams + = cpuMaterialParams.template getRealParams<::Opm::EclMultiplexerApproach::TwoPhase>(); + if (cpuTwoPhaseParams.approach() != ::Opm::EclTwoPhaseApproach::GasWater) { + OPM_THROW(std::logic_error, + "GPU material manager only supports the GasWater two-phase 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) { + OPM_THROW(std::logic_error, + "GPU material manager only supports the Default three-phase 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 materialLawParamsByRegion_ {}; + Storage satnumRegionArray_ {}; +}; + +} // namespace Opm::EclMaterialLaw + +namespace Opm::gpuistl +{ + +#if HAVE_CUDA + +/*! + * \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; + return GpuManagerBuffer(cpu, cpu.satnumRegionArrayStorage().size()); +} + +/*! + * \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 // HAVE_CUDA + +#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/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/EclThermalLawManager.hpp b/opm/material/thermal/EclThermalLawManager.hpp index 9c977bee7f3..9a3e7364923 100644 --- a/opm/material/thermal/EclThermalLawManager.hpp +++ b/opm/material/thermal/EclThermalLawManager.hpp @@ -70,6 +70,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..97c15d61b40 --- /dev/null +++ b/opm/material/thermal/GpuEclThermalLawManager.hpp @@ -0,0 +1,381 @@ +// -*- 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 + +namespace Opm::gpuistl +{ +template +class GpuBuffer; + +template +class GpuView; +} // namespace Opm::gpuistl + +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 { + }; + +#if HAVE_CUDA + template + struct SpecrockSampleHolder { + std::vector<::Opm::gpuistl::GpuBuffer> sampleBuffers {}; + }; +#endif + +} // 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>; + + GpuManager() = default; + + 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). + */ +#if HAVE_CUDA + 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)) + { + } +#endif + + /*! \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. + */ +#if HAVE_CUDA + template = 0> + std::vector<::Opm::gpuistl::GpuBuffer>& specrockSampleBuffers() + { + return this->sampleBuffers; + } +#endif + +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 { + +#if HAVE_CUDA + +/*! + * \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 // HAVE_CUDA + +#endif // OPM_GPU_ECL_THERMAL_LAW_MANAGER_HPP From 8c1ec0ce75d1cc3fdbe662c1817bccf76805b83f Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Thu, 6 Aug 2026 11:07:20 +0200 Subject: [PATCH 02/10] use gpu friendly error-macros --- .../PiecewiseLinearTwoPhaseMaterial.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp index 1a69e4c182c..e2be873630f 100644 --- a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp @@ -29,6 +29,7 @@ #include "PiecewiseLinearTwoPhaseMaterialParams.hpp" +#include #include #include @@ -113,7 +114,7 @@ class PiecewiseLinearTwoPhaseMaterial : public TraitsT */ template OPM_HOST_DEVICE static void saturations(Container& /* values */, const Params& /* params */, const FluidState& /* fs */) - { throw std::logic_error("Not implemented: saturations()"); } + { OPM_THROW(std::logic_error, "Not implemented: saturations()"); } /*! * \brief The relative permeabilities @@ -160,11 +161,11 @@ class PiecewiseLinearTwoPhaseMaterial : public TraitsT */ template OPM_HOST_DEVICE static Evaluation Sw(const Params& /* params */, const FluidState& /* fs */) - { throw std::logic_error("Not implemented: Sw()"); } + { OPM_THROW(std::logic_error, "Not implemented: Sw()"); } template OPM_HOST_DEVICE static Evaluation twoPhaseSatSw(const Params& /* params */, const Evaluation& /* pC */) - { throw std::logic_error("Not implemented: twoPhaseSatSw()"); } + { OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatSw()"); } /*! * \brief Calculate the non-wetting phase saturations depending on From e7cc185ee2664229b6ef1bebee965e5bd1768e4a Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Thu, 6 Aug 2026 13:00:24 +0200 Subject: [PATCH 03/10] refactor speclawparams and tabulated1d --- opm/material/common/Tabulated1DFunction.cpp | 12 +- opm/material/common/Tabulated1DFunction.hpp | 128 +++++++++++++-- .../SatCurveMultiplexerParams.hpp | 1 - .../BlackOilFluidSystem_macrotemplate.hpp | 1 - opm/material/thermal/EclSpecrockLaw.hpp | 3 +- opm/material/thermal/EclSpecrockLawParams.hpp | 154 ++++++------------ .../thermal/GpuEclThermalLawManager.hpp | 25 ++- 7 files changed, 186 insertions(+), 138 deletions(-) diff --git a/opm/material/common/Tabulated1DFunction.cpp b/opm/material/common/Tabulated1DFunction.cpp index 17140d713e6..44ef7b541a3 100644 --- a/opm/material/common/Tabulated1DFunction.cpp +++ b/opm/material/common/Tabulated1DFunction.cpp @@ -28,8 +28,8 @@ namespace Opm { -template -void Tabulated1DFunction:: +template class Storage> +void Tabulated1DFunction:: printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os) const { Scalar x0 = std::min(xi0, xi1); @@ -63,10 +63,10 @@ printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os) const } template void -Tabulated1DFunction::printCSV(double,double, - unsigned,std::ostream&) const; +Tabulated1DFunction::printCSV( + double, double, unsigned, std::ostream&) const; template void -Tabulated1DFunction::printCSV(float,float, - unsigned,std::ostream&) const; +Tabulated1DFunction::printCSV( + float, float, unsigned, std::ostream&) const; } // namespace Opm diff --git a/opm/material/common/Tabulated1DFunction.hpp b/opm/material/common/Tabulated1DFunction.hpp index 7d9521e11c9..1fea415a8f3 100644 --- a/opm/material/common/Tabulated1DFunction.hpp +++ b/opm/material/common/Tabulated1DFunction.hpp @@ -27,7 +27,11 @@ #ifndef OPM_TABULATED_1D_FUNCTION_HPP #define OPM_TABULATED_1D_FUNCTION_HPP +#include #include +#include +#include +#include #include #include @@ -35,6 +39,7 @@ #include #include #include +#include #include namespace Opm { @@ -43,20 +48,45 @@ struct SegmentIndex { std::size_t value; }; +template class Storage = VectorWithDefaultAllocator> +class Tabulated1DFunction; + +#if HAVE_CUDA +namespace gpuistl { + +template +Tabulated1DFunction +copy_to_gpu(const Tabulated1DFunction& cpu); + +template class ContainerT> +Tabulated1DFunction +make_view(Tabulated1DFunction& gpuBuffers); + +} // namespace gpuistl +#endif // HAVE_CUDA + /*! * \brief Implements a linearly interpolated scalar function that depends on one * variable. */ -template +template class Storage> class Tabulated1DFunction { public: + using ValueVector = Storage; + /*! * \brief Default constructor for a piecewise linear function. * * To specfiy the acutal curve, use one of the set() methods. */ - Tabulated1DFunction() + OPM_HOST_DEVICE Tabulated1DFunction() = default; + + OPM_HOST_DEVICE Tabulated1DFunction(ValueVector xValues, + ValueVector yValues) + : xValues_(std::move(xValues)) + , yValues_(std::move(yValues)) {} /*! @@ -234,10 +264,10 @@ class Tabulated1DFunction Scalar xAt(std::size_t i) const { return xValues_[i]; } - const std::vector& xValues() const + const ValueVector& xValues() const { return xValues_; } - const std::vector& yValues() const + const ValueVector& yValues() const { return yValues_; } /*! @@ -263,14 +293,19 @@ class Tabulated1DFunction * failed assertation. */ template - Evaluation eval(const Evaluation& x, bool extrapolate = false) const + OPM_HOST_DEVICE Evaluation eval(const Evaluation& x, bool extrapolate = false) const { - SegmentIndex segIdx = findSegmentIndex(x, extrapolate); - return eval(x, segIdx); + if constexpr (OPM_IS_INSIDE_DEVICE_FUNCTION) { + return evalDevice_(x); + } + else { + SegmentIndex segIdx = findSegmentIndex(x, extrapolate); + return eval(x, segIdx); + } } template - Evaluation eval(const Evaluation& x, SegmentIndex segIdxIn) const + OPM_HOST_DEVICE Evaluation eval(const Evaluation& x, SegmentIndex segIdxIn) const { std::size_t segIdx = segIdxIn.value; Scalar x0 = xValues_[segIdx]; @@ -430,7 +465,7 @@ class Tabulated1DFunction */ void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os) const; - bool operator==(const Tabulated1DFunction& data) const { + bool operator==(const Tabulated1DFunction& data) const { return xValues_ == data.xValues_ && yValues_ == data.yValues_; } @@ -505,6 +540,37 @@ class Tabulated1DFunction } private: + template + OPM_HOST_DEVICE Evaluation evalDevice_(const Evaluation& x) const + { + const std::size_t n = xValues_.size(); + assert(n >= 2); + + std::size_t segIdx = 0; + if (x <= xValues_[1]) { + segIdx = 0; + } + else if (x >= xValues_[n - 2]) { + segIdx = n - 2; + } + else { + std::size_t lowerIdx = 1; + std::size_t upperIdx = n - 2; + while (lowerIdx + 1 < upperIdx) { + const std::size_t pivotIdx = (lowerIdx + upperIdx) / 2; + if (x < xValues_[pivotIdx]) { + upperIdx = pivotIdx; + } + else { + lowerIdx = pivotIdx; + } + } + segIdx = lowerIdx; + } + + return eval(x, SegmentIndex{segIdx}); + } + template Evaluation evalDerivative_(const Evaluation& x, std::size_t segIdx) const { @@ -555,14 +621,14 @@ class Tabulated1DFunction */ struct ComparatorX_ { - explicit ComparatorX_(const std::vector& x) + explicit ComparatorX_(const ValueVector& x) : x_(x) {} bool operator ()(std::size_t idxA, std::size_t idxB) const { return x_.at(idxA) < x_.at(idxB); } - const std::vector& x_; + const ValueVector& x_; }; /*! @@ -583,7 +649,7 @@ class Tabulated1DFunction std::ranges::sort(idxVector, cmp); // reorder the sample points - std::vector tmpX(n), tmpY(n); + ValueVector tmpX(n), tmpY(n); for (std::size_t i = 0; i < idxVector.size(); ++ i) { tmpX[i] = xValues_[idxVector[i]]; tmpY[i] = yValues_[idxVector[i]]; @@ -615,10 +681,44 @@ class Tabulated1DFunction yValues_.resize(nSamples); } - std::vector xValues_; - std::vector yValues_; +#if HAVE_CUDA + template + friend Tabulated1DFunction + gpuistl::copy_to_gpu(const Tabulated1DFunction& cpu); + + template class ContainerT> + friend Tabulated1DFunction + gpuistl::make_view(Tabulated1DFunction& gpuBuffers); +#endif + + ValueVector xValues_; + ValueVector yValues_; }; } // namespace Opm +#if HAVE_CUDA +namespace Opm::gpuistl { + +template +Tabulated1DFunction +copy_to_gpu(const Tabulated1DFunction& cpu) +{ + return Tabulated1DFunction( + GpuBuffer(cpu.xValues()), + GpuBuffer(cpu.yValues())); +} + +template class ContainerT> +Tabulated1DFunction +make_view(Tabulated1DFunction& gpuBuffers) +{ + return Tabulated1DFunction( + make_view(gpuBuffers.xValues_), + make_view(gpuBuffers.yValues_)); +} + +} // namespace Opm::gpuistl +#endif // HAVE_CUDA + #endif diff --git a/opm/material/fluidmatrixinteractions/SatCurveMultiplexerParams.hpp b/opm/material/fluidmatrixinteractions/SatCurveMultiplexerParams.hpp index c48e6da24b3..001109a6881 100644 --- a/opm/material/fluidmatrixinteractions/SatCurveMultiplexerParams.hpp +++ b/opm/material/fluidmatrixinteractions/SatCurveMultiplexerParams.hpp @@ -120,7 +120,6 @@ 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/fluidsystems/BlackOilFluidSystem_macrotemplate.hpp b/opm/material/fluidsystems/BlackOilFluidSystem_macrotemplate.hpp index 452bb970461..3ac5b5e6992 100644 --- a/opm/material/fluidsystems/BlackOilFluidSystem_macrotemplate.hpp +++ b/opm/material/fluidsystems/BlackOilFluidSystem_macrotemplate.hpp @@ -180,7 +180,6 @@ class FLUIDSYSTEM_CLASSNAME : public BaseFluidSystem(T); + // return params.template eval(T); + return params.internalEnergyFunction().eval(T, /*extrapolate=*/true); } }; diff --git a/opm/material/thermal/EclSpecrockLawParams.hpp b/opm/material/thermal/EclSpecrockLawParams.hpp index 6e2f1418ff3..afed81aa8ae 100644 --- a/opm/material/thermal/EclSpecrockLawParams.hpp +++ b/opm/material/thermal/EclSpecrockLawParams.hpp @@ -33,12 +33,14 @@ #include #include +#include #include #include #include #include #include +#include namespace Opm { @@ -68,24 +70,25 @@ namespace Opm { * 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. + * \c Tabulated1DFunction. The function's storage policy is propagated so + * this class can be instantiated as a CPU object, an owning GPU object, + * or a non-owning GPU view usable from a kernel. */ template class Storage> class EclSpecrockLawParams : public EnsureFinalized { + using InternalEnergyFunction = + Tabulated1DFunction; + public: using Scalar = ScalarT; - using ValueVector = Storage; + using ValueVector = typename InternalEnergyFunction::ValueVector; OPM_HOST_DEVICE EclSpecrockLawParams() = default; - OPM_HOST_DEVICE EclSpecrockLawParams(ValueVector temperatureSamples, - ValueVector internalEnergySamples) - : temperatureSamples_(std::move(temperatureSamples)) - , internalEnergySamples_(std::move(internalEnergySamples)) + OPM_HOST_DEVICE explicit EclSpecrockLawParams( + InternalEnergyFunction internalEnergyFunction) + : internalEnergyFunction_(std::move(internalEnergyFunction)) { EnsureFinalized::finalize(); } @@ -111,34 +114,32 @@ class EclSpecrockLawParams : public EnsureFinalized "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]); + const std::size_t n = temperature.size(); + std::vector temperatures(n); + std::vector internalEnergies(n); + Scalar curU = temperature[0] * heatCapacity[0]; for (std::size_t i = 0; i < n; ++i) { - temperatureSamples_[i] = static_cast(temperature[i]); - internalEnergySamples_[i] = curU; + temperatures[i] = temperature[i]; + internalEnergies[i] = curU; - if (i + 1 >= n) { + if (i >= n - 1) { break; } - // 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); + // Integrate the heat capacity from the current sampling point + // to the next one using the trapezoidal rule. + curU += 0.5 * (heatCapacity[i] + heatCapacity[i + 1]) + * (temperature[i + 1] - temperature[i]); } + internalEnergyFunction_.setXYContainers(temperatures, internalEnergies); } /*! - * \brief Set the sample tables directly. Marks the object as finalized. - * - * Available only on the CPU instantiation. + * \brief Set the sample tables directly. Marks the object as finalized. + * + * This compatibility helper is kept for callers that construct + * SPECROCK parameters directly instead of from heat capacities. */ template , @@ -151,84 +152,35 @@ class EclSpecrockLawParams : public EnsureFinalized "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]); - } + internalEnergyFunction_.setXYContainers(temperature, internalEnergy); EnsureFinalized::finalize(); } - OPM_HOST_DEVICE std::size_t numSamples() const - { - return temperatureSamples_.size(); - } + OPM_HOST_DEVICE const InternalEnergyFunction& internalEnergyFunction() const + { return internalEnergyFunction_; } - /*! - * \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 std::size_t numSamples() const + { return internalEnergyFunction_.numSamples(); } - OPM_HOST_DEVICE const ValueVector& temperatureSamples() const - { - EnsureFinalized::check(); - return temperatureSamples_; - } + const ValueVector& temperatureSamples() const + { return internalEnergyFunction_.xValues(); } - OPM_HOST_DEVICE const ValueVector& internalEnergySamples() const - { - EnsureFinalized::check(); - return internalEnergySamples_; - } - - ValueVector& temperatureSamplesMutable() - { - return temperatureSamples_; - } - - ValueVector& internalEnergySamplesMutable() - { - return internalEnergySamples_; - } + const ValueVector& internalEnergySamples() const + { return internalEnergyFunction_.yValues(); } private: - ValueVector temperatureSamples_ {}; - ValueVector internalEnergySamples_ {}; +#if HAVE_CUDA + template + friend ::Opm::EclSpecrockLawParams + Opm::gpuistl::copy_to_gpu(const ::Opm::EclSpecrockLawParams& cpu); + + template class ContainerT> + friend ::Opm::EclSpecrockLawParams + Opm::gpuistl::make_view( + ::Opm::EclSpecrockLawParams& gpuBuffers); +#endif + + InternalEnergyFunction internalEnergyFunction_ {}; }; } // namespace Opm @@ -241,17 +193,15 @@ ::Opm::EclSpecrockLawParams copy_to_gpu(const ::Opm::EclSpecrockLawParams& cpu) { return ::Opm::EclSpecrockLawParams( - GpuBuffer(cpu.temperatureSamples()), - GpuBuffer(cpu.internalEnergySamples())); + copy_to_gpu(cpu.internalEnergyFunction_)); } 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); + return ::Opm::EclSpecrockLawParams( + make_view(gpuBuffers.internalEnergyFunction_)); } } // namespace Opm::gpuistl diff --git a/opm/material/thermal/GpuEclThermalLawManager.hpp b/opm/material/thermal/GpuEclThermalLawManager.hpp index 97c15d61b40..eef76eb3933 100644 --- a/opm/material/thermal/GpuEclThermalLawManager.hpp +++ b/opm/material/thermal/GpuEclThermalLawManager.hpp @@ -94,7 +94,8 @@ namespace detail { #if HAVE_CUDA template struct SpecrockSampleHolder { - std::vector<::Opm::gpuistl::GpuBuffer> sampleBuffers {}; + std::vector<::Opm::EclSpecrockLawParams< + Scalar, ::Opm::gpuistl::GpuBuffer>> sampleBuffers {}; }; #endif @@ -152,7 +153,8 @@ class GpuManager */ #if HAVE_CUDA template = 0> - GpuManager(std::vector<::Opm::gpuistl::GpuBuffer> sampleBuffers, + GpuManager(std::vector<::Opm::EclSpecrockLawParams< + Scalar, ::Opm::gpuistl::GpuBuffer>> sampleBuffers, OuterStorage solidEnergyParams, OuterStorage elementToSolidRegionIdx, OuterStorage thermalConductionParams) @@ -213,7 +215,8 @@ class GpuManager */ #if HAVE_CUDA template = 0> - std::vector<::Opm::gpuistl::GpuBuffer>& specrockSampleBuffers() + std::vector<::Opm::EclSpecrockLawParams< + Scalar, ::Opm::gpuistl::GpuBuffer>>& specrockSampleBuffers() { return this->sampleBuffers; } @@ -328,21 +331,17 @@ copy_to_gpu(const ::Opm::EclThermalLaw::GpuManager& cpu) 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()); + using SolidEnergyLawParamsBuffer = + ::Opm::EclSpecrockLawParams; + std::vector sampleBuffers; + sampleBuffers.reserve(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)); + sampleBuffers.emplace_back(::Opm::gpuistl::copy_to_gpu(cpuRegion)); + hostStaging.emplace_back(::Opm::gpuistl::make_view(sampleBuffers.back())); } return ManagerBuf(std::move(sampleBuffers), From 95df78d3bc80a05e864ae6ce04e6b6f7f5f4c344 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Thu, 6 Aug 2026 13:40:30 +0200 Subject: [PATCH 04/10] reduce diff --- opm/material/thermal/EclSpecrockLaw.hpp | 1 - opm/material/thermal/EclSpecrockLawParams.hpp | 54 +++++++++---------- 2 files changed, 26 insertions(+), 29 deletions(-) diff --git a/opm/material/thermal/EclSpecrockLaw.hpp b/opm/material/thermal/EclSpecrockLaw.hpp index 6595b6e307d..13799c1f1b0 100644 --- a/opm/material/thermal/EclSpecrockLaw.hpp +++ b/opm/material/thermal/EclSpecrockLaw.hpp @@ -56,7 +56,6 @@ class EclSpecrockLaw OPM_HOST_DEVICE static Evaluation solidInternalEnergy(const Params& params, const FluidState& fluidState) { const auto& T = fluidState.temperature(/*phaseIdx=*/0); - // return params.template eval(T); return params.internalEnergyFunction().eval(T, /*extrapolate=*/true); } }; diff --git a/opm/material/thermal/EclSpecrockLawParams.hpp b/opm/material/thermal/EclSpecrockLawParams.hpp index afed81aa8ae..64037004c6b 100644 --- a/opm/material/thermal/EclSpecrockLawParams.hpp +++ b/opm/material/thermal/EclSpecrockLawParams.hpp @@ -86,6 +86,8 @@ class EclSpecrockLawParams : public EnsureFinalized OPM_HOST_DEVICE EclSpecrockLawParams() = default; + OPM_HOST_DEVICE EclSpecrockLawParams(const EclSpecrockLawParams&) = default; + OPM_HOST_DEVICE explicit EclSpecrockLawParams( InternalEnergyFunction internalEnergyFunction) : internalEnergyFunction_(std::move(internalEnergyFunction)) @@ -108,31 +110,30 @@ class EclSpecrockLawParams : public EnsureFinalized void setHeatCapacities(const ContainerT& temperature, const ContainerT& heatCapacity) { - if (temperature.size() != heatCapacity.size()) { - OPM_THROW(std::invalid_argument, - "EclSpecrockLawParams: temperature and heat-capacity arrays must have " - "matching sizes"); - } - - // Integrate the heat capacity to compute the internal energy. - const std::size_t n = temperature.size(); - std::vector temperatures(n); - std::vector internalEnergies(n); - Scalar curU = temperature[0] * heatCapacity[0]; - for (std::size_t i = 0; i < n; ++i) { - temperatures[i] = temperature[i]; - internalEnergies[i] = curU; - - if (i >= n - 1) { + 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) break; - } - // Integrate the heat capacity from the current sampling point - // to the next one using the trapezoidal rule. - curU += 0.5 * (heatCapacity[i] + heatCapacity[i + 1]) - * (temperature[i + 1] - temperature[i]); + // 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); } - internalEnergyFunction_.setXYContainers(temperatures, internalEnergies); + + internalEnergyFunction_.setXYContainers(T, u); } /*! @@ -157,16 +158,13 @@ class EclSpecrockLawParams : public EnsureFinalized } OPM_HOST_DEVICE const InternalEnergyFunction& internalEnergyFunction() const - { return internalEnergyFunction_; } - - OPM_HOST_DEVICE std::size_t numSamples() const - { return internalEnergyFunction_.numSamples(); } + { EnsureFinalized::check(); return internalEnergyFunction_; } const ValueVector& temperatureSamples() const - { return internalEnergyFunction_.xValues(); } + { EnsureFinalized::check(); return internalEnergyFunction_.xValues(); } const ValueVector& internalEnergySamples() const - { return internalEnergyFunction_.yValues(); } + { EnsureFinalized::check(); return internalEnergyFunction_.yValues(); } private: #if HAVE_CUDA From 0788dc31a594cbba05e34bcca4f86c0961c2038a Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Thu, 6 Aug 2026 13:55:09 +0200 Subject: [PATCH 05/10] reduce diff --- opm/material/thermal/EclThconrLaw.hpp | 25 +++++++------------------ 1 file changed, 7 insertions(+), 18 deletions(-) diff --git a/opm/material/thermal/EclThconrLaw.hpp b/opm/material/thermal/EclThconrLaw.hpp index e3802dce54d..99399dd7cf9 100644 --- a/opm/material/thermal/EclThconrLaw.hpp +++ b/opm/material/thermal/EclThconrLaw.hpp @@ -57,26 +57,15 @@ class EclThconrLaw const FluidState& fluidState) { // THCONR + THCONSF approach. - 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 { - gasActive = fluidState.fluidSystem().phaseIsActive(gasPhaseIdx); - } - - if (gasActive) { - const Scalar alpha = params.dTotalThermalConductivity_dSg(); + Scalar lambdaRef = params.referenceTotalThermalConductivity(); + static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx; + if (fluidState.fluidSystem().phaseIsActive(gasPhaseIdx)) { + Scalar alpha = params.dTotalThermalConductivity_dSg(); const Evaluation& Sg = decay(fluidState.saturation(gasPhaseIdx)); - return lambdaRef * (Scalar(1) - alpha * Sg); + return lambdaRef*(1.0 - alpha*Sg); + } else { + return lambdaRef; } - return Evaluation(lambdaRef); } }; From 4355cafaefed8fb05050e7c6d40679c5fb8af36d Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Thu, 6 Aug 2026 13:58:09 +0200 Subject: [PATCH 06/10] reduce diff --- opm/material/thermal/EclThermalLawManager.hpp | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/opm/material/thermal/EclThermalLawManager.hpp b/opm/material/thermal/EclThermalLawManager.hpp index 9a3e7364923..9f1c8e7ffb4 100644 --- a/opm/material/thermal/EclThermalLawManager.hpp +++ b/opm/material/thermal/EclThermalLawManager.hpp @@ -78,22 +78,11 @@ class EclThermalLawManager 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. - */ + /*! \brief Return the element-index to SATNUM-region mapping */ 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. - */ + /*! \brief Return the per-region solid-energy law parameter vector. */ const std::vector& solidEnergyLawParamsVector() const { return solidEnergyLawParams_; } From 2f611a1f723a660ea70aa88af7ea142d2bcc62c7 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Fri, 7 Aug 2026 11:35:09 +0200 Subject: [PATCH 07/10] unify tabulation evaluation for cpu and gpu --- opm/material/common/Tabulated1DFunction.hpp | 82 ++++--------------- opm/material/thermal/EclSpecrockLaw.hpp | 1 + opm/material/thermal/EclSpecrockLawParams.hpp | 6 ++ .../thermal/GpuEclThermalLawManager.hpp | 8 ++ 4 files changed, 32 insertions(+), 65 deletions(-) diff --git a/opm/material/common/Tabulated1DFunction.hpp b/opm/material/common/Tabulated1DFunction.hpp index 1fea415a8f3..c5011765583 100644 --- a/opm/material/common/Tabulated1DFunction.hpp +++ b/opm/material/common/Tabulated1DFunction.hpp @@ -55,10 +55,6 @@ class Tabulated1DFunction; #if HAVE_CUDA namespace gpuistl { -template -Tabulated1DFunction -copy_to_gpu(const Tabulated1DFunction& cpu); - template class ContainerT> Tabulated1DFunction make_view(Tabulated1DFunction& gpuBuffers); @@ -243,7 +239,7 @@ class Tabulated1DFunction /*! * \brief Returns the number of sampling points. */ - std::size_t numSamples() const + OPM_HOST_DEVICE std::size_t numSamples() const { return xValues_.size(); } /*! @@ -280,7 +276,7 @@ class Tabulated1DFunction * \brief Return true iff the given x is in range [x1, xn]. */ template - bool applies(const Evaluation& x) const + OPM_HOST_DEVICE bool applies(const Evaluation& x) const { return xValues_[0] <= x && x <= xValues_[numSamples() - 1]; } /*! @@ -295,13 +291,8 @@ class Tabulated1DFunction template OPM_HOST_DEVICE Evaluation eval(const Evaluation& x, bool extrapolate = false) const { - if constexpr (OPM_IS_INSIDE_DEVICE_FUNCTION) { - return evalDevice_(x); - } - else { - SegmentIndex segIdx = findSegmentIndex(x, extrapolate); - return eval(x, segIdx); - } + SegmentIndex segIdx = findSegmentIndex(x, extrapolate); + return eval(x, segIdx); } template @@ -313,7 +304,6 @@ class Tabulated1DFunction Scalar y0 = yValues_[segIdx]; Scalar y1 = yValues_[segIdx + 1]; - return y0 + (y1 - y0)*(x - x0)/(x1 - x0); } @@ -471,25 +461,18 @@ class Tabulated1DFunction } template - SegmentIndex findSegmentIndex(const Evaluation& x, bool extrapolate = false) const + OPM_HOST_DEVICE SegmentIndex findSegmentIndex(const Evaluation& x, bool extrapolate = false) const { - if (!isfinite(x)) { - throw std::runtime_error("We can not search for extrapolation/interpolation " - "segment in an 1D table for non-finite value " + - std::to_string(getValue(x)) + " ."); - } - if (!extrapolate && !applies(x)) - throw std::logic_error("Trying to evaluate a tabulated function outside of its range"); +#if OPM_IS_INSIDE_HOST_FUNCTION + // relies on std::isfinite() which is not supported in device code + OPM_ERROR_IF(!isfinite(x), "Trying to evaluate a tabulated function at a non-finite value"); +#endif - // we need at least two sampling points! - if (numSamples() < 2) { - throw std::logic_error("We need at least two sampling points to " - "do interpolation/extrapolation, " - "and the table only contains " + - std::to_string(numSamples()) + - " sampling points"); - } + OPM_ERROR_IF(numSamples() < 2, "Trying to evaluate a tabulated function with less than two samples"); + + OPM_ERROR_IF(!extrapolate && !applies(x), + "Trying to evaluate a tabulated function outside of its range"); if (x <= xValues_[1]) return SegmentIndex{0}; @@ -508,6 +491,7 @@ class Tabulated1DFunction } if (xValues_[lowerIdx] > x || x > xValues_[lowerIdx + 1]) { +#if OPM_IS_INSIDE_HOST_FUNCTION std::string msg = "Problematic interpolation/extrapolation " "segment is found for the input value " + std::to_string(Opm::getValue(x)) + @@ -534,43 +518,15 @@ class Tabulated1DFunction msg += "\n"; OpmLog::debug(msg); throw std::runtime_error(msg); +#else + OPM_THROW(std::runtime_error, "Problematic interpolation/extrapolation segment found"); +#endif } return SegmentIndex{lowerIdx}; } } private: - template - OPM_HOST_DEVICE Evaluation evalDevice_(const Evaluation& x) const - { - const std::size_t n = xValues_.size(); - assert(n >= 2); - - std::size_t segIdx = 0; - if (x <= xValues_[1]) { - segIdx = 0; - } - else if (x >= xValues_[n - 2]) { - segIdx = n - 2; - } - else { - std::size_t lowerIdx = 1; - std::size_t upperIdx = n - 2; - while (lowerIdx + 1 < upperIdx) { - const std::size_t pivotIdx = (lowerIdx + upperIdx) / 2; - if (x < xValues_[pivotIdx]) { - upperIdx = pivotIdx; - } - else { - lowerIdx = pivotIdx; - } - } - segIdx = lowerIdx; - } - - return eval(x, SegmentIndex{segIdx}); - } - template Evaluation evalDerivative_(const Evaluation& x, std::size_t segIdx) const { @@ -682,10 +638,6 @@ class Tabulated1DFunction } #if HAVE_CUDA - template - friend Tabulated1DFunction - gpuistl::copy_to_gpu(const Tabulated1DFunction& cpu); - template class ContainerT> friend Tabulated1DFunction gpuistl::make_view(Tabulated1DFunction& gpuBuffers); diff --git a/opm/material/thermal/EclSpecrockLaw.hpp b/opm/material/thermal/EclSpecrockLaw.hpp index 13799c1f1b0..3030abe2363 100644 --- a/opm/material/thermal/EclSpecrockLaw.hpp +++ b/opm/material/thermal/EclSpecrockLaw.hpp @@ -30,6 +30,7 @@ #include "EclSpecrockLawParams.hpp" #include +#include namespace Opm { diff --git a/opm/material/thermal/EclSpecrockLawParams.hpp b/opm/material/thermal/EclSpecrockLawParams.hpp index 64037004c6b..3a060339522 100644 --- a/opm/material/thermal/EclSpecrockLawParams.hpp +++ b/opm/material/thermal/EclSpecrockLawParams.hpp @@ -134,6 +134,7 @@ class EclSpecrockLawParams : public EnsureFinalized } internalEnergyFunction_.setXYContainers(T, u); + EnsureFinalized::finalize(); } /*! @@ -190,6 +191,11 @@ template ::Opm::EclSpecrockLawParams copy_to_gpu(const ::Opm::EclSpecrockLawParams& cpu) { + const auto& cpuFunction = cpu.internalEnergyFunction(); + if (cpuFunction.numSamples() < 2) { + OPM_THROW(std::logic_error, + "SPECROCK GPU upload requires at least two table samples"); + } return ::Opm::EclSpecrockLawParams( copy_to_gpu(cpu.internalEnergyFunction_)); } diff --git a/opm/material/thermal/GpuEclThermalLawManager.hpp b/opm/material/thermal/GpuEclThermalLawManager.hpp index eef76eb3933..46752af1737 100644 --- a/opm/material/thermal/GpuEclThermalLawManager.hpp +++ b/opm/material/thermal/GpuEclThermalLawManager.hpp @@ -170,6 +170,14 @@ class GpuManager OPM_HOST_DEVICE SolidEnergyLawParams solidEnergyLawParams(unsigned elemIdx) const { const int regionIdx = elementToSolidRegionIdx_[elemIdx]; + + OPM_ERROR_IF(elemIdx >= elementToSolidRegionIdx_.size() + || regionIdx < 0, + "Invalid element index or region index in GpuManager::solidEnergyLawParams"); + OPM_ERROR_IF(solidEnergyParams_.size() < 0 + || static_cast(regionIdx) >= solidEnergyParams_.size(), + "Invalid region index in GpuManager::solidEnergyLawParams"); + return solidEnergyParams_[regionIdx]; } From b97723b73dd7eb8c59f689772fae1131ebdaa884 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Fri, 7 Aug 2026 12:55:53 +0200 Subject: [PATCH 08/10] refine GpuEclMatieralLawManager.hpp --- .../GpuEclMaterialLawManager.hpp | 46 ++++--------------- 1 file changed, 8 insertions(+), 38 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp index 4469c7a0527..f3fe3baf225 100644 --- a/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp +++ b/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp @@ -59,15 +59,6 @@ #include #include -namespace Opm::gpuistl -{ -template -class GpuBuffer; - -template -class GpuView; -} // namespace Opm::gpuistl - namespace Opm::EclMaterialLaw { @@ -214,29 +205,9 @@ class GpuManager : private detail::GpuPiecewiseLinearSampleHolder< { } - /*! - * \brief Construct from a CPU \c Opm::EclMaterialLaw::Manager directly - * into device-resident \c GpuBuffer storage. - * - * For each SATNUM region, the 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. - */ /*! \brief Material-law parameters of an active cell. */ - OPM_HOST_DEVICE MaterialLawParams materialLawParams(unsigned elemIdx) const + OPM_HOST_DEVICE const 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 materialLawParamsByRegion_[satnumRegionArray_[elemIdx]]; } @@ -346,16 +317,15 @@ class GpuManager : private detail::GpuPiecewiseLinearSampleHolder< GpuData data; data.satnumRegionArray = buildHostSatnumRegionArray(cpu, numElements); - const auto maxRegion = data.satnumRegionArray.empty() - ? -1 - : *std::max_element(data.satnumRegionArray.begin(), data.satnumRegionArray.end()); - if (maxRegion < 0) { - return data; - } - data.materialLawParamsByRegion.resize(static_cast(maxRegion) + 1u); + auto maxRegion = std::max_element(data.satnumRegionArray.begin(), data.satnumRegionArray.end()); + OPM_ERROR_IF(maxRegion == data.satnumRegionArray.end(), + "Failed to find maximum SATNUM region index in GPU material manager"); + + data.materialLawParamsByRegion.resize(static_cast(*maxRegion) + 1u); std::vector initialized(data.materialLawParamsByRegion.size(), false); std::vector firstElement(data.materialLawParamsByRegion.size(), 0u); - data.sampleBuffers.reserve(data.materialLawParamsByRegion.size() * 12u); + // reserve 6 or 12 sample buffers depending 2 or three-phase + data.sampleBuffers.reserve(data.materialLawParamsByRegion.size() * (6u + 6u * static_cast(isTwoPhase))); auto pushSampleBuffer = [&](const auto& sampleVector) { std::vector hostCopy(sampleVector.begin(), sampleVector.end()); From 7e0b79bb7a405a32b308d596bdbf95c641cff991 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Fri, 7 Aug 2026 13:12:56 +0200 Subject: [PATCH 09/10] set correct attribution --- opm/material/thermal/GpuEclThermalLawManager.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/material/thermal/GpuEclThermalLawManager.hpp b/opm/material/thermal/GpuEclThermalLawManager.hpp index 46752af1737..bdfd705f394 100644 --- a/opm/material/thermal/GpuEclThermalLawManager.hpp +++ b/opm/material/thermal/GpuEclThermalLawManager.hpp @@ -1,7 +1,7 @@ // -*- 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 + Copyright 2026 Equinor ASA This file is part of the Open Porous Media project (OPM). From 512139fdc5c457200ecf5463b1a6f88bd5725248 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Fri, 18 Sep 2026 10:54:53 +0200 Subject: [PATCH 10/10] resolving PR feedback --- .../fluidmatrixinteractions/GpuEclMaterialLawManager.hpp | 3 +++ opm/material/thermal/EclSpecrockLawParams.hpp | 2 ++ opm/material/thermal/GpuEclThermalLawManager.hpp | 9 +++++---- 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp index f3fe3baf225..e9583886013 100644 --- a/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp +++ b/opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp @@ -85,6 +85,9 @@ namespace detail params.template getRealParams< ::Opm::SatCurveMultiplexerApproach::PiecewiseLinear>(); }) { + OPM_ERROR_IF(params.approach() + != ::Opm::SatCurveMultiplexerApproach::PiecewiseLinear, + "GPU material manager only supports piecewise-linear saturation curves"); return params .template getRealParams<::Opm::SatCurveMultiplexerApproach::PiecewiseLinear>(); } else { diff --git a/opm/material/thermal/EclSpecrockLawParams.hpp b/opm/material/thermal/EclSpecrockLawParams.hpp index 3a060339522..06428fd56ba 100644 --- a/opm/material/thermal/EclSpecrockLawParams.hpp +++ b/opm/material/thermal/EclSpecrockLawParams.hpp @@ -88,6 +88,8 @@ class EclSpecrockLawParams : public EnsureFinalized OPM_HOST_DEVICE EclSpecrockLawParams(const EclSpecrockLawParams&) = default; + OPM_HOST_DEVICE EclSpecrockLawParams(EclSpecrockLawParams&&) = default; + OPM_HOST_DEVICE explicit EclSpecrockLawParams( InternalEnergyFunction internalEnergyFunction) : internalEnergyFunction_(std::move(internalEnergyFunction)) diff --git a/opm/material/thermal/GpuEclThermalLawManager.hpp b/opm/material/thermal/GpuEclThermalLawManager.hpp index bdfd705f394..d90f2624123 100644 --- a/opm/material/thermal/GpuEclThermalLawManager.hpp +++ b/opm/material/thermal/GpuEclThermalLawManager.hpp @@ -169,12 +169,13 @@ class GpuManager /*! \brief Solid-energy law parameters for an active cell. */ OPM_HOST_DEVICE SolidEnergyLawParams solidEnergyLawParams(unsigned elemIdx) const { - const int regionIdx = elementToSolidRegionIdx_[elemIdx]; - OPM_ERROR_IF(elemIdx >= elementToSolidRegionIdx_.size() - || regionIdx < 0, + OPM_ERROR_IF(elemIdx >= elementToSolidRegionIdx_.size(), "Invalid element index or region index in GpuManager::solidEnergyLawParams"); - OPM_ERROR_IF(solidEnergyParams_.size() < 0 + + const int regionIdx = elementToSolidRegionIdx_[elemIdx]; + + OPM_ERROR_IF(regionIdx < 0 || static_cast(regionIdx) >= solidEnergyParams_.size(), "Invalid region index in GpuManager::solidEnergyLawParams");