GpuFlowProblem.hpp

Go to the documentation of this file.

// -*- 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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_GPU_FLOW_PROBLEM_HPP
#define OPM_GPU_FLOW_PROBLEM_HPP
 
#include <opm/common/utility/VectorWithDefaultAllocator.hpp>
#include <opm/common/utility/gpuDecorators.hpp>
#include <opm/common/ErrorMacros.hpp>
 
#include <opm/models/discretization/common/linearizationtype.hh>
 
#include <opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp>
#include <opm/material/thermal/GpuEclThermalLawManager.hpp>
 
#include <concepts>
#include <cstddef>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace Opm::gpuistl
{
template <typename T>
class GpuBuffer;
template <typename T>
class GpuView;
} // namespace Opm::gpuistl
 
namespace Opm {
 
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 {
 
inline ::Opm::NoThermalLawManager
copy_to_gpu(const ::Opm::NoThermalLawManager& /*cpu*/)
{
    return ::Opm::NoThermalLawManager{};
}
 
inline ::Opm::NoThermalLawManager
make_view(::Opm::NoThermalLawManager& /*buf*/)
{
    return ::Opm::NoThermalLawManager{};
}
 
} // namespace Opm::gpuistl
 
namespace Opm {
 
template <class ScalarT,
          class MaterialLawManagerT,
          template <class> 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;
 
    static constexpr bool hasThermal
        = !std::is_same_v<typename EclThermalLawManager::FluidSystem, void>;
 
    struct ModelView
    {
        struct LinearizerView
        {
            OPM_HOST_DEVICE LinearizationType getLinearizationType() const
            {
                return LinearizationType{};
            }
        };
 
        OPM_HOST_DEVICE LinearizerView linearizer() const
        {
            return LinearizerView{};
        }
    };
 
    GpuFlowProblem() = default;
 
