From 795cf27e951ad6a040d68efe6e1af21ca9944743 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Thu, 30 Apr 2026 21:48:03 +0000 Subject: [PATCH] GpuFlowProblem and GpuEcl{Material,Thermal}LawManager Adds the GPU-side problem and material/thermal law manager headers (GpuFlowProblem, GpuEclMaterialLawManager, GpuEclThermalLawManager) plus the matching GPU element-context header (fvbaseelementcontextgpu.hh) and the FlowGasWaterEnergyTypeTag.hpp typetag. Wires them through FlowProblemParameters and adds a Boost test (test_gpu_ecl_thermal_law_manager.cu) backed by tests/very_simple_deck.DATA. The new test is gated on HIP or CUDA >= 13.1 in CMakeLists_files.cmake. --- CMakeLists.txt | 2 + CMakeLists_files.cmake | 9 + .../common/fvbaseelementcontextgpu.hh | 309 +++++++++ .../flow/FlowGasWaterEnergyTypeTag.hpp | 99 +++ opm/simulators/flow/FlowProblemParameters.cpp | 18 + opm/simulators/flow/FlowProblemParameters.hpp | 9 + .../flow/GpuEclMaterialLawManager.hpp | 461 ++++++++++++++ .../flow/GpuEclThermalLawManager.hpp | 367 +++++++++++ opm/simulators/flow/GpuFlowProblem.hpp | 587 ++++++++++++++++++ .../test_gpu_ecl_thermal_law_manager.cu | 210 +++++++ tests/very_simple_deck.DATA | 53 ++ 11 files changed, 2124 insertions(+) create mode 100644 opm/models/discretization/common/fvbaseelementcontextgpu.hh create mode 100644 opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp create mode 100644 opm/simulators/flow/GpuEclMaterialLawManager.hpp create mode 100644 opm/simulators/flow/GpuEclThermalLawManager.hpp create mode 100644 opm/simulators/flow/GpuFlowProblem.hpp create mode 100644 tests/gpuistl/test_gpu_ecl_thermal_law_manager.cu create mode 100644 tests/very_simple_deck.DATA diff --git a/CMakeLists.txt b/CMakeLists.txt index 18f9395e5f9..c4988db2f0c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -366,6 +366,7 @@ macro(opm-simulators_sources_hook) set_source_files_properties( tests/gpuistl/test_blackoilfluidstategpu.cu tests/gpuistl/test_gpu_ad.cu + tests/gpuistl/test_gpu_ecl_thermal_law_manager.cu tests/gpuistl/test_gpu_linear_two_phase_material.cu tests/gpuistl/test_gpuPvt.cu tests/gpuistl/test_gpuBlackOilFluidSystem.cu @@ -490,6 +491,7 @@ macro(opm-simulators_tests_hook) cuVector_operations deviceBlockOperations gpu_ad + gpu_ecl_thermal_law_manager gpu_linear_two_phase_material gpu_resources gpu_safe_call diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 6075ba58577..b01a139711a 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -340,6 +340,10 @@ if(CUDA_FOUND OR hip_FOUND) # HEADERS + list(APPEND PUBLIC_HEADER_FILES opm/models/discretization/common/fvbaseelementcontextgpu.hh) + list(APPEND PUBLIC_HEADER_FILES opm/simulators/flow/GpuEclMaterialLawManager.hpp) + list(APPEND PUBLIC_HEADER_FILES opm/simulators/flow/GpuEclThermalLawManager.hpp) + list(APPEND PUBLIC_HEADER_FILES opm/simulators/flow/GpuFlowProblem.hpp) ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/autotuner.hpp) ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/coloringAndReorderingUtils.hpp) ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/gpu_safe_call.hpp) @@ -570,6 +574,9 @@ if(CUDA_FOUND OR hip_FOUND) ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_safe_conversion.cpp) ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_solver_adapter.cpp) ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_gpu_ad.cu) + if(hip_FOUND OR CUDA_VERSION VERSION_GREATER_EQUAL 13.1) + ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_gpu_ecl_thermal_law_manager.cu) + endif() ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_gpu_linear_two_phase_material.cu) ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_gpuPvt.cu) ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_gpu_smart_pointers.cu) @@ -707,6 +714,7 @@ list (APPEND TEST_DATA_FILES tests/data/test_stokes2c.dgf tests/data/test_stokes2cni.dgf tests/data/waterair.dgf + tests/very_simple_deck.DATA ) @@ -969,6 +977,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/FIPContainer.hpp opm/simulators/flow/FlowBaseProblemProperties.hpp opm/simulators/flow/FlowBaseVanguard.hpp + opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp opm/simulators/flow/FlowGenericProblem.hpp opm/simulators/flow/FlowGenericProblem_impl.hpp opm/simulators/flow/FlowGenericVanguard.hpp diff --git a/opm/models/discretization/common/fvbaseelementcontextgpu.hh b/opm/models/discretization/common/fvbaseelementcontextgpu.hh new file mode 100644 index 00000000000..f7237ea9eb0 --- /dev/null +++ b/opm/models/discretization/common/fvbaseelementcontextgpu.hh @@ -0,0 +1,309 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief GPU-compatible stub of FvBaseElementContext. + * + * All methods are annotated with OPM_HOST_DEVICE and return dummy/default + * values. This class is intended as a placeholder so that GPU kernels that + * receive an element context can compile; it does not implement any real + * behaviour. + * + * The heavy Dune/grid infrastructure present in FvBaseElementContext is + * intentionally absent here because those types are not usable inside a + * CUDA/HIP kernel. + */ +#ifndef OPM_FV_BASE_ELEMENT_CONTEXT_GPU_HH +#define OPM_FV_BASE_ELEMENT_CONTEXT_GPU_HH + +#include +#include +#include + +#include + +namespace Opm +{ + +/*! + * \ingroup FiniteVolumeDiscretizations + * + * \brief GPU-compatible stub mirroring the public interface of FvBaseElementContext. + * + * Uses the same TypeTag-based template parameter as FvBaseElementContext; all + * concrete types are extracted via GetPropType at compile time. + * + * \tparam TypeTag The type tag from which Scalar, PrimaryVariables, etc. are derived. + */ +template +class FvBaseElementContextGpu +{ +public: + using Scalar = GetPropType; + using PrimaryVariables = GetPropType; + using IntensiveQuantities = GetPropType; + using ExtensiveQuantities = GetPropType; + using Problem = GetPropType; + using Model = GetPropType; + using GradientCalculator = GetPropType; + + // ----------------------------------------------------------------------- + // Construction / assignment + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE FvBaseElementContextGpu() = default; + OPM_HOST_DEVICE ~FvBaseElementContextGpu() = default; + OPM_HOST_DEVICE FvBaseElementContextGpu(const FvBaseElementContextGpu&) = default; + OPM_HOST_DEVICE FvBaseElementContextGpu(FvBaseElementContextGpu&&) noexcept = default; + OPM_HOST_DEVICE FvBaseElementContextGpu& operator=(const FvBaseElementContextGpu&) = default; + OPM_HOST_DEVICE FvBaseElementContextGpu& operator=(FvBaseElementContextGpu&&) noexcept + = default; + + // ----------------------------------------------------------------------- + // Stencil / element update stubs + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE void updateAll() + { + } + OPM_HOST_DEVICE void updateStencil() + { + } + OPM_HOST_DEVICE void updatePrimaryStencil() + { + } + OPM_HOST_DEVICE void updateStencilTopology() + { + } + OPM_HOST_DEVICE void updateAllIntensiveQuantities() + { + } + OPM_HOST_DEVICE void updateIntensiveQuantities(unsigned /*timeIdx*/) + { + } + OPM_HOST_DEVICE void updatePrimaryIntensiveQuantities(unsigned /*timeIdx*/) + { + } + OPM_HOST_DEVICE void updateIntensiveQuantities(const PrimaryVariables& /*priVars*/, + unsigned /*dofIdx*/, + unsigned /*timeIdx*/) + { + } + OPM_HOST_DEVICE void updateAllExtensiveQuantities() + { + } + OPM_HOST_DEVICE void updateExtensiveQuantities(unsigned /*timeIdx*/) + { + } + + // ----------------------------------------------------------------------- + // Focus DOF + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE void setFocusDofIndex(unsigned dofIdx) + { + focusDofIdx_ = static_cast(dofIdx); + } + OPM_HOST_DEVICE unsigned focusDofIndex() const + { + return focusDofIdx_; + } + + // ----------------------------------------------------------------------- + // Linearization type + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE LinearizationType linearizationType() const + { + return LinearizationType {}; + } + + // ----------------------------------------------------------------------- + // Problem / model accessors — return references to dummy stored objects + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE const Problem& problem() const + { + return problem_; + } + OPM_HOST_DEVICE const Model& model() const + { + return model_; + } + + // ----------------------------------------------------------------------- + // DOF / face counts + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE std::size_t numDof(unsigned /*timeIdx*/) const + { + return 0; + } + OPM_HOST_DEVICE std::size_t numPrimaryDof(unsigned /*timeIdx*/) const + { + return 0; + } + OPM_HOST_DEVICE std::size_t numInteriorFaces(unsigned /*timeIdx*/) const + { + return 0; + } + OPM_HOST_DEVICE std::size_t numBoundaryFaces(unsigned /*timeIdx*/) const + { + return 0; + } + + // ----------------------------------------------------------------------- + // Global space index / volume + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE unsigned globalSpaceIndex(unsigned /*dofIdx*/, unsigned /*timeIdx*/) const + { + return 0; + } + + OPM_HOST_DEVICE Scalar dofVolume(unsigned /*dofIdx*/, unsigned /*timeIdx*/) const + { + return Scalar {0}; + } + + OPM_HOST_DEVICE Scalar dofTotalVolume(unsigned /*dofIdx*/, unsigned /*timeIdx*/) const + { + return Scalar {0}; + } + + // ----------------------------------------------------------------------- + // Boundary + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE bool onBoundary() const + { + return false; + } + + // ----------------------------------------------------------------------- + // Intensive quantities + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE const IntensiveQuantities& intensiveQuantities(unsigned /*dofIdx*/, + unsigned /*timeIdx*/) const + { + return intensiveQuantitiesStashed_; + } + + OPM_HOST_DEVICE IntensiveQuantities& intensiveQuantities(unsigned /*dofIdx*/, + unsigned /*timeIdx*/) + { + return intensiveQuantitiesStashed_; + } + + OPM_HOST_DEVICE const IntensiveQuantities* thermodynamicHint(unsigned /*dofIdx*/, + unsigned /*timeIdx*/) const + { + return nullptr; + } + + // ----------------------------------------------------------------------- + // Primary variables + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE const PrimaryVariables& primaryVars(unsigned /*dofIdx*/, + unsigned /*timeIdx*/) const + { + return priVarsStashed_; + } + + // ----------------------------------------------------------------------- + // Stash / restore + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE bool haveStashedIntensiveQuantities() const + { + return stashedDofIdx_ != -1; + } + + OPM_HOST_DEVICE int stashedDofIdx() const + { + return stashedDofIdx_; + } + + OPM_HOST_DEVICE void stashIntensiveQuantities(unsigned dofIdx) + { + stashedDofIdx_ = static_cast(dofIdx); + } + + OPM_HOST_DEVICE void restoreIntensiveQuantities(unsigned /*dofIdx*/) + { + stashedDofIdx_ = -1; + } + + // ----------------------------------------------------------------------- + // Gradient calculator + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE const GradientCalculator& gradientCalculator() const + { + return gradientCalculator_; + } + + // ----------------------------------------------------------------------- + // Extensive quantities + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE const ExtensiveQuantities& extensiveQuantities(unsigned /*fluxIdx*/, + unsigned /*timeIdx*/) const + { + return extensiveQuantitiesStashed_; + } + + // ----------------------------------------------------------------------- + // Storage cache flag + // ----------------------------------------------------------------------- + + OPM_HOST_DEVICE bool enableStorageCache() const + { + return enableStorageCache_; + } + OPM_HOST_DEVICE void setEnableStorageCache(bool yesno) + { + enableStorageCache_ = yesno; + } + +private: + // Dummy stored objects returned by reference accessors. + // On device these are default-constructed and carry no meaningful data. + Problem problem_ {}; + Model model_ {}; + GradientCalculator gradientCalculator_ {}; + IntensiveQuantities intensiveQuantitiesStashed_ {}; + ExtensiveQuantities extensiveQuantitiesStashed_ {}; + PrimaryVariables priVarsStashed_ {}; + + int stashedDofIdx_ {-1}; + int focusDofIdx_ {-1}; + bool enableStorageCache_ {false}; +}; + +} // namespace Opm + +#endif // OPM_FV_BASE_ELEMENT_CONTEXT_GPU_HH diff --git a/opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp b/opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp new file mode 100644 index 00000000000..1c91e08b05f --- /dev/null +++ b/opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp @@ -0,0 +1,99 @@ +// -*- 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 + * + * Shared TypeTag declaration for the gas+water+energy (CO2STORE-style) + * blackoil flow problem. Originally inlined into + * \c flow/flow_gaswater_energy.cpp; extracted into a header so that other + * translation units (e.g. the GPU intensive-quantities dispatcher) can + * obtain matching property bindings via \c GetPropType. + */ +#ifndef OPM_FLOW_GASWATER_ENERGY_TYPETAG_HPP +#define OPM_FLOW_GASWATER_ENERGY_TYPETAG_HPP + +#include + +#include +#include +#include + +#include + +namespace Opm::Properties { +namespace TTag { +struct FlowGasWaterEnergyProblem { + using InheritsFrom = std::tuple; +}; +} // namespace TTag + +template +struct Linearizer { using type = TpfaLinearizer; }; + +template +struct LocalResidual { using type = BlackOilLocalResidualTPFA; }; + +template +struct EnableDiffusion { static constexpr bool value = true; }; + +template +struct EnableDispersion { static constexpr bool value = true; }; + +template +struct EnergyModuleType +{ static constexpr EnergyModules value = EnergyModules::FullyImplicitThermal; }; + +template +struct EnableDisgasInWater { + static constexpr bool value = true; +}; + +template +struct EnableVapwat { + static constexpr bool value = true; +}; + +//! The indices required by the model +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because + // this leads to cyclic definitions of some properties. + using BaseTypeTag = TTag::FlowProblem; + using FluidSystem = GetPropType; + static constexpr EnergyModules energyModuleType = getPropValue(); + static constexpr int numEnergyVars = (energyModuleType == EnergyModules::FullyImplicitThermal) ? 1 : 0; +public: + using type = BlackOilTwoPhaseIndices(), + getPropValue(), + getPropValue(), + numEnergyVars, + getPropValue(), + getPropValue(), + /*PVOffset=*/0, + /*disabledCompIdx=*/FluidSystem::oilCompIdx, + getPropValue()>; +}; + +} // namespace Opm::Properties + +#endif // OPM_FLOW_GASWATER_ENERGY_TYPETAG_HPP diff --git a/opm/simulators/flow/FlowProblemParameters.cpp b/opm/simulators/flow/FlowProblemParameters.cpp index 2f5f29668da..1daa0a09b2a 100644 --- a/opm/simulators/flow/FlowProblemParameters.cpp +++ b/opm/simulators/flow/FlowProblemParameters.cpp @@ -81,6 +81,24 @@ void registerFlowProblemParameters() ("Conserve inner energy and not enthalpy " "even if THERMAL is used."); +#if HAVE_CUDA && OPM_HAVE_GPU_BLACKOIL_INTENSIVE_QUANTITIES_DISPATCHER + Parameters::Register + ("Experimental: compute BlackOilIntensiveQuantities on the GPU " + "via the GpuBlackoilIntensiveQuantitiesDispatcher. Only takes " + "effect for CO2STORE-compatible TypeTags; ignored otherwise."); +#else + // Always register so that the run fails with an explanatory error + // message (rather than an "unknown parameter" message) when this binary + // has been built without the GPU intensive-quantities dispatcher (i.e. + // without CUDA/HIP support, or with a CUDA toolkit older than 13.1) + // but the user still requests the GPU dispatcher. + Parameters::Register + ("Experimental: compute BlackOilIntensiveQuantities on the GPU. " + "This binary was built without the GPU intensive-quantities " + "dispatcher (no CUDA/HIP support, or CUDA<13.1), so enabling " + "this option will cause the run to fail at start-up."); +#endif + // By default, stop it after the universe will probably have stopped // to exist. (the ECL problem will finish the simulation explicitly // after it simulated the last episode specified in the deck.) diff --git a/opm/simulators/flow/FlowProblemParameters.hpp b/opm/simulators/flow/FlowProblemParameters.hpp index 9908026b63e..0c2c307723d 100644 --- a/opm/simulators/flow/FlowProblemParameters.hpp +++ b/opm/simulators/flow/FlowProblemParameters.hpp @@ -69,6 +69,15 @@ struct UseHybridNewton { static constexpr bool value = false; }; // Conserve inner energy instead of enthalpy even if THERMAL is used struct ConserveInnerEnergyThermal { static constexpr bool value = false; }; +// Experimental: route the per-element BlackOilIntensiveQuantities update +// through the GPU dispatcher instead of computing it on the CPU. Only takes +// effect for CO2STORE-compatible TypeTags (gas+water, no oil). The parameter +// is always declared so that the run can fail with an explanatory error +// message when the option is requested in a binary that was built without +// CUDA/HIP support (or without a sufficiently new toolkit), or when the +// active TypeTag is incompatible with the GPU dispatcher. +struct ExperimentalComputePropertiesOnGpu { static constexpr bool value = false; }; + } // namespace Opm::Parameters namespace Opm { diff --git a/opm/simulators/flow/GpuEclMaterialLawManager.hpp b/opm/simulators/flow/GpuEclMaterialLawManager.hpp new file mode 100644 index 00000000000..46611f4c688 --- /dev/null +++ b/opm/simulators/flow/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/simulators/flow/GpuEclThermalLawManager.hpp b/opm/simulators/flow/GpuEclThermalLawManager.hpp new file mode 100644 index 00000000000..bc87a3c9884 --- /dev/null +++ b/opm/simulators/flow/GpuEclThermalLawManager.hpp @@ -0,0 +1,367 @@ +// -*- 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 + * + * 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. + * + * Iterates over the active cells, querying the CPU problem for each + * cell's \c SolidEnergyLawParams (multiplexer) and + * \c ThermalConductionLawParams (multiplexer). The per-cell SPECROCK + * multiplexer references are deduplicated by address to recover the + * underlying region index (so the resulting GPU manager stores one + * \c SolidEnergyLawParams per region rather than one per cell). + * + * Throws via \c OPM_THROW for any solid-energy approach other than + * SPECROCK or any thermal-conduction approach other than THCONR. + */ +template +GpuManager +buildCpuManagerFromFlowProblem(const CpuFlowProblemT& cpu, std::size_t numElements) +{ + using ManagerCpu = GpuManager; + using SolidEnergyLawParams = typename ManagerCpu::SolidEnergyLawParams; + using ThermalConductionLawParams = typename ManagerCpu::ThermalConductionLawParams; + + std::vector regionMultiplexerAddress; + std::vector hostElementToSolidRegionIdx(numElements); + std::vector hostSolidEnergyParams; + std::vector hostThermalConductionParams(numElements); + + for (std::size_t i = 0; i < numElements; ++i) { + const auto& cpuSolidMultiplexer + = cpu.solidEnergyLawParams(static_cast(i), 0u); + if (cpuSolidMultiplexer.solidEnergyApproach() + != ::Opm::EclSolidEnergyApproach::Specrock) { + OPM_THROW(std::logic_error, + "Opm::EclThermalLaw::GpuManager only supports the SPECROCK " + "solid-energy approach."); + } + const void* address = static_cast(&cpuSolidMultiplexer); + int regionIdx = -1; + for (std::size_t r = 0; r < regionMultiplexerAddress.size(); ++r) { + if (regionMultiplexerAddress[r] == address) { + regionIdx = static_cast(r); + break; + } + } + if (regionIdx < 0) { + const auto& cpuSpecrock = cpuSolidMultiplexer.template getRealParams< + ::Opm::EclSolidEnergyApproach::Specrock>(); + const auto& xValues = cpuSpecrock.temperatureSamples(); + const auto& yValues = cpuSpecrock.internalEnergySamples(); + std::vector temperatureSamples(xValues.begin(), xValues.end()); + std::vector internalEnergySamples(yValues.begin(), yValues.end()); + SolidEnergyLawParams params; + params.setSamples(temperatureSamples, internalEnergySamples); + regionIdx = static_cast(hostSolidEnergyParams.size()); + hostSolidEnergyParams.emplace_back(std::move(params)); + regionMultiplexerAddress.push_back(address); + } + hostElementToSolidRegionIdx[i] = regionIdx; + + const auto& cpuThermalMultiplexer + = cpu.thermalConductionLawParams(static_cast(i), 0u); + if (cpuThermalMultiplexer.thermalConductionApproach() + != ::Opm::EclThermalConductionApproach::Thconr) { + OPM_THROW(std::logic_error, + "Opm::EclThermalLaw::GpuManager only supports the THCONR " + "thermal-conduction approach."); + } + const auto& cpuThconr = cpuThermalMultiplexer.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; + + 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 diff --git a/opm/simulators/flow/GpuFlowProblem.hpp b/opm/simulators/flow/GpuFlowProblem.hpp new file mode 100644 index 00000000000..73d88ee2f3a --- /dev/null +++ b/opm/simulators/flow/GpuFlowProblem.hpp @@ -0,0 +1,587 @@ +// -*- 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::FlowProblem. + * + * Only exposes the subset of the FlowProblem / FlowProblemBlackoil interface + * that BlackOilIntensiveQuantities needs in order to call ::update(). Anything + * that is not used by the intensive quantities (well model, aquifers, source + * terms, boundary conditions, ...) is intentionally absent. + */ +#ifndef OPM_GPU_FLOW_PROBLEM_HPP +#define OPM_GPU_FLOW_PROBLEM_HPP + +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +namespace Opm::gpuistl +{ +template +class GpuBuffer; +template +class GpuView; +} // namespace Opm::gpuistl + +namespace Opm { + +/*! + * \brief No-op thermal-law manager used as the default 4th template arg of + * \c GpuFlowProblem. + * + * Carries no data and exposes the minimum surface needed by the + * \c GpuFlowProblem builder/copy/view machinery: a \c FluidSystem alias + * fixed to \c void (so the \c hasThermal compile-time switch evaluates + * to false), trivial empty params types, and trivial accessors that + * return default-constructed values. + * + * Intentionally not callable on the device when actually queried: it is + * only used by non-thermal \c GpuFlowProblem instantiations whose dispatch + * paths never call \c solidEnergyLawParams or \c thermalConductionLawParams. + */ +struct NoThermalLawManager +{ + using FluidSystem = void; + struct SolidEnergyLawParams {}; + struct ThermalConductionLawParams {}; + + OPM_HOST_DEVICE SolidEnergyLawParams solidEnergyLawParams(unsigned /*elemIdx*/) const + { + return SolidEnergyLawParams{}; + } + OPM_HOST_DEVICE ThermalConductionLawParams + thermalConductionLawParams(unsigned /*elemIdx*/) const + { + return ThermalConductionLawParams{}; + } +}; + +} // namespace Opm + +namespace Opm::gpuistl { + +/*! \brief copy_to_gpu overload for the no-op thermal manager. */ +inline ::Opm::NoThermalLawManager +copy_to_gpu(const ::Opm::NoThermalLawManager& /*cpu*/) +{ + return ::Opm::NoThermalLawManager{}; +} + +/*! \brief make_view overload for the no-op thermal manager. */ +inline ::Opm::NoThermalLawManager +make_view(::Opm::NoThermalLawManager& /*buf*/) +{ + return ::Opm::NoThermalLawManager{}; +} + +} // namespace Opm::gpuistl + +namespace Opm { + +/*! + * \brief Minimal, GPU-compatible problem class. + * + * Provides exactly the interface required by BlackOilIntensiveQuantities to + * compute a primary-variable update on the device. All cell-wise rock / + * fluid bookkeeping data is kept in templated Storage so the same class can + * be used as a CPU object (default \c std::vector) or as a GPU object + * (\c GpuBuffer / \c GpuView). + * + * The material law manager is held by value (so it is trivially copyable + * into a CUDA kernel together with the rest of the problem state). + * + * \tparam ScalarT Floating point type used for rock data. + * \tparam MaterialLawManagerT GPU material law manager type (e.g. + * \c Opm::EclMaterialLaw::GpuManager). + * \tparam Storage Container template used for per-cell data. + * \tparam ThermalLawManagerT GPU thermal-law manager type (e.g. + * \c Opm::EclThermalLaw::GpuManager). Defaults + * to \c NoThermalLawManager, a no-op stub used + * by non-thermal callers. + */ +template class Storage = ::Opm::VectorWithDefaultAllocator, + class ThermalLawManagerT = ::Opm::NoThermalLawManager> +class GpuFlowProblem +{ +public: + using Scalar = ScalarT; + using EclMaterialLawManager = MaterialLawManagerT; + using MaterialLawParams = typename EclMaterialLawManager::MaterialLawParams; + using EclThermalLawManager = ThermalLawManagerT; + using SolidEnergyLawParams = typename EclThermalLawManager::SolidEnergyLawParams; + using ThermalConductionLawParams = typename EclThermalLawManager::ThermalConductionLawParams; + + /*! \brief True iff the thermal-law manager carries a real fluid-system + * tag (anything other than \c void) and therefore the + * thermal-extraction code paths should be enabled. */ + static constexpr bool hasThermal + = !std::is_same_v; + + /*! \brief Trivial nested model() helper that satisfies the + * \c problem.model().linearizer().getLinearizationType() chain + * that BlackOilIntensiveQuantities walks. */ + struct ModelView + { + struct LinearizerView + { + OPM_HOST_DEVICE LinearizationType getLinearizationType() const + { + return LinearizationType{}; + } + }; + + OPM_HOST_DEVICE LinearizerView linearizer() const + { + return LinearizerView{}; + } + }; + + OPM_HOST_DEVICE GpuFlowProblem() = default; + + OPM_HOST_DEVICE GpuFlowProblem(EclMaterialLawManager materialLawManager, + Storage porosity, + Storage rockCompressibility, + Storage rockReferencePressure, + Storage maxOilSaturation, + Storage maxOilVaporizationFactor, + Storage maxGasDissolutionFactor, + EclThermalLawManager thermalLawManager, + Storage rockFraction, + Storage pvtRegionIndex) + : materialLawManager_(std::move(materialLawManager)) + , porosity_(std::move(porosity)) + , rockCompressibility_(std::move(rockCompressibility)) + , rockReferencePressure_(std::move(rockReferencePressure)) + , maxOilSaturation_(std::move(maxOilSaturation)) + , maxOilVaporizationFactor_(std::move(maxOilVaporizationFactor)) + , maxGasDissolutionFactor_(std::move(maxGasDissolutionFactor)) + , thermalLawManager_(std::move(thermalLawManager)) + , rockFraction_(std::move(rockFraction)) + , pvtRegionIndex_(std::move(pvtRegionIndex)) + { + } + + /*! + * \brief Construct from any CPU \c FlowProblem (or \c FlowProblemBlackoil). + * + * Iterates over all grid degrees of freedom and copies the rock / + * mixing-control fields used by \c BlackOilIntensiveQuantities into local + * vectors. The inner material-law manager is constructed from + * \c *cpu.materialLawManager(), so the GPU manager type held by this + * problem must in turn be constructible from + * \c Opm::EclMaterialLaw::Manager (i.e. use CPU vector storage and + * piecewise-linear vector-storage two-phase laws). + * + * Only enabled when this problem itself uses CPU vector storage. + */ + template , + std::enable_if_t< + std::is_same_v>, + int> = 0> + explicit GpuFlowProblem(const CpuProblem& cpu) + : materialLawManager_(*cpu.materialLawManager(), cpu.model().numGridDof()) + { + const std::size_t n = cpu.model().numGridDof(); + porosity_.resize(n); + rockCompressibility_.resize(n); + rockReferencePressure_.resize(n); + maxOilSaturation_.resize(n); + maxOilVaporizationFactor_.resize(n); + maxGasDissolutionFactor_.resize(n); + for (std::size_t i = 0; i < n; ++i) { + const unsigned u = static_cast(i); + porosity_[i] = cpu.porosity(u, 0u); + rockCompressibility_[i] = cpu.rockCompressibility(u); + rockReferencePressure_[i] = cpu.rockReferencePressure(u); + maxOilSaturation_[i] = cpu.maxOilSaturation(u); + maxOilVaporizationFactor_[i] = cpu.maxOilVaporizationFactor(0u, u); + maxGasDissolutionFactor_[i] = cpu.maxGasDissolutionFactor(0u, u); + } + if constexpr (hasThermal) { + using FluidSystemTag = typename EclThermalLawManager::FluidSystem; + thermalLawManager_ + = ::Opm::EclThermalLaw::buildCpuManagerFromFlowProblem( + cpu, n); + rockFraction_.resize(n); + pvtRegionIndex_.resize(n); + for (std::size_t i = 0; i < n; ++i) { + const unsigned u = static_cast(i); + rockFraction_[i] = cpu.rockFraction(u, 0u); + pvtRegionIndex_[i] = static_cast(cpu.pvtRegionIndex(u)); + } + } + } + + /*! + * \brief Construct directly into device-resident \c GpuBuffer storage + * from any CPU \c FlowProblem (or \c FlowProblemBlackoil). + * + * Internally builds the per-cell rock/mixing data on the host as + * \c std::vector, then uploads each as a single + * \c GpuBuffer. The inner GPU material law manager is constructed + * from \c *cpu.materialLawManager() via its own GpuBuffer-storage + * constructor, which uploads the per-cell sample arrays. + * + * Only enabled when this problem itself uses \c GpuBuffer storage. + */ + template , + std::enable_if_t< + std::is_same_v>, + int> = 0> + explicit GpuFlowProblem(const CpuProblem& cpu) + : materialLawManager_(*cpu.materialLawManager(), cpu.model().numGridDof()) + , porosity_(extractRockField(cpu, [](const CpuProblem& p, unsigned u) { + return Scalar(p.porosity(u, 0u)); + })) + , rockCompressibility_(extractRockField(cpu, [](const CpuProblem& p, unsigned u) { + return Scalar(p.rockCompressibility(u)); + })) + , rockReferencePressure_(extractRockField(cpu, [](const CpuProblem& p, unsigned u) { + return Scalar(p.rockReferencePressure(u)); + })) + , maxOilSaturation_(extractRockField(cpu, [](const CpuProblem& p, unsigned u) { + return Scalar(p.maxOilSaturation(u)); + })) + , maxOilVaporizationFactor_(extractRockField(cpu, [](const CpuProblem& p, unsigned u) { + return Scalar(p.maxOilVaporizationFactor(0u, u)); + })) + , maxGasDissolutionFactor_(extractRockField(cpu, [](const CpuProblem& p, unsigned u) { + return Scalar(p.maxGasDissolutionFactor(0u, u)); + })) + , thermalLawManager_(buildGpuThermalManager(cpu)) + , rockFraction_(extractThermalScalarField(cpu, [](const CpuProblem& p, unsigned u) { + return Scalar(p.rockFraction(u, 0u)); + })) + , pvtRegionIndex_(extractThermalIntField(cpu, [](const CpuProblem& p, unsigned u) { + return static_cast(p.pvtRegionIndex(u)); + })) + { + } + + OPM_HOST_DEVICE ModelView model() const + { + return ModelView{}; + } + + OPM_HOST_DEVICE int satnumRegionIndex(std::size_t elemIdx) const + { + return materialLawManager_.satnumRegionIdx(static_cast(elemIdx)); + } + + OPM_HOST_DEVICE MaterialLawParams materialLawParams(std::size_t elemIdx) const + { + return materialLawManager_.materialLawParams(static_cast(elemIdx)); + } + + OPM_HOST_DEVICE Scalar rockCompressibility(std::size_t elemIdx) const + { + return rockCompressibility_.size() == 0 ? Scalar(0) : rockCompressibility_[elemIdx]; + } + + OPM_HOST_DEVICE Scalar rockReferencePressure(std::size_t elemIdx) const + { + return rockReferencePressure_.size() == 0 ? Scalar(0) : rockReferencePressure_[elemIdx]; + } + + OPM_HOST_DEVICE Scalar porosity(std::size_t elemIdx, unsigned /*timeIdx*/) const + { + return porosity_.size() == 0 ? Scalar(0) : porosity_[elemIdx]; + } + + OPM_HOST_DEVICE Scalar maxOilVaporizationFactor(unsigned /*timeIdx*/, std::size_t elemIdx) const + { + return maxOilVaporizationFactor_.size() == 0 ? Scalar(0) : maxOilVaporizationFactor_[elemIdx]; + } + + OPM_HOST_DEVICE Scalar maxGasDissolutionFactor(unsigned /*timeIdx*/, std::size_t elemIdx) const + { + return maxGasDissolutionFactor_.size() == 0 ? Scalar(0) : maxGasDissolutionFactor_[elemIdx]; + } + + OPM_HOST_DEVICE Scalar maxOilSaturation(std::size_t elemIdx) const + { + return maxOilSaturation_.size() == 0 ? Scalar(0) : maxOilSaturation_[elemIdx]; + } + + /*! \brief Per-cell PVT region index. Returns 0 when the GpuFlowProblem + * instantiation has no thermal support. */ + OPM_HOST_DEVICE unsigned pvtRegionIndex(std::size_t elemIdx) const + { + return pvtRegionIndex_.size() == 0 + ? 0u + : static_cast(pvtRegionIndex_[elemIdx]); + } + + /*! \brief Per-cell rock fraction (1 - effective porosity). Returns 0 + * when the GpuFlowProblem instantiation has no thermal + * support. */ + OPM_HOST_DEVICE Scalar rockFraction(std::size_t elemIdx, unsigned /*timeIdx*/) const + { + return rockFraction_.size() == 0 ? Scalar(0) : rockFraction_[elemIdx]; + } + + /*! \brief Solid-energy law parameters for a single cell. Forwards to + * the embedded thermal-law manager. */ + OPM_HOST_DEVICE SolidEnergyLawParams + solidEnergyLawParams(std::size_t elemIdx, unsigned /*timeIdx*/) const + { + return thermalLawManager_.solidEnergyLawParams(static_cast(elemIdx)); + } + + /*! \brief Thermal-conduction law parameters for a single cell. + * Forwards to the embedded thermal-law manager. */ + OPM_HOST_DEVICE ThermalConductionLawParams + thermalConductionLawParams(std::size_t elemIdx, unsigned /*timeIdx*/) const + { + return thermalLawManager_.thermalConductionLawParams(static_cast(elemIdx)); + } + + /*! \brief Default rock-pore-volume multiplier (1 if no compressibility). */ + template + OPM_HOST_DEVICE Evaluation rockCompPoroMultiplier(const auto& /*intQuants*/, std::size_t /*elemIdx*/) const + { + return Evaluation(1.0); + } + + /*! \brief Default rock-trans multiplier (1 if no compressibility). */ + template + OPM_HOST_DEVICE Evaluation rockCompTransMultiplier(const auto& /*intQuants*/, std::size_t /*elemIdx*/) const + { + return Evaluation(1.0); + } + + /*! + * \brief Update the relative permeabilities of all phases for a single + * cell, in the same way as the CPU \c FlowProblem does. + * + * Calls \c MaterialLaw::relativePermeabilities on the cell's + * \c MaterialLawParams. The \c BlackOilIntensiveQuantities later divides + * the result by the phase viscosity to obtain the mobility. + * + * Directional relative permeabilities are not supported by the GPU + * problem, so the \p dirMob output is left untouched. + */ + template + OPM_HOST_DEVICE void updateRelperms(auto& mobility, + auto& /*dirMob*/, + FluidState& fluidState, + unsigned globalSpaceIdx) const + { + using ContainerT = std::decay_t; + const auto materialParams = materialLawParams(globalSpaceIdx); + EclMaterialLawManager::MaterialLaw::template relativePermeabilities< + ContainerT, FluidState, Args...>(mobility, materialParams, fluidState); + } + + /*! \name Direct accessors for serialization to the GPU. */ + //!\{ + const EclMaterialLawManager& materialLawManager() const { return materialLawManager_; } + EclMaterialLawManager& materialLawManager() { return materialLawManager_; } + const Storage& porosityStorage() const { return porosity_; } + Storage& porosityStorage() { return porosity_; } + const Storage& rockCompressibilityStorage() const { return rockCompressibility_; } + Storage& rockCompressibilityStorage() { return rockCompressibility_; } + const Storage& rockReferencePressureStorage() const { return rockReferencePressure_; } + Storage& rockReferencePressureStorage() { return rockReferencePressure_; } + const Storage& maxOilSaturationStorage() const { return maxOilSaturation_; } + Storage& maxOilSaturationStorage() { return maxOilSaturation_; } + const Storage& maxOilVaporizationFactorStorage() const { return maxOilVaporizationFactor_; } + Storage& maxOilVaporizationFactorStorage() { return maxOilVaporizationFactor_; } + const Storage& maxGasDissolutionFactorStorage() const { return maxGasDissolutionFactor_; } + Storage& maxGasDissolutionFactorStorage() { return maxGasDissolutionFactor_; } + const EclThermalLawManager& thermalLawManager() const { return thermalLawManager_; } + EclThermalLawManager& thermalLawManager() { return thermalLawManager_; } + const Storage& rockFractionStorage() const { return rockFraction_; } + Storage& rockFractionStorage() { return rockFraction_; } + const Storage& pvtRegionIndexStorage() const { return pvtRegionIndex_; } + Storage& pvtRegionIndexStorage() { return pvtRegionIndex_; } + //!\} + +private: + template + static Storage extractRockField(const CpuProblem& cpu, F f) + { + const std::size_t n = cpu.model().numGridDof(); + std::vector v(n); + for (std::size_t i = 0; i < n; ++i) { + v[i] = f(cpu, static_cast(i)); + } + if constexpr (std::is_same_v, + ::Opm::VectorWithDefaultAllocator>) { + return Storage(v.begin(), v.end()); + } else { + return Storage(v); + } + } + + /*! \brief Build the per-cell rockFraction storage; empty when the + * GpuFlowProblem instantiation has no thermal support. */ + template + static Storage extractThermalScalarField(const CpuProblem& cpu, F f) + { + if constexpr (!hasThermal) { + return Storage{}; + } else { + return extractRockField(cpu, f); + } + } + + /*! \brief Build the per-cell pvtRegionIndex storage; empty when the + * GpuFlowProblem instantiation has no thermal support. */ + template + static Storage extractThermalIntField(const CpuProblem& cpu, F f) + { + if constexpr (!hasThermal) { + return Storage{}; + } else { + const std::size_t n = cpu.model().numGridDof(); + std::vector v(n); + for (std::size_t i = 0; i < n; ++i) { + v[i] = f(cpu, static_cast(i)); + } + if constexpr (std::is_same_v, + ::Opm::VectorWithDefaultAllocator>) { + return Storage(v.begin(), v.end()); + } else { + return Storage(v); + } + } + } + + /*! \brief Build the GPU thermal-law manager from a CPU FlowProblem. + * Returns a default-constructed (empty) manager when this + * GpuFlowProblem instantiation has no thermal support. */ + template + static EclThermalLawManager buildGpuThermalManager(const CpuProblem& cpu) + { + if constexpr (!hasThermal) { + return EclThermalLawManager{}; + } else { + using FluidSystemTag = typename EclThermalLawManager::FluidSystem; + auto cpuMgr + = ::Opm::EclThermalLaw::buildCpuManagerFromFlowProblem( + cpu, cpu.model().numGridDof()); + return ::Opm::gpuistl::copy_to_gpu(cpuMgr); + } + } + + EclMaterialLawManager materialLawManager_{}; + Storage porosity_{}; + Storage rockCompressibility_{}; + Storage rockReferencePressure_{}; + Storage maxOilSaturation_{}; + Storage maxOilVaporizationFactor_{}; + Storage maxGasDissolutionFactor_{}; + EclThermalLawManager thermalLawManager_{}; + Storage rockFraction_{}; + Storage pvtRegionIndex_{}; +}; + +} // namespace Opm + +namespace Opm::gpuistl { + +/*! + * \brief Copy a CPU GpuFlowProblem to GPU-resident GpuBuffer storage. + * + * The material law manager is also moved to the GPU (via copy_to_gpu). + */ +template +auto copy_to_gpu( + const ::Opm::GpuFlowProblem& cpu) +{ + using GpuMaterialLawManagerBuffer + = decltype(::Opm::gpuistl::copy_to_gpu(cpu.materialLawManager())); + using GpuThermalLawManagerBuffer + = decltype(::Opm::gpuistl::copy_to_gpu(cpu.thermalLawManager())); + using GpuProblemType = ::Opm::GpuFlowProblem; + + return GpuProblemType(::Opm::gpuistl::copy_to_gpu(cpu.materialLawManager()), + GpuBuffer(cpu.porosityStorage()), + GpuBuffer(cpu.rockCompressibilityStorage()), + GpuBuffer(cpu.rockReferencePressureStorage()), + GpuBuffer(cpu.maxOilSaturationStorage()), + GpuBuffer(cpu.maxOilVaporizationFactorStorage()), + GpuBuffer(cpu.maxGasDissolutionFactorStorage()), + ::Opm::gpuistl::copy_to_gpu(cpu.thermalLawManager()), + GpuBuffer(cpu.rockFractionStorage()), + GpuBuffer(cpu.pvtRegionIndexStorage())); +} + +/*! + * \brief Make a non-owning GpuView based GpuFlowProblem from an owning + * GpuBuffer based GpuFlowProblem. + */ +template +auto make_view(::Opm::GpuFlowProblem& buf) +{ + using GpuMaterialLawManagerView + = decltype(::Opm::gpuistl::make_view(buf.materialLawManager())); + using GpuThermalLawManagerView + = decltype(::Opm::gpuistl::make_view(buf.thermalLawManager())); + using GpuProblemView = ::Opm::GpuFlowProblem; + + auto toView = [](auto& storage) { + using T = typename std::decay_t::value_type; + return GpuView(storage.data(), storage.size()); + }; + + return GpuProblemView(::Opm::gpuistl::make_view(buf.materialLawManager()), + toView(buf.porosityStorage()), + toView(buf.rockCompressibilityStorage()), + toView(buf.rockReferencePressureStorage()), + toView(buf.maxOilSaturationStorage()), + toView(buf.maxOilVaporizationFactorStorage()), + toView(buf.maxGasDissolutionFactorStorage()), + ::Opm::gpuistl::make_view(buf.thermalLawManager()), + toView(buf.rockFractionStorage()), + toView(buf.pvtRegionIndexStorage())); +} + +} // namespace Opm::gpuistl + +#endif // OPM_GPU_FLOW_PROBLEM_HPP diff --git a/tests/gpuistl/test_gpu_ecl_thermal_law_manager.cu b/tests/gpuistl/test_gpu_ecl_thermal_law_manager.cu new file mode 100644 index 00000000000..534e9d1abc2 --- /dev/null +++ b/tests/gpuistl/test_gpu_ecl_thermal_law_manager.cu @@ -0,0 +1,210 @@ +/* + 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 . +*/ +#include + +#define BOOST_TEST_MODULE TestGpuEclThermalLawManager + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +namespace { + +using Scalar = double; + +// A minimal stand-in fluid system. The thermal-conduction law only ever +// asks the fluid system whether the gas phase is active and what its +// phase index is. Both are compile-time constants. +struct DummyFluidSystem { + static constexpr int numPhases = 3; + static constexpr int gasPhaseIdx = 2; + OPM_HOST_DEVICE static constexpr bool phaseIsActive(int phaseIdx) + { + return phaseIdx == gasPhaseIdx; + } +}; + +// A minimal stand-in fluid state with just enough surface area to feed +// EclSpecrockLaw / EclThconrLaw / their GPU counterparts. +template +struct DummyFluidState { + using ValueType = ScalarT; + ScalarT temperature_ = 0; + ScalarT gasSaturation_ = 0; + + OPM_HOST_DEVICE ScalarT temperature(int /*phaseIdx*/) const { return temperature_; } + OPM_HOST_DEVICE ScalarT saturation(int /*phaseIdx*/) const { return gasSaturation_; } + OPM_HOST_DEVICE DummyFluidSystem fluidSystem() const { return DummyFluidSystem{}; } +}; + +using GpuSpecrockParamsBuf = Opm::EclSpecrockLawParams; +using GpuSpecrockParamsView = Opm::EclSpecrockLawParams; +using GpuSpecrockLawView = Opm::EclSpecrockLaw; +using GpuThconrLaw = Opm::EclThconrLaw; + +using ManagerCpu = Opm::EclThermalLaw::GpuManager; +using ManagerBuf + = Opm::EclThermalLaw::GpuManager; +using ManagerView + = Opm::EclThermalLaw::GpuManager; + +// Kernel: per-cell, evaluate solid internal energy from the law manager +// at the supplied temperature, and total thermal conductivity at the +// supplied gas saturation. +__global__ void evaluateThermalKernel(ManagerView view, + const Scalar* temperatures, + const Scalar* gasSaturations, + Scalar* outSolidEnergy, + Scalar* outConductivity, + std::size_t numCells) +{ + const std::size_t i = + static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (i >= numCells) { + return; + } + DummyFluidState fs; + fs.temperature_ = temperatures[i]; + fs.gasSaturation_ = gasSaturations[i]; + + const auto solidEnergyParams = view.solidEnergyLawParams(static_cast(i)); + outSolidEnergy[i] = + GpuSpecrockLawView::solidInternalEnergy(solidEnergyParams, fs); + + const auto thermalConductionParams = + view.thermalConductionLawParams(static_cast(i)); + outConductivity[i] = + GpuThconrLaw::thermalConductivity(thermalConductionParams, fs); +} + +} // namespace + +BOOST_AUTO_TEST_CASE(SpecrockAndThconrEvaluateMatchesCpu) +{ + // -- Build a synthetic two-region SPECROCK plus per-cell THCONR setup + // on the host. Cells 0..2 belong to region 0, cells 3..5 to region 1. + constexpr std::size_t numCells = 6; + + const std::vector regionTemperatures0 = { 280.0, 320.0, 360.0, 400.0 }; + const std::vector regionEnergies0 = { 0.0, 4.0e7, 9.0e7, 1.5e8 }; + const std::vector regionTemperatures1 = { 280.0, 350.0, 420.0 }; + const std::vector regionEnergies1 = { 0.0, 3.5e7, 8.4e7 }; + + // Build a host-side CPU manager directly from the basic constructor. + Opm::EclSpecrockLawParams region0Params; + region0Params.setSamples(regionTemperatures0, regionEnergies0); + Opm::EclSpecrockLawParams region1Params; + region1Params.setSamples(regionTemperatures1, regionEnergies1); + + std::vector> cpuSolidEnergyParams; + cpuSolidEnergyParams.push_back(region0Params); + cpuSolidEnergyParams.push_back(region1Params); + + std::vector cpuElementToRegion(numCells); + for (std::size_t i = 0; i < numCells; ++i) { + cpuElementToRegion[i] = (i < 3) ? 0 : 1; + } + + std::vector> cpuThermalConductionParams(numCells); + for (std::size_t i = 0; i < numCells; ++i) { + cpuThermalConductionParams[i].setReferenceTotalThermalConductivity( + 100.0 + 5.0 * static_cast(i)); + cpuThermalConductionParams[i].setDTotalThermalConductivity_dSg( + 0.20 + 0.01 * static_cast(i)); + cpuThermalConductionParams[i].finalize(); + } + + ManagerCpu cpuManager(std::move(cpuSolidEnergyParams), + std::move(cpuElementToRegion), + std::move(cpuThermalConductionParams)); + + // -- Upload to the GPU and create a non-owning view. + auto bufManager = Opm::gpuistl::copy_to_gpu(cpuManager); + auto viewManager = Opm::gpuistl::make_view(bufManager); + + // -- Pick a per-cell test temperature and gas saturation each. + std::array hostTemperatures = { + 290.0, 310.0, 355.0, 285.0, 360.0, 410.0 + }; + std::array hostGasSaturations = { + 0.0, 0.1, 0.25, 0.5, 0.75, 0.99 + }; + + // -- Run the kernel. + auto dTOwner = Opm::gpuistl::make_gpu_unique_ptr_array(numCells); + auto dSgOwner = Opm::gpuistl::make_gpu_unique_ptr_array(numCells); + auto dEnergyOwner = Opm::gpuistl::make_gpu_unique_ptr_array(numCells); + auto dConductivityOwner = Opm::gpuistl::make_gpu_unique_ptr_array(numCells); + Scalar* dT = dTOwner.get(); + Scalar* dSg = dSgOwner.get(); + Scalar* dEnergy = dEnergyOwner.get(); + Scalar* dConductivity = dConductivityOwner.get(); + OPM_GPU_SAFE_CALL(cudaMemcpy(dT, hostTemperatures.data(), + numCells * sizeof(Scalar), cudaMemcpyHostToDevice)); + OPM_GPU_SAFE_CALL(cudaMemcpy(dSg, hostGasSaturations.data(), + numCells * sizeof(Scalar), cudaMemcpyHostToDevice)); + + evaluateThermalKernel<<<1, 32>>>(viewManager, dT, dSg, dEnergy, dConductivity, numCells); + OPM_GPU_SAFE_CALL(cudaDeviceSynchronize()); + OPM_GPU_SAFE_CALL(cudaGetLastError()); + + std::array hostEnergyResult {}; + std::array hostConductivityResult {}; + OPM_GPU_SAFE_CALL(cudaMemcpy(hostEnergyResult.data(), dEnergy, + numCells * sizeof(Scalar), cudaMemcpyDeviceToHost)); + OPM_GPU_SAFE_CALL(cudaMemcpy(hostConductivityResult.data(), dConductivity, + numCells * sizeof(Scalar), cudaMemcpyDeviceToHost)); + + // -- Compare against the CPU manager. + for (std::size_t i = 0; i < numCells; ++i) { + DummyFluidState fs; + fs.temperature_ = hostTemperatures[i]; + fs.gasSaturation_ = hostGasSaturations[i]; + + const auto& cpuSolidParams = + cpuManager.solidEnergyLawParams(static_cast(i)); + const Scalar cpuEnergy = + Opm::EclSpecrockLaw::solidInternalEnergy(cpuSolidParams, fs); + BOOST_CHECK_CLOSE(cpuEnergy, hostEnergyResult[i], 1e-9); + + const auto& cpuConductionParams = + cpuManager.thermalConductionLawParams(static_cast(i)); + const Scalar cpuConductivity = + GpuThconrLaw::thermalConductivity(cpuConductionParams, fs); + BOOST_CHECK_CLOSE(cpuConductivity, hostConductivityResult[i], 1e-9); + } +} diff --git a/tests/very_simple_deck.DATA b/tests/very_simple_deck.DATA new file mode 100644 index 00000000000..8c0623914ab --- /dev/null +++ b/tests/very_simple_deck.DATA @@ -0,0 +1,53 @@ +-- =============== RUNSPEC +RUNSPEC +DIMENS +3 3 3 / +EQLDIMS +/ +TABDIMS +/ +WATER +GAS +CO2STORE +METRIC + +-- =============== GRID +GRID +GRIDFILE +0 0 / +DX +27*1 / +DY +27*1 / +DZ +27*1 / +TOPS +9*0 / +PERMX +27*1013.25 / +PORO +27*0.25 / +COPY +PERMX PERMY / +PERMX PERMZ / +/ + +-- =============== PROPS +PROPS +SGWFN +0.000000E+00 0.000000E+00 1.000000E+00 3.060000E-02 +1.000000E+00 1.000000E+00 0.000000E+00 3.060000E-01 / + +-- =============== SOLUTION +SOLUTION +RPTRST +'BASIC=0' / +EQUIL +0 300 100 0 0 0 1 1 0 / + +-- =============== SCHEDULE +SCHEDULE +RPTRST +'BASIC=0' / +TSTEP +1 /