    GpuFlowProblem(EclMaterialLawManager materialLawManager,
                   Storage<Scalar> porosity,
                   Storage<Scalar> rockCompressibility,
                   Storage<Scalar> rockReferencePressure,
                   Storage<Scalar> maxOilSaturation,
                   Storage<Scalar> maxOilVaporizationFactor,
                   Storage<Scalar> maxGasDissolutionFactor,
                   EclThermalLawManager thermalLawManager,
                   Storage<Scalar> rockFraction,
                   Storage<int> pvtRegionIndex,
                   bool usesDefaultRockCompaction)
        : 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))
    {
        if (!usesDefaultRockCompaction) {
            OPM_THROW(std::logic_error,
                      "GPU FlowProblem does not support non-default rock compaction multipliers");
        }
    }
 
    template <class CpuProblem>
        requires std::same_as<Storage<Scalar>, ::Opm::VectorWithDefaultAllocator<Scalar>>
    explicit GpuFlowProblem(const CpuProblem& cpu)
        : materialLawManager_(cpu.materialLawManager(), cpu.model().numGridDof())
    {
        if (!cpu.usesDefaultRockCompaction()) {
            OPM_THROW(std::logic_error,
                      "GPU FlowProblem does not support non-default rock compaction multipliers");
        }
        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<unsigned>(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<Scalar, FluidSystemTag>(
                    cpu, n);
            rockFraction_.resize(n);
            pvtRegionIndex_.resize(n);
            for (std::size_t i = 0; i < n; ++i) {
                const unsigned u = static_cast<unsigned>(i);
                rockFraction_[i] = cpu.rockFraction(u, 0u);
                pvtRegionIndex_[i] = static_cast<int>(cpu.pvtRegionIndex(u));
            }
        }
    }
 
    template <class CpuProblem>
        requires std::same_as<Storage<Scalar>, ::Opm::gpuistl::GpuBuffer<Scalar>>
    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<int>(p.pvtRegionIndex(u));
          }))
    {
        if (!cpu.usesDefaultRockCompaction()) {
            OPM_THROW(std::logic_error,
                      "GPU FlowProblem does not support non-default rock compaction multipliers");
        }
    }
 
    OPM_HOST_DEVICE ModelView model() const
    {
        return ModelView{};
    }
 
    OPM_HOST_DEVICE int satnumRegionIndex(std::size_t elemIdx) const
    {
        return materialLawManager_.satnumRegionIdx(static_cast<unsigned>(elemIdx));
    }
 
    OPM_HOST_DEVICE MaterialLawParams materialLawParams(std::size_t elemIdx) const
    {
        return materialLawManager_.materialLawParams(static_cast<unsigned>(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];
    }
 
    OPM_HOST_DEVICE unsigned pvtRegionIndex(std::size_t elemIdx) const
    {
        return pvtRegionIndex_.size() == 0
                   ? 0u
                   : static_cast<unsigned>(pvtRegionIndex_[elemIdx]);
    }
 
    OPM_HOST_DEVICE Scalar rockFraction(std::size_t elemIdx, unsigned /*timeIdx*/) const
    {
        return rockFraction_.size() == 0 ? Scalar(0) : rockFraction_[elemIdx];
    }
 
    OPM_HOST_DEVICE SolidEnergyLawParams
    solidEnergyLawParams(std::size_t elemIdx, unsigned /*timeIdx*/) const
    {
        return thermalLawManager_.solidEnergyLawParams(static_cast<unsigned>(elemIdx));
    }
 
    OPM_HOST_DEVICE ThermalConductionLawParams
    thermalConductionLawParams(std::size_t elemIdx, unsigned /*timeIdx*/) const
    {
        return thermalLawManager_.thermalConductionLawParams(static_cast<unsigned>(elemIdx));
    }
 
    template <class Evaluation>
    OPM_HOST_DEVICE Evaluation rockCompPoroMultiplier(const auto& /*intQuants*/, std::size_t /*elemIdx*/) const
    {
        // Ctor guarantees no rock compaction, this is therefore always correct
        return Evaluation(1.0);
    }
 
    template <class Evaluation>
    OPM_HOST_DEVICE Evaluation rockCompTransMultiplier(const auto& /*intQuants*/, std::size_t /*elemIdx*/) const
    {
        // Ctor guarantees no rock compaction, this is therefore always correct
        return Evaluation(1.0);
    }
 
    template <class FluidState, class... Args>
    OPM_HOST_DEVICE void updateRelperms(auto& mobility,
                                        auto& /*dirMob*/,
                                        FluidState& fluidState,
                                        unsigned globalSpaceIdx) const
    {
        using ContainerT = std::decay_t<decltype(mobility)>;
        const auto materialParams = materialLawParams(globalSpaceIdx);
        EclMaterialLawManager::MaterialLaw::template relativePermeabilities<
            ContainerT, FluidState, Args...>(mobility, materialParams, fluidState);
    }
 
    const EclMaterialLawManager& materialLawManager() const
    { return materialLawManager_; }
 
    EclMaterialLawManager& materialLawManager()
    { return materialLawManager_; }
 
    const Storage<Scalar>& porosityStorage() const
    { return porosity_; }
 
    Storage<Scalar>& porosityStorage()
    { return porosity_; }
 
    const Storage<Scalar>& rockCompressibilityStorage() const
    { return rockCompressibility_; }
 
    Storage<Scalar>& rockCompressibilityStorage()
    { return rockCompressibility_; }
 
    const Storage<Scalar>& rockReferencePressureStorage() const
    { return rockReferencePressure_; }
 
    Storage<Scalar>& rockReferencePressureStorage()
    { return rockReferencePressure_; }
 
    const Storage<Scalar>& maxOilSaturationStorage() const
    { return maxOilSaturation_; }
 
    Storage<Scalar>& maxOilSaturationStorage()
    { return maxOilSaturation_; }
 
    const Storage<Scalar>& maxOilVaporizationFactorStorage() const
    { return maxOilVaporizationFactor_; }
 
    Storage<Scalar>& maxOilVaporizationFactorStorage()
    { return maxOilVaporizationFactor_; }
 
    const Storage<Scalar>& maxGasDissolutionFactorStorage() const
    { return maxGasDissolutionFactor_; }
 
    Storage<Scalar>& maxGasDissolutionFactorStorage()
    { return maxGasDissolutionFactor_; }
 
    const EclThermalLawManager& thermalLawManager() const
    { return thermalLawManager_; }
 
    EclThermalLawManager& thermalLawManager()
    { return thermalLawManager_; }
 
    const Storage<Scalar>& rockFractionStorage() const
    { return rockFraction_; }
 
    Storage<Scalar>& rockFractionStorage()
    { return rockFraction_; }
 
    const Storage<int>& pvtRegionIndexStorage() const
    { return pvtRegionIndex_; }
 
    Storage<int>& pvtRegionIndexStorage()
    { return pvtRegionIndex_; }
 
private:
    template <class CpuProblem, class F>
    static Storage<Scalar> extractRockField(const CpuProblem& cpu, F f)
    {
        const std::size_t n = cpu.model().numGridDof();
        std::vector<Scalar> v(n);
        for (std::size_t i = 0; i < n; ++i) {
            v[i] = f(cpu, static_cast<unsigned>(i));
        }
        if constexpr (std::is_same_v<Storage<Scalar>,
                                     ::Opm::VectorWithDefaultAllocator<Scalar>>) {
            return Storage<Scalar>(v.begin(), v.end());
        } else {
            return Storage<Scalar>(v);
        }
    }
 
    template <class CpuProblem, class F>
    static Storage<Scalar> extractThermalScalarField(const CpuProblem& cpu, F f)
    {
        if constexpr (!hasThermal) {
            return Storage<Scalar>{};
        } else {
            return extractRockField(cpu, f);
        }
    }
 
    template <class CpuProblem, class F>
    static Storage<int> extractThermalIntField(const CpuProblem& cpu, F f)
    {
        if constexpr (!hasThermal) {
            return Storage<int>{};
        } else {
            const std::size_t n = cpu.model().numGridDof();
            std::vector<int> v(n);
            for (std::size_t i = 0; i < n; ++i) {
                v[i] = f(cpu, static_cast<unsigned>(i));
            }
            if constexpr (std::is_same_v<Storage<int>,
                                         ::Opm::VectorWithDefaultAllocator<int>>) {
                return Storage<int>(v.begin(), v.end());
            } else {
                return Storage<int>(v);
            }
        }
    }
 
    template <class CpuProblem>
    static EclThermalLawManager buildGpuThermalManager(const CpuProblem& cpu)
    {
        if constexpr (!hasThermal) {
            return EclThermalLawManager{};
        } else {
            using FluidSystemTag = typename EclThermalLawManager::FluidSystem;
            auto cpuMgr
                = ::Opm::EclThermalLaw::buildCpuManagerFromFlowProblem<Scalar, FluidSystemTag>(
                    cpu, cpu.model().numGridDof());
            return ::Opm::gpuistl::copy_to_gpu(cpuMgr);
        }
    }
 
    EclMaterialLawManager materialLawManager_{};
    Storage<Scalar> porosity_{};
    Storage<Scalar> rockCompressibility_{};
    Storage<Scalar> rockReferencePressure_{};
    Storage<Scalar> maxOilSaturation_{};
    Storage<Scalar> maxOilVaporizationFactor_{};
    Storage<Scalar> maxGasDissolutionFactor_{};
    EclThermalLawManager thermalLawManager_{};
    Storage<Scalar> rockFraction_{};
    Storage<int> pvtRegionIndex_{};
};
 
} // namespace Opm
 
namespace Opm::gpuistl {
 
template <class ScalarT, class CpuMaterialLawManager, class CpuThermalLawManager>
auto copy_to_gpu(
    const ::Opm::GpuFlowProblem<ScalarT,
                                CpuMaterialLawManager,
                                ::Opm::VectorWithDefaultAllocator,
                                CpuThermalLawManager>& 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<ScalarT,
                                                 GpuMaterialLawManagerBuffer,
                                                 GpuBuffer,
                                                 GpuThermalLawManagerBuffer>;
 
    return GpuProblemType(::Opm::gpuistl::copy_to_gpu(cpu.materialLawManager()),
                          GpuBuffer<ScalarT>(cpu.porosityStorage()),
                          GpuBuffer<ScalarT>(cpu.rockCompressibilityStorage()),
                          GpuBuffer<ScalarT>(cpu.rockReferencePressureStorage()),
                          GpuBuffer<ScalarT>(cpu.maxOilSaturationStorage()),
                          GpuBuffer<ScalarT>(cpu.maxOilVaporizationFactorStorage()),
                          GpuBuffer<ScalarT>(cpu.maxGasDissolutionFactorStorage()),
                          ::Opm::gpuistl::copy_to_gpu(cpu.thermalLawManager()),
                          GpuBuffer<ScalarT>(cpu.rockFractionStorage()),
                          GpuBuffer<int>(cpu.pvtRegionIndexStorage()),
                          true);
}
 
template <class ScalarT, class GpuBufferMaterialLawManager, class GpuBufferThermalLawManager>
auto make_view(::Opm::GpuFlowProblem<ScalarT,
                                     GpuBufferMaterialLawManager,
                                     GpuBuffer,
                                     GpuBufferThermalLawManager>& buf)
{
    using GpuMaterialLawManagerView
        = decltype(::Opm::gpuistl::make_view(buf.materialLawManager()));
    using GpuThermalLawManagerView
        = decltype(::Opm::gpuistl::make_view(buf.thermalLawManager()));
    using GpuProblemView = ::Opm::GpuFlowProblem<ScalarT,
                                                 GpuMaterialLawManagerView,
                                                 GpuView,
                                                 GpuThermalLawManagerView>;
 
    auto toView = [](auto& storage) {
        using T = typename std::decay_t<decltype(storage)>::value_type;
        return GpuView<T>(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()),
                          true);
}
 
} // namespace Opm::gpuistl
 
#endif // OPM_GPU_FLOW_PROBLEM_HPP