diff --git a/CMakeLists.txt b/CMakeLists.txt
index 976f72086d3..3eb30947ca0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -365,8 +365,10 @@ macro(opm-simulators_sources_hook)
     # cuda warns when constxpr functions are used in kernels.
     # since the entire stl is constexpr ..., we enable relaxed flag
     set_source_files_properties(
+      opm/simulators/linalg/gpuistl/GpuBlackoilIntensiveQuantitiesDispatcher.cu
       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
@@ -735,6 +737,7 @@ macro(opm-simulators_targets_hook)
       SOURCES
         flow/flow_${OBJ}.cpp
     )
+
     simulators_add_target_options(TARGET flow_lib${OBJ})
     list(APPEND FLOW_TGTS $<TARGET_OBJECTS:flow_lib${OBJ}>)
     list(APPEND FLOWMODELS_PREFIXED flow_${OBJ})
@@ -775,6 +778,58 @@ macro(opm-simulators_targets_hook)
       opmsimulators
   )
 
+  # Flow GPU executables handled separately to avoid making flow dependent on CUDA/HIP
+  set(FLOW_GPU_MODELS
+    gpu_gaswater_thermal
+  )
+
+  if(CONVERT_CUDA_TO_HIP OR CUDA_FOUND)
+
+    if(CONVERT_CUDA_TO_HIP)
+      set(EXTENSION "hip")
+    else()
+      set(EXTENSION "cu")
+    endif()
+
+    foreach(OBJ ${FLOW_GPU_MODELS})
+
+      opm_add_library(
+        TARGET
+          flow_lib${OBJ}
+        TYPE
+          OBJECT
+        SOURCES
+          flow/flow_${OBJ}.${EXTENSION}
+      )
+
+      simulators_add_target_options(TARGET flow_lib${OBJ})
+      list(APPEND FLOW_GPU_TGTS $<TARGET_OBJECTS:flow_lib${OBJ}>)
+      list(APPEND FLOWMODELS_PREFIXED flow_${OBJ})
+
+      opm_add_executable(
+        TARGET
+          flow_${OBJ}
+        SOURCES
+          flow/flow_${OBJ}_main.cpp
+          $<TARGET_OBJECTS:moduleVersion>
+          $<TARGET_OBJECTS:flow_lib${OBJ}>
+        LIBRARIES
+          opmsimulators
+      )
+    endforeach()
+
+    opm_add_executable(
+      TARGET
+        flow_gpu
+      SOURCES
+        flow/flow_gpu_main.cpp
+        ${FLOW_GPU_TGTS}
+        $<TARGET_OBJECTS:moduleVersion>
+      LIBRARIES
+        opmsimulators
+    )
+  endif()
+
   if(OPM_ENABLE_PYTHON)
     set_target_properties(
       flow_libblackoil
diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index cc4b74daa90..1588e6233c7 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -338,9 +338,12 @@ if(CUDA_FOUND OR hip_FOUND)
   ADD_CUDA_OR_HIP_FILE(MAIN_SOURCE_FILES opm/simulators/linalg detail/FlexibleSolverWrapper.cpp)
   ADD_CUDA_OR_HIP_FILE(MAIN_SOURCE_FILES opm/simulators/linalg FlexibleSolver_gpu_instantiate.cpp)
   ADD_CUDA_OR_HIP_FILE(MAIN_SOURCE_FILES opm/simulators/linalg PreconditionerFactory_gpu_instantiate.cpp)
+  ADD_CUDA_OR_HIP_FILE(MAIN_SOURCE_FILES opm/simulators/linalg GpuBlackoilIntensiveQuantitiesDispatcher.cu)
 
 
   # HEADERS
+  ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg GpuBlackoilIntensiveQuantitiesDispatcher.hpp)
+  ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg GpuFlowGasWaterEnergyTypeTags.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)
@@ -567,6 +570,7 @@ 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)
+  ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_gpu_ecl_thermal_law_manager.cu)
   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)
@@ -584,6 +588,7 @@ if(CUDA_FOUND OR hip_FOUND)
     ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_primary_variables_gpu.cu)
   endif()
   ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_MiniMatrix.cu)
+  ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_blackoilintensivequantities_gpu.cu)
   # Boost < 1.75 + nvcc = trouble in this test
   if(Boost_VERSION VERSION_GREATER 1.74)
     ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_MiniVector.cu)
@@ -592,6 +597,42 @@ if(CUDA_FOUND OR hip_FOUND)
   if(MPI_FOUND)
     ADD_CUDA_OR_HIP_FILE(TEST_SOURCE_FILES tests test_GpuOwnerOverlapCopy.cpp)
   endif()
+
+  # for loop providing the flag --expt-relaxed-constexpr to fix some cuda issues with constexpr
+  if(NOT CONVERT_CUDA_TO_HIP)
+    set(CU_FILES_NEEDING_RELAXED_CONSTEXPR
+      tests/gpuistl/test_gpu_ad.cu
+      tests/gpuistl/test_gpu_linear_two_phase_material.cu
+      tests/gpuistl/test_gpuPvt.cu
+      tests/gpuistl/test_gpuBlackOilFluidSystem.cu
+      tests/gpuistl/test_GpuSparseMatrix.cu
+      tests/gpuistl/test_GpuSparseTable.cu
+      tests/gpuistl/test_blackoilfluidstategpu.cu
+      flow/flow_gpu_gaswater_thermal.cu
+    )
+
+    foreach(file ${CU_FILES_NEEDING_RELAXED_CONSTEXPR})
+      set_source_files_properties(${file}
+        PROPERTIES
+        COMPILE_FLAGS
+          "--expt-relaxed-constexpr"
+      )
+    endforeach()
+
+    set(CU_FILES_NEEDING_FPERMISSIVE
+      tests/gpuistl/test_primary_variables_gpu.cu
+    )
+
+    foreach(file ${CU_FILES_NEEDING_FPERMISSIVE})
+      # Certain structures in OPM requires the -fpermissive flag to compile with nvcc,
+      # this enables this for the specific files
+      set_source_files_properties(${file}
+        PROPERTIES
+        COMPILE_FLAGS
+          "-fpermissive --expt-relaxed-constexpr -Xcompiler=-fpermissive"
+      )
+    endforeach()
+  endif()
 endif()
 
 if(USE_GPU_BRIDGE)
@@ -788,6 +829,9 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/models/discretization/common/linearizationtype.hh
   opm/models/discretization/common/restrictprolong.hh
   opm/models/discretization/common/tpfalinearizer.hh
+  opm/models/discretization/common/tpfalinearizergpukernels.hh
+  opm/models/discretization/common/tpfalinearizergpuparams.hh
+  opm/models/discretization/common/tpfalinearizerstructs.hh
   opm/models/discretization/common/tpsalinearizer.hpp
   opm/models/discretization/ecfv/ecfvbaseoutputmodule.hh
   opm/models/discretization/ecfv/ecfvdiscretization.hh
@@ -963,6 +1007,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/FacePropertiesTPSA_impl.hpp
   opm/simulators/flow/FemCpGridCompat.hpp
   opm/simulators/flow/FIBlackoilModel.hpp
+  opm/simulators/flow/SimpleFIBlackOilModel.hpp
   opm/simulators/flow/FIPContainer.hpp
   opm/simulators/flow/FlowBaseProblemProperties.hpp
   opm/simulators/flow/FlowBaseVanguard.hpp
@@ -1019,6 +1064,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/SimulatorSerializer.hpp
   opm/simulators/flow/SolutionContainers.hpp
   opm/simulators/flow/SubDomain.hpp
+  opm/simulators/flow/ThermalGasWaterFlowProblem.hpp
   opm/simulators/flow/TTagFlowProblemTPFA.hpp
   opm/simulators/flow/TTagFlowProblemTPSA.hpp
   opm/simulators/flow/TTagFlowProblemGasWater.hpp
diff --git a/flow/flow_gpu.hpp b/flow/flow_gpu.hpp
new file mode 100644
index 00000000000..89f6b7442d1
--- /dev/null
+++ b/flow/flow_gpu.hpp
@@ -0,0 +1,79 @@
+/*
+  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 FLOW_GPU_HPP
+#define FLOW_GPU_HPP
+
+#include <opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp>
+#include <opm/simulators/flow/SimpleFIBlackOilModel.hpp>
+/*
+    This file extracts typetag declarations that must be present in both the .CU and .HIP
+   executables for Flow to avoid double maintenance.
+*/
+
+namespace Opm
+{
+
+namespace Properties
+{
+    namespace TTag
+    {
+        // FlowGasWaterEnergyProblemGPU is declared in FlowGasWaterEnergyTypeTag.hpp
+        // (InheritsFrom = FlowGasWaterEnergyProblem).  The template below maps it
+        // to FlowGasWaterEnergyProblemGPUTrue<Storage> for GPU-storage variants.
+
+        template <template <class> class Storage>
+        struct FlowGasWaterEnergyProblemGPUTrue {
+            using InheritsFrom = std::tuple<FlowGasWaterEnergyProblemGPU>;
+        };
+
+        template <template <class> class Storage>
+        struct to_gpu_type<FlowGasWaterEnergyProblemGPU, Storage> {
+            using type = FlowGasWaterEnergyProblemGPUTrue<Storage>;
+        };
+    } // namespace TTag
+
+    template <class TypeTag>
+    struct RunAssemblyOnGpu<TypeTag, TTag::FlowGasWaterEnergyProblemGPU> {
+        static constexpr bool value = true;
+    };
+
+    template <class TypeTag>
+    struct Scalar<TypeTag, TTag::FlowGasWaterEnergyProblemGPU> {
+        using type = double;
+    };
+
+    template <class TypeTag>
+    struct GpuFIBlackOilModel<TypeTag, TTag::FlowGasWaterEnergyProblemGPU> {
+        using type = SimpleFIBlackOilModel<TypeTag>;
+    };
+
+    template <class TypeTag, template <class> class Storage>
+    struct FluidSystem<TypeTag, TTag::FlowGasWaterEnergyProblemGPUTrue<Storage>> {
+        using type = Opm::
+            BlackOilFluidSystemNonStatic<double, Opm::BlackOilDefaultFluidSystemIndices, Storage>;
+    };
+} // namespace Properties
+
+//! \brief Main function used in flow binary.
+int flowGasWaterEnergyMainGPU(int argc, char** argv, bool outputCout, bool outputFiles);
+
+//! \brief Main function used in flow_gaswater binary.
+int flowGasWaterEnergyMainGPUStandalone(int argc, char** argv);
+
+} // namespace Opm
+
+#endif // FLOW_GPU_HPP
diff --git a/flow/flow_gpu_gaswater_thermal.cu b/flow/flow_gpu_gaswater_thermal.cu
new file mode 100644
index 00000000000..05868b6a5f0
--- /dev/null
+++ b/flow/flow_gpu_gaswater_thermal.cu
@@ -0,0 +1,70 @@
+/*
+  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/>.
+*/
+#include "config.h"
+
+#undef HAVE_DUNE_ALUGRID // Disable ALUGRID if it is available as nvcc struggles with it
+#undef HAVE_DUNE_FEM // Disable DUNE-FEM if it is available as nvcc struggles with it
+
+// For now flow_gpu is developed to support SPE11 simulations, that is 2-phase gas-water flow with
+// thermal effects The goal is to support both property-evaluation and matrix assembly on the GPU,
+// in addition to the linear solver which already works on the GPU. The CPU would still have to
+// manage well contributions, although there are no wells in SPE11, only source-terms.
+
+#include <flow/flow_gpu.hpp>
+
+#include <opm/common/utility/gpuDecorators.hpp>
+
+#include <opm/material/common/ResetLocale.hpp>
+#include <opm/material/thermal/EnergyModuleType.hpp>
+#include <opm/models/blackoil/blackoiltwophaseindices.hh>
+
+#include <opm/grid/CpGrid.hpp>
+#include <opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp>
+#include <opm/simulators/flow/Main.hpp>
+#include <opm/simulators/flow/SimpleFIBlackOilModel.hpp>
+#include <opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp>
+
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
+#include <opm/models/discretization/common/tpfalinearizer.hh>
+
+namespace Opm
+{
+
+int
+flowGasWaterEnergyMainGPU(int argc, char** argv, bool outputCout, bool outputFiles)
+{
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    resetLocale();
+
+    FlowMain<Properties::TTag::FlowGasWaterEnergyProblemGPU> mainfunc {
+        argc, argv, outputCout, outputFiles};
+    return mainfunc.execute();
+}
+
+int
+flowGasWaterEnergyMainGPUStandalone(int argc, char** argv)
+{
+    using TypeTag = Properties::TTag::FlowGasWaterEnergyProblemGPU;
+    auto mainObject = std::make_unique<::Opm::Main>(argc, argv);
+    auto ret = mainObject->runStatic<TypeTag>();
+    // Destruct mainObject as the destructor calls MPI_Finalize!
+    mainObject.reset();
+    return ret;
+}
+
+} // namespace Opm
diff --git a/flow/flow_gpu_gaswater_thermal.hip b/flow/flow_gpu_gaswater_thermal.hip
new file mode 100644
index 00000000000..2c676de3556
--- /dev/null
+++ b/flow/flow_gpu_gaswater_thermal.hip
@@ -0,0 +1,67 @@
+/*
+  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/>.
+*/
+#include "config.h"
+
+// For now flow_gpu is developed to support SPE11 simulations, that is 2-phase gas-water flow with
+// thermal effects The goal is to support both property-evaluation and matrix assembly on the GPU,
+// in addition to the linear solver which already works on the GPU. The CPU would still have to
+// manage well contributions, although there are no wells in SPE11, only source-terms.
+
+#include <flow/flow_gpu.hpp>
+
+#include <opm/common/utility/gpuDecorators.hpp>
+
+#include <opm/material/common/ResetLocale.hpp>
+#include <opm/material/thermal/EnergyModuleType.hpp>
+
+#include <opm/grid/CpGrid.hpp>
+#include <opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp>
+#include <opm/simulators/flow/Main.hpp>
+#include <opm/simulators/flow/SimpleFIBlackOilModel.hpp>
+#include <opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp>
+
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
+#include <opm/models/blackoil/blackoiltwophaseindices.hh>
+#include <opm/models/discretization/common/tpfalinearizer.hh>
+
+namespace Opm
+{
+
+int
+flowGasWaterEnergyMainGPU(int argc, char** argv, bool outputCout, bool outputFiles)
+{
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    resetLocale();
+
+    FlowMain<Properties::TTag::FlowGasWaterEnergyProblemGPU> mainfunc {
+        argc, argv, outputCout, outputFiles};
+    return mainfunc.execute();
+}
+
+int
+flowGasWaterEnergyMainGPUStandalone(int argc, char** argv)
+{
+    using TypeTag = Properties::TTag::FlowGasWaterEnergyProblemGPU;
+    auto mainObject = std::make_unique<::Opm::Main>(argc, argv);
+    auto ret = mainObject->runStatic<TypeTag>();
+    // Destruct mainObject as the destructor calls MPI_Finalize!
+    mainObject.reset();
+    return ret;
+}
+
+} // namespace Opm
diff --git a/flow/flow_gpu_gaswater_thermal_main.cpp b/flow/flow_gpu_gaswater_thermal_main.cpp
new file mode 100644
index 00000000000..4664876220e
--- /dev/null
+++ b/flow/flow_gpu_gaswater_thermal_main.cpp
@@ -0,0 +1,25 @@
+/*
+  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/>.
+*/
+#include "config.h"
+#include <flow/flow_gpu.hpp>
+
+
+int
+main(int argc, char** argv)
+{
+    return Opm::flowGasWaterEnergyMainGPUStandalone(argc, argv);
+}
diff --git a/flow/flow_gpu_main.cpp b/flow/flow_gpu_main.cpp
new file mode 100644
index 00000000000..4664876220e
--- /dev/null
+++ b/flow/flow_gpu_main.cpp
@@ -0,0 +1,25 @@
+/*
+  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/>.
+*/
+#include "config.h"
+#include <flow/flow_gpu.hpp>
+
+
+int
+main(int argc, char** argv)
+{
+    return Opm::flowGasWaterEnergyMainGPUStandalone(argc, argv);
+}
diff --git a/opm/models/blackoil/blackoilconvectivemixingmodule.hh b/opm/models/blackoil/blackoilconvectivemixingmodule.hh
index 7a8b87ff29c..8d3442627e9 100644
--- a/opm/models/blackoil/blackoilconvectivemixingmodule.hh
+++ b/opm/models/blackoil/blackoilconvectivemixingmodule.hh
@@ -35,6 +35,7 @@
 
 #include <opm/material/common/MathToolbox.hpp>
 #include <opm/material/common/Valgrind.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 
 #include <opm/models/common/multiphasebaseproperties.hh>
 #include <opm/models/blackoil/blackoilenergymodules.hh>
@@ -164,10 +165,8 @@ class BlackOilConvectiveMixingModule<TypeTag, /*enableConvectiveMixing=*/true>
     using Toolbox = MathToolbox<Evaluation>;
     using ConvectiveMixingModuleParamT = ConvectiveMixingModuleParam<Scalar>;
 
-    enum { conti0EqIdx = Indices::conti0EqIdx };
-    enum { dimWorld = GridView::dimensionworld };
-    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
-    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
+    static constexpr int conti0EqIdx = Indices::conti0EqIdx;
+    static constexpr int dimWorld = GridView::dimensionworld;
     static constexpr bool enableFullyImplicitThermal = (getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal);
     static constexpr unsigned contiEnergyEqIdx = Indices::contiEnergyEqIdx;
 
@@ -200,6 +199,10 @@ public:
                                                  const IntensiveQuantities& intQuantsEx,
                                                  const unsigned phaseIdx,
                                                  const CMMParam& info) {
+        // Local constexpr aliases avoid nvcc's failure to resolve class-scope
+        // static constexpr members in device code.
+        constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
+        constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
 
         if (info.active_.empty()) {
             return;
@@ -214,9 +217,9 @@ public:
         }
 
         const auto& liquidPhaseIdx =
-            fsys.phaseIsActive(fsys.waterPhaseIdx)
-               ? fsys.waterPhaseIdx
-               : fsys.oilPhaseIdx;
+            fsys.phaseIsActive(waterPhaseIdx)
+               ? waterPhaseIdx
+               : oilPhaseIdx;
 
         // Compute avg density based on pure water
         const auto& t_in = intQuantsIn.fluidState().temperature(liquidPhaseIdx);
@@ -309,6 +312,11 @@ public:
                                                         const Scalar faceArea,
                                                         const CMMParam& info)
     {
+        // Local constexpr aliases avoid nvcc's failure to resolve class-scope
+        // static constexpr members in device code.
+        constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
+        constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
+
         const FluidSystem& fsys = intQuantsIn.getFluidSystem();
 
         if (info.active_.empty()) {
@@ -320,9 +328,9 @@ public:
         }
 
         const auto& liquidPhaseIdx =
-            fsys.phaseIsActive(fsys.waterPhaseIdx)
-                ? fsys.waterPhaseIdx
-                : fsys.oilPhaseIdx;
+            fsys.phaseIsActive(waterPhaseIdx)
+                ? waterPhaseIdx
+                : oilPhaseIdx;
 
         // interiour
         const auto& t_in = intQuantsIn.fluidState().temperature(liquidPhaseIdx);
@@ -331,14 +339,14 @@ public:
         const auto& salt_in = intQuantsIn.fluidState().saltSaturation();
 
         const auto bLiquidSatIn =
-            fsys.phaseIsActive(fsys.waterPhaseIdx)
+            fsys.phaseIsActive(waterPhaseIdx)
                 ? fsys.waterPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(),
                                                                        t_in, p_in, rssat_in, salt_in)
                 : fsys.oilPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(),
                                                                      t_in, p_in, rssat_in);
 
         const auto& densityLiquidIn =
-            fsys.phaseIsActive(fsys.waterPhaseIdx)
+            fsys.phaseIsActive(waterPhaseIdx)
                 ? fsys.waterPvt().waterReferenceDensity(intQuantsIn.pvtRegionIndex())
                 : fsys.oilPvt().oilReferenceDensity(intQuantsIn.pvtRegionIndex());
 
@@ -354,14 +362,14 @@ public:
         const auto rssat_ex = Opm::getValue(intQuantsEx.saturatedDissolutionFactor());
         const auto salt_ex = Opm::getValue(intQuantsEx.fluidState().saltSaturation());
         const auto bLiquidSatEx =
-            fsys.phaseIsActive(fsys.waterPhaseIdx)
+            fsys.phaseIsActive(waterPhaseIdx)
                 ? fsys.waterPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(),
                                                                        t_ex, p_ex, rssat_ex, salt_ex)
                 : fsys.oilPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(),
                                                                      t_ex, p_ex, rssat_ex);
 
         const auto& densityLiquidEx =
-            fsys.phaseIsActive(fsys.waterPhaseIdx)
+            fsys.phaseIsActive(waterPhaseIdx)
                 ? fsys.waterPvt().waterReferenceDensity(intQuantsEx.pvtRegionIndex())
                 : fsys.oilPvt().oilReferenceDensity(intQuantsEx.pvtRegionIndex());
 
@@ -389,7 +397,7 @@ public:
             const auto& rssat_up = (upIdx == interiorDofIdx) ? rssat_in : rssat_ex;
             unsigned globalUpIndex = (upIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
             const auto& Rsup =
-                fsys.phaseIsActive(fsys.waterPhaseIdx)
+                fsys.phaseIsActive(waterPhaseIdx)
                     ? up.fluidState().Rsw()
                     : up.fluidState().Rs();
 
diff --git a/opm/models/blackoil/blackoildiffusionmodule.hh b/opm/models/blackoil/blackoildiffusionmodule.hh
index 51db5942f01..28eacd4281d 100644
--- a/opm/models/blackoil/blackoildiffusionmodule.hh
+++ b/opm/models/blackoil/blackoildiffusionmodule.hh
@@ -362,6 +362,12 @@ class BlackOilDiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/false>
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
 
 public:
+    BlackOilDiffusionIntensiveQuantities() = default;
+
+    template <class OtherTypeTag>
+    BlackOilDiffusionIntensiveQuantities(const BlackOilDiffusionIntensiveQuantities<OtherTypeTag, false>&)
+    {}
+
     /*!
      * \brief Returns the tortuousity of the sub-domain of a fluid
      *        phase in the porous medium.
@@ -441,6 +447,12 @@ public:
         , diffusionCoefficient_(diffusionCoefficient)
     {}
 
+    template <class OtherTypeTag>
+    BlackOilDiffusionIntensiveQuantities(const BlackOilDiffusionIntensiveQuantities<OtherTypeTag, true>& other)
+        : tortuosity_(other.tortuosity())
+        , diffusionCoefficient_(other.diffusionCoefficients())
+    {}
+
     BlackOilDiffusionIntensiveQuantities&
     operator=(const BlackOilDiffusionIntensiveQuantities& rhs)
     {
@@ -462,6 +474,12 @@ public:
     OPM_HOST_DEVICE Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
     { return diffusionCoefficient_[phaseIdx][compIdx]; }
 
+    /*!
+     * \brief Returns all the diffusion coefficients
+     */
+    OPM_HOST_DEVICE const std::array<std::array<Evaluation, numComponents>, numPhases>& diffusionCoefficients() const
+    { return diffusionCoefficient_; }
+
     /*!
      * \brief Returns the tortuousity of the sub-domain of a fluid
      *        phase in the porous medium.
@@ -469,6 +487,12 @@ public:
     OPM_HOST_DEVICE Evaluation tortuosity(unsigned phaseIdx) const
     { return tortuosity_[phaseIdx]; }
 
+    /*!
+     * \brief Returns all the tortuosities
+     */
+    OPM_HOST_DEVICE const std::array<Evaluation, numPhases>& tortuosities() const
+    { return tortuosity_; }
+
     /*!
      * \brief Returns the effective molecular diffusion coefficient of
      *        the porous medium for a component in a phase.
@@ -695,9 +719,9 @@ protected:
     }
 
 public:
-    static void update(EvaluationArray& effectiveDiffusionCoefficient,
-                       const IntensiveQuantities& intQuantsInside,
-                       const IntensiveQuantities& intQuantsOutside)
+    OPM_HOST_DEVICE static void update(EvaluationArray& effectiveDiffusionCoefficient,
+                                       const IntensiveQuantities& intQuantsInside,
+                                       const IntensiveQuantities& intQuantsOutside)
     {
         const FluidSystem& fsys = intQuantsInside.getFluidSystem();
 
diff --git a/opm/models/blackoil/blackoilenergymodules.hh b/opm/models/blackoil/blackoilenergymodules.hh
index e2cce9cd2a5..0d58f850d34 100644
--- a/opm/models/blackoil/blackoilenergymodules.hh
+++ b/opm/models/blackoil/blackoilenergymodules.hh
@@ -380,6 +380,14 @@ public:
     {
     }
 
+    template <class OtherTypeTag>
+    BlackOilEnergyIntensiveQuantities(
+        const BlackOilEnergyIntensiveQuantities<OtherTypeTag, EnergyModules::FullyImplicitThermal>& other)
+        : rockInternalEnergy_(other.rockInternalEnergy())
+        , totalThermalConductivity_(other.totalThermalConductivity())
+        , rockFraction_(other.rockFraction())
+    {}
+
     BlackOilEnergyIntensiveQuantities() = default;
 
     /*!
diff --git a/opm/models/blackoil/blackoilextbomodules.hh b/opm/models/blackoil/blackoilextbomodules.hh
index 0a0cce88532..804ab00cadb 100644
--- a/opm/models/blackoil/blackoilextbomodules.hh
+++ b/opm/models/blackoil/blackoilextbomodules.hh
@@ -187,7 +187,7 @@ public:
                 storage[contiZfracEqIdx - 1] += regWghtFactor*Toolbox::template decay<LhsEval>(intQuants.zFraction());
             }
             else {
-                throw std::runtime_error("Only component conservation in terms of surface volumes is implemented. ");
+                OPM_THROW(std::runtime_error, "Only component conservation in terms of surface volumes is implemented. ");
             }
         }
     }
diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh
index c249bdd3b9f..20a840aa85b 100644
--- a/opm/models/blackoil/blackoilintensivequantities.hh
+++ b/opm/models/blackoil/blackoilintensivequantities.hh
@@ -175,9 +175,8 @@ public:
             fluidState_.setRsw(0.0);
         }
     }
-    BlackOilIntensiveQuantities(const BlackOilIntensiveQuantities& other) = default;
 
-    BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) = default;
+    OPM_HOST_DEVICE BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) = default;
 
     template<class OtherTypeTag>
     friend class BlackOilIntensiveQuantities;
@@ -186,26 +185,28 @@ public:
     template<class OtherTypeTag>
     explicit BlackOilIntensiveQuantities(
         const BlackOilIntensiveQuantities<OtherTypeTag>& other, const FluidSystem& fsystem)
-        : fluidState_(other.fluidState_.withOtherFluidSystem(fsystem))
-        , BlackOilEnergyIntensiveQuantities<TypeTag, energyModuleType>(
-            other.rockInternalEnergy_, other.totalThermalConductivity_, other.rockFraction_)
-        , BlackOilDiffusionIntensiveQuantities<TypeTag, enableDiffusion>(
-            other.tortuosities(), other.diffusionCoefficients())
+        : fluidState_(other.fluidState_.template withOtherFluidSystem<FluidSystem>(fsystem))
+        , BlackOilEnergyIntensiveQuantities<TypeTag, energyModuleType>(other)
+        , BlackOilDiffusionIntensiveQuantities<TypeTag, enableDiffusion>(other)
         , referencePorosity_(other.referencePorosity_)
         , porosity_(other.porosity_)
         , rockCompTransMultiplier_(other.rockCompTransMultiplier_)
         , mobility_(other.mobility_)
         , dirMob_(/*NOT YET SUPPORTED ON GPU*/)
     {
-        static_assert(!enableSolvent);
-        static_assert(!enableExtbo);
-        static_assert(!enablePolymer);
-        static_assert(!enableFoam);
-        static_assert(!enableMICP);
-        static_assert(!enableBrine);
-        static_assert(!enableDispersion);
+        // Changing diffusion or dispersion support across a TypeTag conversion would
+        // silently produce incorrect results (the missing quantities would be left at
+        // their default-initialised values). Catch such mismatches at compile time.
+        static_assert(getPropValue<OtherTypeTag, Properties::EnableDiffusion>() == enableDiffusion,
+                      "Cannot convert BlackOilIntensiveQuantities between TypeTags with "
+                      "different EnableDiffusion settings.");
+        static_assert(getPropValue<OtherTypeTag, Properties::EnableDispersion>() == enableDispersion,
+                      "Cannot convert BlackOilIntensiveQuantities between TypeTags with "
+                      "different EnableDispersion settings.");
     }
 
+    BlackOilIntensiveQuantities(const BlackOilIntensiveQuantities& other) = default;
+
     /**
      * \brief Create a copy of this intensive quantities object
      *
@@ -221,6 +222,45 @@ public:
         return newIntQuants;
     }
 
+    /*!
+     * \brief Field-by-field overlay of the BlackOil intensive-quantity values from
+     *        another \c BlackOilIntensiveQuantities instantiation onto this one.
+     *
+     * Used by the experimental GPU intensive-quantities dispatcher to write the
+     * GPU-computed result onto a CPU-side \c IntensiveQuantities even when the
+     * two TypeTags are not value-compatible (e.g. when the CPU TypeTag enables
+     * dispersion or other modules that the GPU TypeTag does not). Only stored
+     * fields that the dispatcher actually computes are copied; the
+     * \c mobility_ field is left untouched (the GPU relperm path is currently
+     * known to return zero, see \c GpuBlackoilIntensiveQuantitiesDispatcher).
+     */
+    template<class OtherTypeTag>
+    OPM_HOST_DEVICE void overlayBlackOilFieldsFrom(
+        const BlackOilIntensiveQuantities<OtherTypeTag>& other)
+    {
+        // Full fluid-state copy (handles every stored field, including
+        // Rs/Rv/Rsw/Rvw, salt concentration, solvent fields, ...). This
+        // is intentionally chosen over a hand-picked field list so that
+        // the GPU dispatcher can be used as the *only* fill of the IQ
+        // cache (no CPU pre-pass required).
+        fluidState_.assign(other.fluidState_);
+
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                continue;
+            }
+            mobility_[phaseIdx] = other.mobility_[phaseIdx];
+        }
+        if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal) {
+            this->rockInternalEnergy_ = other.rockInternalEnergy_;
+            this->totalThermalConductivity_ = other.totalThermalConductivity_;
+            this->rockFraction_ = other.rockFraction_;
+        }
+        porosity_ = other.porosity_;
+        referencePorosity_ = other.referencePorosity_;
+        rockCompTransMultiplier_ = other.rockCompTransMultiplier_;
+    }
+
     OPM_HOST_DEVICE void updateTempSalt(const Problem& problem,
                                         const PrimaryVariables& priVars,
                                         const unsigned globalSpaceIdx,
@@ -430,11 +470,15 @@ public:
             }
             else {
                 if (getFluidSystem().enableDissolvedGas()) { // Add So > 0? i.e. if only water set rs = 0)
-                    const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
-                        getFluidSystem().saturatedDissolutionFactor(fluidState_,
-                                                                oilPhaseIdx,
-                                                                pvtRegionIdx,
-                                                                SoMax);
+                    Evaluation RsSat;
+                    if constexpr (enableExtbo) {
+                        RsSat = asImp_().rs();
+                    } else {
+                        RsSat = getFluidSystem().saturatedDissolutionFactor(fluidState_,
+                                                                            oilPhaseIdx,
+                                                                            pvtRegionIdx,
+                                                                            SoMax);
+                    }
                     fluidState_.setRs(min(RsMax, RsSat));
                 }
                 else {
@@ -448,11 +492,15 @@ public:
             }
             else {
                 if (getFluidSystem().enableVaporizedOil() ) { // Add Sg > 0? i.e. if only water set rv = 0)
-                    const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
-                        getFluidSystem().saturatedDissolutionFactor(fluidState_,
-                                                                gasPhaseIdx,
-                                                                pvtRegionIdx,
-                                                                SoMax);
+                    Evaluation RvSat;
+                    if constexpr (enableExtbo) {
+                        RvSat = asImp_().rv();
+                    } else {
+                        RvSat = getFluidSystem().saturatedDissolutionFactor(fluidState_,
+                                                                            gasPhaseIdx,
+                                                                            pvtRegionIdx,
+                                                                            SoMax);
+                    }
                     fluidState_.setRv(min(RvMax, RvSat));
                 }
                 else {
@@ -514,11 +562,16 @@ public:
             const auto [b, mu] = getFluidSystem().inverseFormationVolumeFactorAndViscosity(fluidState_, phaseIdx, pvtRegionIdx);
             fluidState_.setInvB(phaseIdx, b);
             for (int i = 0; i < nmobilities; ++i) {
-                if (enableExtbo && phaseIdx == oilPhaseIdx) {
-                    (*mobilities[i])[phaseIdx] /= asImp_().oilViscosity();
-                }
-                else if (enableExtbo && phaseIdx == gasPhaseIdx) {
-                    (*mobilities[i])[phaseIdx] /= asImp_().gasViscosity();
+                if constexpr (enableExtbo) {
+                    if (phaseIdx == oilPhaseIdx) {
+                        (*mobilities[i])[phaseIdx] /= asImp_().oilViscosity();
+                    }
+                    else if (phaseIdx == gasPhaseIdx) {
+                        (*mobilities[i])[phaseIdx] /= asImp_().gasViscosity();
+                    }
+                    else {
+                        (*mobilities[i])[phaseIdx] /= mu;
+                    }
                 }
                 else {
                     (*mobilities[i])[phaseIdx] /= mu;
@@ -845,7 +898,7 @@ public:
                     case Dir::ZPlus:
                         return dirMob_->getArray(2)[phaseIdx];
                     default:
-                        throw std::runtime_error("Unexpected face direction");
+                        OPM_THROW(std::runtime_error, "Unexpected face direction");
                 }
             }
             else{
@@ -906,7 +959,7 @@ public:
             return BrineIntQua::permFactor();
         }
         else {
-            throw std::logic_error("permFactor() called but salt precipitation or bioeffects are disabled");
+            OPM_THROW(std::logic_error, "permFactor() called but salt precipitation and bioeffects are disabled");
         }
     }
 
diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh
index 498678950b7..5955690d184 100644
--- a/opm/models/blackoil/blackoillocalresidualtpfa.hh
+++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh
@@ -47,9 +47,11 @@
 #include <opm/models/blackoil/blackoilpolymermodules.hh>
 #include <opm/models/blackoil/blackoilproperties.hh>
 #include <opm/models/blackoil/blackoilsolventmodules.hh>
+#include <opm/models/blackoil/blackoilmoduleparams.hh>
 
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/utility/gpuDecorators.hpp>
+#include <opm/common/utility/gpuistl_if_available.hpp>
 
 #include <array>
 #include <cassert>
@@ -573,49 +575,6 @@ public:
         }
     }
 
-    template <class BoundaryConditionData, class RateVectorLocal, class LocalProblem>
-    OPM_HOST_DEVICE static void computeBoundaryFlux(RateVectorLocal& bdyFlux,
-                                                    const LocalProblem& problem,
-                                                    const BoundaryConditionData& bdyInfo,
-                                                    const IntensiveQuantities& insideIntQuants,
-                                                    unsigned globalSpaceIdx)
-    {
-#if OPM_IS_INSIDE_HOST_FUNCTION
-        switch (bdyInfo.type) {
-        case BCType::NONE:
-            bdyFlux = 0.0;
-            break;
-        case BCType::RATE:
-            computeBoundaryFluxRate(bdyFlux, bdyInfo);
-            break;
-        case BCType::FREE:
-        case BCType::DIRICHLET:
-            computeBoundaryFluxFree(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
-            break;
-        case BCType::THERMAL:
-            computeBoundaryThermal(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
-            break;
-        default:
-            throw std::logic_error("Unknown boundary condition type "
-                                   + std::to_string(static_cast<int>(bdyInfo.type))
-                                   + " in computeBoundaryFlux().");
-        }
-#else // TODO: support all boundary conditions on GPU as well to unify this code
-        switch (bdyInfo.type) {
-        case BCType::NONE:
-            bdyFlux = 0.0;
-            break;
-        case BCType::THERMAL:
-            computeBoundaryThermal(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
-            break;
-        default:
-            OPM_THROW(std::logic_error,
-                      "Boundary condition type " + std::to_string(static_cast<int>(bdyInfo.type))
-                          + " is not supported for GPU fluid systems in computeBoundaryFlux().");
-        }
-#endif
-    }
-
     template <class BoundaryConditionData>
     static void computeBoundaryFluxRate(RateVector& bdyFlux, const BoundaryConditionData& bdyInfo)
     {
@@ -750,17 +709,12 @@ public:
         if constexpr (enableFullyImplicitThermal) {
             Evaluation heatFlux;
             // avoid overload of functions with same numeber of elements in eclproblem
-
             Scalar alpha;
-            if constexpr (runAssemblyOnGpu) {
-                // This path is currently only intended for the SimplifiedBlackoilModel for GPUs
-                // which currently does not aim to reproduce the full problem object on the GPU.
+            if constexpr (!std::is_empty_v<GetPropType<TypeTag, Properties::FluidSystem>>) {
                 alpha = problem.getAlpha(globalSpaceIdx, bdyInfo.boundaryFaceIndex);
             } else {
-                alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(
-                    globalSpaceIdx, bdyInfo.boundaryFaceIndex);
+                alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(globalSpaceIdx, bdyInfo.boundaryFaceIndex);
             }
-
             unsigned inIdx = 0; // dummy
             // always calculated with derivatives of this cell
             EnergyModule::ExtensiveQuantities::updateEnergyBoundary(heatFlux,
@@ -780,6 +734,49 @@ public:
 #endif
     }
 
+    template <class BoundaryConditionData, class RateVectorLocal, class LocalProblem>
+    OPM_HOST_DEVICE static void computeBoundaryFlux(RateVectorLocal& bdyFlux,
+                                                    const LocalProblem& problem,
+                                                    const BoundaryConditionData& bdyInfo,
+                                                    const IntensiveQuantities& insideIntQuants,
+                                                    unsigned globalSpaceIdx)
+    {
+        if constexpr (std::is_empty_v<GetPropType<TypeTag, Properties::FluidSystem>>) {
+            switch (bdyInfo.type) {
+            case BCType::NONE:
+                bdyFlux = 0.0;
+                break;
+            case BCType::RATE:
+                computeBoundaryFluxRate(bdyFlux, bdyInfo);
+                break;
+            case BCType::FREE:
+            case BCType::DIRICHLET:
+                computeBoundaryFluxFree(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
+                break;
+            case BCType::THERMAL:
+                computeBoundaryThermal(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
+                break;
+            default:
+                OPM_THROW(std::logic_error, "Unknown boundary condition type "
+                              + std::to_string(static_cast<int>(bdyInfo.type))
+                              + " in computeBoundaryFlux().");
+            }
+        } else { // Non-static fluid system used in GPU assembly
+            switch (bdyInfo.type) {
+            case BCType::NONE:
+                bdyFlux = 0.0;
+                break;
+            case BCType::THERMAL:
+                computeBoundaryThermal(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
+                break;
+            default:
+                OPM_THROW(std::logic_error,
+                          "Boundary condition type " + std::to_string(static_cast<int>(bdyInfo.type))
+                              + " is not supported for GPU fluid systems in computeBoundaryFlux().");
+            }
+        }
+    }
+
     static void computeSource(RateVector& source,
                               const Problem& problem,
                               const IntensiveQuantities& insideIntQuants,
@@ -935,18 +932,15 @@ public:
     }
 
     /*!
-     * \brief Helper function to convert the mass-related parts of a Dune::FieldVector
-     *        that stores conservation quantities in terms of "surface-volume" to the
-     *        conservation quantities used by the model.
-     *
-     * Convenience overload for CPU code that uses the static FluidSystem. Delegates to
-     * the FsysType overload below, constructing a default FluidSystem instance.
+     * Wrapper of adaptMassConservationQuantities_ that does not take a fluid system instance.
      */
-    template <class Scalar>
-    static void adaptMassConservationQuantities_(Dune::FieldVector<Scalar, numEq>& container,
-                                                 unsigned pvtRegionIdx)
+    template <class ScalarVector>
+    OPM_HOST_DEVICE static void adaptMassConservationQuantities_(ScalarVector& container,
+                                                                 unsigned pvtRegionIdx)
     {
-        adaptMassConservationQuantities_(container, pvtRegionIdx, FluidSystem {});
+        // Delegate to the generic overload using a default-constructed static FluidSystem
+        // instance. Valid because the static FluidSystem is stateless (std::is_empty_v).
+        adaptMassConservationQuantities_(container, pvtRegionIdx, FluidSystem{});
     }
 
     /*!
diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh
index c324d4ad6b2..76e937b830d 100644
--- a/opm/models/discretization/common/fvbasediscretization.hh
+++ b/opm/models/discretization/common/fvbasediscretization.hh
@@ -660,6 +660,9 @@ public:
         return &intensiveQuantityCache_[timeIdx][globalIdx];
     }
 
+    const auto& intensiveQuantityCache() const
+    { return intensiveQuantityCache_; }
+
     /*!
      * \brief Update the intensive quantity cache for a entity on the grid at given time.
      *
@@ -1204,6 +1207,9 @@ public:
     bool isLocalDof(unsigned globalIdx) const
     { return isLocalDof_[globalIdx]; }
 
+    size_t numDof() const
+    { return asImp_().numGridDof(); }
+
     /*!
      * \brief Returns the volume \f$\mathrm{[m^3]}\f$ of the whole grid which represents
      *        the spatial domain.
diff --git a/opm/models/discretization/common/fvbaseproperties.hh b/opm/models/discretization/common/fvbaseproperties.hh
index 4f3f6dd3ad5..c54afe6e43e 100644
--- a/opm/models/discretization/common/fvbaseproperties.hh
+++ b/opm/models/discretization/common/fvbaseproperties.hh
@@ -46,6 +46,12 @@ namespace TTag {
 struct FvBaseDiscretization
 { using InheritsFrom = std::tuple<FvBaseNewtonMethod, ImplicitModel>; };
 
+template<class TypeTag, template<class> class Storage>
+struct to_gpu_type {};
+
+template<class TypeTag, template <class> class Storage>
+using to_gpu_type_t = typename to_gpu_type<TypeTag, Storage>::type;
+
 } // namespace TTag
 
 
@@ -131,6 +137,11 @@ struct PrimaryVariables { using type = UndefinedProperty; };
 //! The secondary variables within a sub-control volume
 template<class TypeTag, class MyTypeTag>
 struct IntensiveQuantities { using type = UndefinedProperty; };
+
+//! The blackoilmodel that is tailored for GPU compatibility
+template<class TypeTag, class MyTypeTag>
+struct GpuFIBlackOilModel { using type = UndefinedProperty; };
+
 //! The discretization specific part of the intensive quantities
 
 template<class TypeTag, class MyTypeTag>
diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh
index 24e2ec0f9f0..ff4ddd0bcaf 100644
--- a/opm/models/discretization/common/tpfalinearizer.hh
+++ b/opm/models/discretization/common/tpfalinearizer.hh
@@ -47,8 +47,13 @@
 #include <opm/models/discretization/common/baseauxiliarymodule.hh>
 #include <opm/models/discretization/common/fvbaseproperties.hh>
 #include <opm/models/discretization/common/linearizationtype.hh>
+#include <opm/models/discretization/common/tpfalinearizerstructs.hh>
+
 #include <opm/simulators/linalg/exportSystem.hpp>
 
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystemNonStatic.hpp>
+
 #include <cassert>
 #include <cstddef>
 #include <exception>   // current_exception, rethrow_exception
@@ -61,20 +66,30 @@
 #include <unordered_map>
 #include <vector>
 
+#if HAVE_OPENMP
+#include <omp.h>
+#endif
+
 #include <fmt/format.h>
 
 #include <opm/common/utility/gpuDecorators.hpp>
+#include <opm/common/utility/gpuistl_if_available.hpp>
+#include <opm/common/utility/pointerArithmetic.hpp>
 #if HAVE_CUDA
+#include <opm/models/discretization/common/tpfalinearizergpukernels.hh>
+#include <opm/simulators/flow/SimpleFIBlackOilModel.hpp>
+#include <opm/simulators/flow/ThermalGasWaterFlowProblem.hpp>
+#include <opm/models/discretization/common/tpfalinearizergpuparams.hh>
 #if USE_HIP
-#include <opm/simulators/linalg/gpuistl_hip/GpuBuffer.hpp>
-#include <opm/simulators/linalg/gpuistl_hip/GpuView.hpp>
+#include <opm/simulators/linalg/gpuistl_hip/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/gpuistl_hip/MiniMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl_hip/MiniVector.hpp>
+#include <opm/simulators/linalg/gpuistl_hip/detail/gpusparse_matrix_operations.hpp>
 #else
-#include <opm/simulators/linalg/gpuistl/GpuBuffer.hpp>
-#include <opm/simulators/linalg/gpuistl/GpuView.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/gpuistl/MiniMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/MiniVector.hpp>
+#include <opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp>
 #endif
 #endif
 
@@ -90,91 +105,6 @@ namespace Opm {
 template<class TypeTag>
 class EcfvDiscretization;
 
-// Moved these structs out of the class to make them visible in the GPU code.
-template<class Storage = std::vector<int>>
-struct FullDomain
-{
-    Storage cells;
-};
-
-#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
-    inline FullDomain<gpuistl::GpuBuffer<int>> copy_to_gpu(FullDomain<> CPUDomain)
-    {
-        if (CPUDomain.cells.size() == 0) {
-            OPM_THROW(std::runtime_error, "Cannot copy empty full domain to GPU.");
-        }
-        return FullDomain<gpuistl::GpuBuffer<int>>{
-            gpuistl::GpuBuffer<int>(CPUDomain.cells)
-        };
-    };
-
-    inline FullDomain<gpuistl::GpuView<int>> make_view(FullDomain<gpuistl::GpuBuffer<int>>& buffer)
-    {
-        if (buffer.cells.size() == 0) {
-            OPM_THROW(std::runtime_error, "Cannot make view of empty full domain buffer.");
-        }
-        return FullDomain<gpuistl::GpuView<int>>{
-            gpuistl::make_view(buffer.cells)
-        };
-    };
-#endif
-
-template <class ResidualNBInfoType,class BlockType>
-struct NeighborInfoStruct
-{
-    unsigned int neighbor;
-    ResidualNBInfoType res_nbinfo;
-    BlockType* matBlockAddress;
-
-    template <class OtherBlockType>
-    NeighborInfoStruct(const NeighborInfoStruct<ResidualNBInfoType,OtherBlockType>& other)
-        : neighbor(other.neighbor)
-        , res_nbinfo(other.res_nbinfo)
-        , matBlockAddress(nullptr)
-    {
-        if (other.matBlockAddress) {
-            matBlockAddress = reinterpret_cast<BlockType*>(other.matBlockAddress);
-        }
-    }
-
-    template <class PtrType>
-    NeighborInfoStruct(unsigned int n, const ResidualNBInfoType& r, PtrType ptr)
-        : neighbor(n)
-        , res_nbinfo(r)
-        , matBlockAddress(static_cast<BlockType*>(ptr))
-    {
-    }
-
-    // Add a default constructor
-    NeighborInfoStruct()
-        : neighbor(0)
-        , res_nbinfo()
-        , matBlockAddress(nullptr)
-    {
-    }
-};
-
-#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
-namespace  gpuistl {
-    template<class MiniMatrixType, class MatrixBlockType, class ResidualNBInfoType>
-    auto copy_to_gpu(const SparseTable<NeighborInfoStruct<ResidualNBInfoType, MatrixBlockType>>& cpuNeighborInfoTable)
-    {
-        // Convert the DUNE FieldMatrix/MatrixBlock to MiniMatrix types
-        using StructWithMinimatrix = NeighborInfoStruct<ResidualNBInfoType, MiniMatrixType>;
-        std::vector<StructWithMinimatrix> minimatrices(cpuNeighborInfoTable.dataSize());
-        size_t idx = 0;
-        for (auto e : cpuNeighborInfoTable.dataStorage()) {
-            minimatrices[idx++] = StructWithMinimatrix(e);
-        }
-
-        return SparseTable<StructWithMinimatrix, gpuistl::GpuBuffer>(
-            gpuistl::GpuBuffer<StructWithMinimatrix>(minimatrices),
-            gpuistl::GpuBuffer<int>(cpuNeighborInfoTable.rowStarts())
-        );
-    }
-}
-#endif
-
 /*!
  * \ingroup FiniteVolumeDiscretizations
  *
@@ -210,17 +140,29 @@ class TpfaLinearizer
 
     using Vector = GlobalEqVector;
 
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
     enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
     enum { historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>() };
     enum { dimWorld = GridView::dimensionworld };
 
-    using MatrixBlock = typename SparseMatrixAdapter::MatrixBlock;
-    using VectorBlock = Dune::FieldVector<Scalar, numEq>;
-    using ADVectorBlock = GetPropType<TypeTag, Properties::RateVector>;
+    constexpr static bool runAssemblyOnGpu = getPropValue<TypeTag, Properties::RunAssemblyOnGpu>();
 
-#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
-    using MatrixBlockGPU = gpuistl::MiniMatrix<Scalar, numEq * numEq>;
+    using MatrixBlockCPU = typename SparseMatrixAdapter::MatrixBlock;
+    using VectorBlockCPU = Dune::FieldVector<Scalar, numEq>;
+    using ADVectorBlockCPU = GetPropType<TypeTag, Properties::RateVector>;
+
+#if HAVE_CUDA
+    using MatrixBlockGPU = gpuistl::MiniMatrix<Scalar, numEq>;
     using VectorBlockGPU = gpuistl::MiniVector<Scalar, numEq>;
+    using ADVectorBlockGPU = gpuistl::MiniVector<Evaluation, numEq>;
+    using MatrixBlock = std::conditional_t<runAssemblyOnGpu, MatrixBlockGPU, MatrixBlockCPU>;
+    using VectorBlock = std::conditional_t<runAssemblyOnGpu, VectorBlockGPU, VectorBlockCPU>;
+    using ADVectorBlock = std::conditional_t<runAssemblyOnGpu, ADVectorBlockGPU, ADVectorBlockCPU>;
+#else
+    using MatrixBlock = MatrixBlockCPU;
+    using VectorBlock = VectorBlockCPU;
+    using ADVectorBlock = ADVectorBlockCPU;
 #endif
 
     static constexpr bool linearizeNonLocalElements =
@@ -484,7 +426,7 @@ public:
             const auto [type, massrateAD] = problem_().boundaryCondition(bdyInfo.cell, bdyInfo.dir);
 
             // Strip the unnecessary (and zero anyway) derivatives off massrate.
-            VectorBlock massrate(0.0);
+            VectorBlockCPU massrate(0.0);
             for (std::size_t ii = 0; ii < massrate.size(); ++ii) {
                 massrate[ii] = massrateAD[ii].value();
             }
@@ -570,7 +512,7 @@ private:
         std::vector<NeighborSet> sparsityPattern(model.numTotalDof());
         const Scalar gravity = problem_().gravity()[dimWorld - 1];
         unsigned numCells = model.numTotalDof();
-        neighborInfo_.reserve(numCells, 6 * numCells);
+        neighborInfo_.reserve(numCells, 6 * numCells); // Expect ~6 neighbors per cell
         std::vector<NeighborInfoCPU> loc_nbinfo;
         for (const auto& elem : elements(gridView_())) {
             stencil.update(elem);
@@ -633,18 +575,18 @@ private:
                         }
                         const auto [type, massrateAD] = problem_().boundaryCondition(myIdx, dir_id);
                         // Strip the unnecessary (and zero anyway) derivatives off massrate.
-                        VectorBlock massrate(0.0);
+                        VectorBlockCPU massrate(0.0);
                         for (std::size_t ii = 0; ii < massrate.size(); ++ii) {
                             massrate[ii] = massrateAD[ii].value();
                         }
                         const auto& exFluidState = problem_().boundaryFluidState(myIdx, dir_id);
-                        BoundaryConditionData bcdata{type,
-                                                     massrate,
-                                                     exFluidState.pvtRegionIndex(),
-                                                     bfIndex,
-                                                     bf.area(),
-                                                     bf.integrationPos()[dimWorld - 1],
-                                                     exFluidState};
+                        BoundaryConditionDataCPU bcdata {type,
+                                                         massrate,
+                                                         exFluidState.pvtRegionIndex(),
+                                                         bfIndex,
+                                                         bf.area(),
+                                                         bf.integrationPos()[dimWorld - 1],
+                                                         exFluidState};
                         boundaryInfo_.push_back({myIdx, dir_id, bfIndex, bcdata});
                     }
                 }
@@ -671,6 +613,13 @@ private:
             }
         }
 
+#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+        gpuJacobian_.reset(new gpuistl::GpuSparseMatrixWrapper<Scalar>(
+            gpuistl::GpuSparseMatrixWrapper<Scalar>::fromMatrix(jacobian_->istlMatrix())));
+        gpuBufferDiagMatAddress_.reset(new gpuistl::GpuBuffer<MatrixBlockGPU*>(
+            gpuistl::detail::getDiagPtrsTyped<MatrixBlockGPU>(*gpuJacobian_)));
+#endif
+
         // Create dummy full domain.
         fullDomain_.cells.resize(numCells);
         std::iota(fullDomain_.cells.begin(), fullDomain_.cells.end(), 0);
@@ -682,6 +631,10 @@ private:
         residual_ = 0.0;
         // zero all matrix entries
         jacobian_->clear();
+
+#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+        gpuJacobian_->setToZero();
+#endif
     }
 
     // Initialize the flows, flores, and velocity sparse tables
@@ -810,7 +763,9 @@ private:
     }
 
 public:
-    void setResAndJacobi(VectorBlock& res, MatrixBlock& bMat, const ADVectorBlock& resid) const
+    template <class VectorBlockType, class MatrixBlockType, class ADVectorBlockType>
+    OPM_HOST_DEVICE static void
+    setResAndJacobi(VectorBlockType& res, MatrixBlockType& bMat, const ADVectorBlockType& resid)
     {
         for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
             res[eqIdx] = resid[eqIdx].value();
@@ -891,7 +846,7 @@ public:
                 continue;
             }
 
-            ADVectorBlock adres(0.0);
+            ADVectorBlockCPU adres(0.0);
             const unsigned globI = bdyInfo.cell;
             const auto& nbInfos = neighborInfo_[globI];
             const IntensiveQuantities& insideIntQuants = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
@@ -911,6 +866,8 @@ private:
     template <class SubDomainType>
     void linearize_(const SubDomainType& domain)
     {
+        constexpr bool run_assembly_on_gpu = getPropValue<TypeTag, Properties::RunAssemblyOnGpu>();
+
         // This check should be removed once this is addressed by
         // for example storing the previous timesteps' values for
         // rsmax (for DRSDT) and similar.
@@ -920,164 +877,421 @@ private:
             }
         }
 
+        const double invDt = 1.0 / simulator_().timeStepSize();
+
         OPM_TIMEBLOCK(linearize);
 
         // We do not call resetSystem_() here, since that will set
         // the full system to zero, not just our part.
         // Instead, that must be called before starting the linearization.
         const bool dispersionActive = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion();
-        const auto& blockVelocity = simulator_().problem().eclWriter().outputModule().getFlows().blockVelocity();
         const unsigned int numCells = domain.cells.size();
+        const bool onFullDomain = (numCells == model_().numTotalDof());
+
+        if constexpr (!run_assembly_on_gpu) {
+            linearize_parallelization_wrapper<run_assembly_on_gpu, LocalResidual>(
+                numCells /*numCells*/,
+                domain,
+                neighborInfo_,
+                diagMatAddress_,
+                residual_,
+                model_(),
+                invDt,
+                dispersionActive,
+                onFullDomain,
+                problem_());
+
+            linearize_bc<IntensiveQuantities, Model, LocalResidual>(
+                diagMatAddress_, residual_, boundaryInfo_);
+        } else {
+#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+            // Make sure we have can have the domain on the GPU.
+            if constexpr (std::is_same_v<SubDomainType, FullDomain<>>) {
+
+                using GpuParams = TpfaLinearizerGpuParams<TypeTag>;
+
+                int constexpr blockSize = 256; // Experimentally this is a good value for multiple
+                                               // GPUs. Autotune this later.
+
+                GpuParams gpuParams(domain,
+                                    neighborInfo_,
+                                    *gpuJacobian_,
+                                    *gpuBufferDiagMatAddress_,
+                                    jacobian_->istlMatrix(),
+                                    residual_,
+                                    boundaryInfo_,
+                                    model_(),
+                                    problem_(),
+                                    numCells);
+
+                linearize_parallelization_wrapper<run_assembly_on_gpu,
+                                                  typename GpuParams::LocalResidualGPU>(
+                    numCells,
+                    gpuParams.domainView(),
+                    gpuParams.neighborInfoView(),
+                    gpuParams.diagMatAddressView(),
+                    gpuParams.residualView(),
+                    gpuParams.modelView(),
+                    invDt,
+                    dispersionActive,
+                    onFullDomain,
+                    gpuParams.flowProblemView());
+
+                if (gpuParams.boundaryInfoSize() > 0) {
+                    auto boundaryInfoView = gpuParams.boundaryInfoView();
+                    linearize_bc_threadsafe<TpfaLinearizer<TypeTag>,
+                                            typename GpuParams::GPUBOIQ,
+                                            decltype(gpuParams.modelView()),
+                                            typename GpuParams::LocalResidualGPU>
+                        <<<((gpuParams.boundaryInfoSize() + blockSize - 1) / blockSize),
+                           blockSize>>>(gpuParams.diagMatAddressView(),
+                                        gpuParams.residualView(),
+                                        boundaryInfoView,
+                                        gpuParams.modelView(),
+                                        gpuParams.flowProblemView());
+                }
 
-        // Fetch timestepsize used later in accumulation term.
-        const double dt = simulator_().timeStepSize();
+                // The memory copies here are synchronous and in the default stream,
+                // guaranteeing that the GPU kernels have completed.
+                gpuParams.copyResidualToHost(residual_, numCells);
+                gpuParams.copyJacobianToHost(*jacobian_, *gpuJacobian_);
+            } else {
+                OPM_THROW(std::logic_error, "Only FullDomain is supported on GPU");
+            }
+#else
+            OPM_THROW(std::logic_error,
+                      "Trying to run GPU assembly without compiling with GPU support");
+#endif
+        }
 
+        // Handle source terms separately as we want the functionality for CPU and GPU cases
+        // but for now we cannot handle this inside gpu kernels due to the use of problem_()
+        linearize_source_terms(numCells, domain);
+    }
+
+    template <bool useGPU,
+              class LocalResidualT,
+              class ModelClass,
+              class DiagPtrType,
+              class DomainType,
+              class NeighborSparseTable,
+              class ResidualType,
+              class ProblemT>
+    void linearize_parallelization_wrapper(const unsigned int numCells,
+                                           const DomainType& domain,
+                                           const NeighborSparseTable& neighborInfo,
+                                           DiagPtrType& diagMatAddress,
+                                           ResidualType& residual,
+                                           const ModelClass& model,
+                                           Scalar invLocDT,
+                                           bool dispersionActive,
+                                           bool onFullDomain,
+                                           const ProblemT& problem)
+    {
+        if constexpr (useGPU) {
+            static_assert(!enableBioeffects && "Bioeffects not yet supported on GPU");
+            assert(!dispersionActive && "Dispersion not yet supported on GPU");
+#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+            int constexpr blockSize = 256;
+            kernel_linearize<TpfaLinearizer<TypeTag>,
+                             ModelClass,
+                             LocalResidualT,
+                             DiagPtrType,
+                             Scalar,
+                             DomainType,
+                             NeighborSparseTable,
+                             ResidualType,
+                             ProblemT>
+                <<<((numCells + blockSize - 1) / blockSize), blockSize>>>(numCells,
+                                                                          domain,
+                                                                          neighborInfo,
+                                                                          diagMatAddress,
+                                                                          residual,
+                                                                          model,
+                                                                          invLocDT,
+                                                                          onFullDomain,
+                                                                          problem);
+#else
+            OPM_THROW(std::runtime_error, "Trying to run GPU code without GPU support");
+#endif
+        } else {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+            for (unsigned ii = 0; ii < numCells; ++ii) {
+                linearize_cell<false, LocalResidual>(ii,
+                                                     domain,
+                                                     neighborInfo,
+                                                     diagMatAddress,
+                                                     residual,
+                                                     model,
+                                                     velocityInfo_,
+                                                     invLocDT,
+                                                     problem,
+                                                     onFullDomain);
+            }
+        }
+    }
+
+    template <class SubDomainType>
+    void linearize_source_terms(unsigned int numCells, const SubDomainType& domain)
+    {
 #ifdef _OPENMP
 #pragma omp parallel for
 #endif
         for (unsigned ii = 0; ii < numCells; ++ii) {
             OPM_TIMEBLOCK_LOCAL(linearizationForEachCell, Subsystem::Assembly);
             const unsigned globI = domain.cells[ii];
-            const auto& nbInfos = neighborInfo_[globI];
-            VectorBlock res(0.0);
-            MatrixBlock bMat(0.0);
-            ADVectorBlock adres(0.0);
-            ADVectorBlock darcyFlux(0.0);
+            VectorBlockCPU res(0.0);
+            MatrixBlockCPU bMat(0.0);
+            ADVectorBlockCPU adres(0.0);
             const IntensiveQuantities& intQuantsIn = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
+            const double volume = model_().dofTotalVolume(globI);
 
-            // Flux term.
-            {
-                OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachCell, Subsystem::Assembly);
-                short loc = 0;
-                for (const auto& nbInfo : nbInfos) {
-                    OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachFace, Subsystem::Assembly);
-                    const unsigned globJ = nbInfo.neighbor;
-                    assert(globJ != globI);
-                    res = 0.0;
-                    bMat = 0.0;
-                    adres = 0.0;
-                    darcyFlux = 0.0;
-                    const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0);
-                    LocalResidual::computeFlux(adres, darcyFlux, globI, globJ, intQuantsIn, intQuantsEx,
-                                               nbInfo.res_nbinfo, problem_().moduleParams());
-                    adres *= nbInfo.res_nbinfo.faceArea;
-                    if (dispersionActive || enableBioeffects) {
+            // Cell-wise source terms.
+            // This will include well sources if SeparateSparseSourceTerms is false.
+            res = 0.0;
+            bMat = 0.0;
+            adres = 0.0;
+            if (separateSparseSourceTerms_) {
+                LocalResidual::computeSourceDense(adres, problem_(), intQuantsIn, globI, 0);
+            } else {
+                LocalResidual::computeSource(adres, problem_(), intQuantsIn, globI, 0);
+            }
+            adres *= -volume;
+            setResAndJacobi(res, bMat, adres);
+            residual_[globI] += res;
+            // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
+            *diagMatAddress_[globI] += bMat;
+        } // end of loop for cell globI.
+
+        // Add sparse source terms. For now only wells.
+        if (separateSparseSourceTerms_) {
+            problem_().wellModel().addReservoirSourceTerms(residual_, diagMatAddress_);
+        }
+    }
+
+public:
+    template <bool useGPU,
+              class LocalResidualT,
+              class ProblemType,
+              class VelocityInfoType,
+              class ModelClass,
+              class DiagPtrType,
+              class DomainType,
+              class NeighborSparseTable,
+              class ResidualType>
+    OPM_HOST_DEVICE static void linearize_cell(const unsigned int ii,
+                                               const DomainType& domain,
+                                               const NeighborSparseTable& neighborInfo,
+                                               const DiagPtrType& diagMatAddress,
+                                               ResidualType& residual,
+                                               const ModelClass& model,
+                                               VelocityInfoType& velocityInfo,
+                                               const Scalar invLocDT,
+                                               const ProblemType& problem,
+                                               const bool onFullDomain)
+    {
+#if OPM_IS_INSIDE_HOST_FUNCTION
+        OPM_TIMEBLOCK_LOCAL(linearizationForEachCell, Subsystem::Assembly);
+#endif
+        const unsigned globI = domain.cells[ii];
+        const auto& nbInfos = neighborInfo[globI];
+        VectorBlock res(0.0);
+        MatrixBlock bMat(0.0);
+        ADVectorBlock adres(0.0);
+        ADVectorBlock darcyFlux(0.0);
+        const auto& intQuantsIn = model.intensiveQuantities(globI, /*timeIdx*/ 0);
+
+        // Flux term.
+        {
+#if OPM_IS_INSIDE_HOST_FUNCTION
+            OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachCell, Subsystem::Assembly);
+#endif
+            short loc = 0;
+            for (const auto& nbInfo : nbInfos) {
+                const unsigned globJ = nbInfo.neighbor;
+                res = 0.0;
+                bMat = 0.0;
+                adres = 0.0;
+                darcyFlux = 0.0;
+
+                const auto& intQuantsEx = model.intensiveQuantities(globJ, /*timeIdx*/ 0);
+
+                LocalResidualT::computeFlux(adres,
+                                            darcyFlux,
+                                            globI,
+                                            globJ,
+                                            intQuantsIn,
+                                            intQuantsEx,
+                                            nbInfo.res_nbinfo,
+                                            problem.moduleParams());
+
+                adres *= nbInfo.res_nbinfo.faceArea;
+
+                // currently not supported on GPU
+                if constexpr (!useGPU) {
+                    if (velocityInfo.rowSize(globI) > 0) {
                         for (unsigned phaseIdx = 0; phaseIdx < numEq; ++phaseIdx) {
-                            velocityInfo_[globI][loc].velocity[phaseIdx] =
-                                darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea;
-                        }
-                    }
-                    else if (!blockVelocity.empty()) {
-                        if (std::ranges::binary_search(blockVelocity,
-                                                       simulator_().vanguard().cartesianIndex(globI))) {
-                            for (unsigned phaseIdx = 0; phaseIdx < numEq; ++phaseIdx) {
-                                velocityInfo_[globI][loc].velocity[phaseIdx] =
-                                    darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea;
-                            }
+                            velocityInfo[globI][loc].velocity[phaseIdx]
+                                = darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea;
                         }
                     }
-                    setResAndJacobi(res, bMat, adres);
-                    residual_[globI] += res;
-                    //SparseAdapter syntax:  jacobian_->addToBlock(globI, globI, bMat);
-                    *diagMatAddress_[globI] += bMat;
-                    bMat *= -1.0;
-                    //SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat);
-                    *nbInfo.matBlockAddress += bMat;
-                    ++loc;
                 }
-            }
 
-            // Accumulation term.
-            const double volume = model_().dofTotalVolume(globI);
-            const Scalar storefac = volume / dt;
-            adres = 0.0;
-            {
-                OPM_TIMEBLOCK_LOCAL(computeStorage, Subsystem::Assembly);
-                LocalResidual::template computeStorage<Evaluation>(adres, intQuantsIn);
+                setResAndJacobi(res, bMat, adres);
+                residual[globI] += res;
+                // SparseAdapter syntax:  jacobian_->addToBlock(globI, globI, bMat);
+                *diagMatAddress[globI] += bMat;
+                bMat *= -1.0;
+                // SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat);
+                *nbInfo.matBlockAddress += bMat;
+                ++loc;
             }
-            setResAndJacobi(res, bMat, adres);
+        }
+
+        // Accumulation term.
+        adres = 0.0;
+        {
+            LocalResidualT::template computeStorage<Evaluation>(adres, intQuantsIn);
+        }
+        setResAndJacobi(res, bMat, adres);
+
+        if constexpr (!useGPU) { // Cached storage not enabled for GPU
             // Either use cached storage term, or compute it on the fly.
-            if (model_().enableStorageCache()) {
+            if (model.enableStorageCache()) {
                 // The cached storage for timeIdx 0 (current time) is not
                 // used, but after storage cache is shifted at the end of the
                 // timestep, it will become cached storage for timeIdx 1.
-                model_().updateCachedStorage(globI, /*timeIdx=*/0, res);
+                model.updateCachedStorage(globI, /*timeIdx=*/0, res);
                 // We should not update the storage cache here for NLDD local solves.
                 // This will reset the start-of-step storage to incorrect numbers when
                 // we do local solves, where the iteration number will start from 0,
                 // but the starting state may not be identical to the start-of-step state.
                 // Note that a full assembly must be done before local solves
                 // otherwise this will be left un-updated.
-                if (problem_().iterationContext().isFirstGlobalIteration()) {
+                if (problem.iterationContext().isFirstGlobalIteration()) {
                     // Need to update the storage cache.
-                    if (problem_().recycleFirstIterationStorage()) {
+                    if (problem.recycleFirstIterationStorage()) {
                         // Assumes nothing have changed in the system which
                         // affects masses calculated from primary variables.
-                        model_().updateCachedStorage(globI, /*timeIdx=*/1, res);
-                    }
-                    else {
-                        Dune::FieldVector<Scalar, numEq> tmp;
-                        const IntensiveQuantities intQuantOld = model_().intensiveQuantities(globI, 1);
-                        LocalResidual::template computeStorage<Scalar>(tmp, intQuantOld);
-                        model_().updateCachedStorage(globI, /*timeIdx=*/1, tmp);
+                        if (onFullDomain) {
+                            // This is to avoid resetting the start-of-step storage
+                            // to incorrect numbers when we do local solves, where the iteration
+                            // number will start from 0, but the starting state may not be identical
+                            // to the start-of-step state.
+                            // Note that a full assembly must be done before local solves
+                            // otherwise this will be left un-updated.
+                            model.updateCachedStorage(globI, 1, res);
+                        }
+                    } else {
+                        VectorBlock tmp;
+                        const auto& intQuantOld = model.intensiveQuantities(globI, 1);
+                        LocalResidualT::template computeStorage<Scalar>(tmp, intQuantOld);
+                        model.updateCachedStorage(globI, 1, tmp);
                     }
                 }
-                res -= model_().cachedStorage(globI, 1);
-            }
-            else {
+                res -= model.cachedStorage(globI, 1);
+            } else {
+#if OPM_IS_INSIDE_HOST_FUNCTION
                 OPM_TIMEBLOCK_LOCAL(computeStorage0, Subsystem::Assembly);
-                Dune::FieldVector<Scalar, numEq> tmp;
-                const IntensiveQuantities intQuantOld = model_().intensiveQuantities(globI, 1);
-                LocalResidual::template computeStorage<Scalar>(tmp, intQuantOld);
+#endif
+                VectorBlock tmp;
+                const auto& intQuantOld = model.intensiveQuantities(globI, 1);
+                LocalResidualT::template computeStorage<Scalar>(tmp, intQuantOld);
                 // assume volume do not change
                 res -= tmp;
             }
-            res *= storefac;
-            bMat *= storefac;
-            residual_[globI] += res;
-            //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
-            *diagMatAddress_[globI] += bMat;
+        } else {
+            VectorBlock tmp;
+            const auto& intQuantOld = model.intensiveQuantities(globI, 1);
+            LocalResidualT::template computeStorage<Scalar>(tmp, intQuantOld);
+            // assume volume do not change
+            res -= tmp;
+        }
 
-            // Cell-wise source terms.
-            // This will include well sources if SeparateSparseSourceTerms is false.
-            res = 0.0;
-            bMat = 0.0;
-            adres = 0.0;
-            if (separateSparseSourceTerms_) {
-                LocalResidual::computeSourceDense(adres, problem_(), intQuantsIn, globI, 0);
+        const Scalar volume = model.dofTotalVolume(globI);
+        const Scalar storefac = volume * invLocDT; // volume / dt;
+        res *= storefac;
+        bMat *= storefac;
+        residual[globI] += res;
+        // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
+        *diagMatAddress[globI] += bMat;
+    }
+
+#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+    template <class LocalIntensiveQuantities,
+              class LocalResidualT,
+              class ModelT,
+              class DiagPtrType,
+              class ResidualT,
+              class BoundaryInfoT,
+              class ProblemT>
+    OPM_HOST_DEVICE static void linearize_bc_threadsafe_single_cell(DiagPtrType diagMatAdress,
+                                                               ResidualT residual,
+                                                               const BoundaryInfoT boundaryInfoElement,
+                                                               ModelT model,
+                                                               ProblemT problem)
+    {
+        using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+        constexpr int numEq = getPropValue<TypeTag, Properties::NumEq>();
+        if (boundaryInfoElement.bcdata.type != BCType::NONE) {
+            VectorBlock res(.0);
+            MatrixBlock bMat(0.0);
+            ADVectorBlock adres(0.0);
+            const unsigned globI = boundaryInfoElement.cell;
+            const LocalIntensiveQuantities& insideIntQuants
+                = model.intensiveQuantities(globI, /*timeIdx*/ 0);
+            if constexpr (!std::is_empty_v<GetPropType<TypeTag, Properties::FluidSystem>>) {
+                LocalResidualT::computeBoundaryFlux(
+                    adres, problem, boundaryInfoElement.bcdata, insideIntQuants, globI);
             }
-            else {
-                LocalResidual::computeSource(adres, problem_(), intQuantsIn, globI, 0);
+            adres *= boundaryInfoElement.bcdata.faceArea;
+            TpfaLinearizer<TypeTag>::setResAndJacobi(res, bMat, adres);
+            auto* residualPtr = &(residual.data()[globI]);
+            for (int i = 0; i < numEq; ++i) {
+                atomicAdd(&((*residualPtr)[i]), res[i]);
+            }
+            auto* matPtr = diagMatAdress[globI];
+            for (int row = 0; row < bMat.size(); ++row) {
+                for (int col = 0; col < bMat.size(); ++col) {
+                    Scalar* elemPtr = &((*matPtr)[row][col]);
+                    atomicAdd(elemPtr, bMat[row][col]);
+                }
             }
-            adres *= -volume;
-            setResAndJacobi(res, bMat, adres);
-            residual_[globI] += res;
-            //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
-            *diagMatAddress_[globI] += bMat;
-        } // end of loop for cell globI.
-
-        // Add sparse source terms. For now only wells.
-        if (separateSparseSourceTerms_) {
-            problem_().wellModel().addReservoirSourceTerms(residual_, diagMatAddress_);
         }
+    }
+#endif
 
+private:
+    template <class LocalIntensiveQuantities,
+              class ModelClass,
+              class LocalResidualT,
+              class DiagPtrType,
+              class ResidualType,
+              class BoundaryInfoT>
+    void linearize_bc(DiagPtrType& diagMatAdress,
+                      ResidualType& residual,
+                      const BoundaryInfoT& boundaryInfo)
+    {
         // Boundary terms. Only looping over cells with nontrivial bcs.
-        for (const auto& bdyInfo : boundaryInfo_) {
+        for (const auto& bdyInfo : boundaryInfo) {
             if (bdyInfo.bcdata.type == BCType::NONE) {
                 continue;
             }
-
-            VectorBlock res(0.0);
+            VectorBlock res(.0);
             MatrixBlock bMat(0.0);
             ADVectorBlock adres(0.0);
             const unsigned globI = bdyInfo.cell;
-            const IntensiveQuantities& insideIntQuants = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
-            LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, insideIntQuants, globI);
+            const LocalIntensiveQuantities& insideIntQuants
+                = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
+            LocalResidual::computeBoundaryFlux(
+                adres, problem_(), bdyInfo.bcdata, insideIntQuants, globI);
             adres *= bdyInfo.bcdata.faceArea;
             setResAndJacobi(res, bMat, adres);
-            residual_[globI] += res;
+            residual[globI] += res;
             ////SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
-            *diagMatAddress_[globI] += bMat;
+            *diagMatAdress[globI] += bMat;
         }
     }
 
@@ -1107,6 +1321,10 @@ private:
 
     // the jacobian matrix
     std::unique_ptr<SparseMatrixAdapter> jacobian_{};
+#if HAVE_CUDA
+    std::unique_ptr<gpuistl::GpuSparseMatrixWrapper<Scalar>> gpuJacobian_;
+    std::unique_ptr<gpuistl::GpuBuffer<MatrixBlockGPU*>> gpuBufferDiagMatAddress_;
+#endif
 
     // the right-hand side
     GlobalEqVector residual_;
@@ -1114,10 +1332,10 @@ private:
     LinearizationType linearizationType_{};
 
     using ResidualNBInfo = typename LocalResidual::ResidualNBInfo;
-    using NeighborInfoCPU = NeighborInfoStruct<ResidualNBInfo, MatrixBlock>;
+    using NeighborInfoCPU = NeighborInfoStruct<ResidualNBInfo, MatrixBlockCPU>;
 
     SparseTable<NeighborInfoCPU> neighborInfo_{};
-    std::vector<MatrixBlock*> diagMatAddress_{};
+    std::vector<MatrixBlockCPU*> diagMatAddress_ {};
 
     struct FlowInfo
     {
@@ -1136,25 +1354,12 @@ private:
     SparseTable<VelocityInfo> velocityInfo_;
 
     using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState;
-    struct BoundaryConditionData
-    {
-        BCType type;
-        VectorBlock massRate;
-        unsigned pvtRegionIdx;
-        unsigned boundaryFaceIndex;
-        double faceArea;
-        double faceZCoord;
-        ScalarFluidState exFluidState;
-    };
 
-    struct BoundaryInfo
-    {
-        unsigned int cell;
-        int dir;
-        unsigned int bfIndex;
-        BoundaryConditionData bcdata;
-    };
-    std::vector<BoundaryInfo> boundaryInfo_;
+    using BoundaryConditionDataCPU = BoundaryConditionData<VectorBlockCPU, ScalarFluidState>;
+
+    using BoundaryInfoCPU = BoundaryInfo<BoundaryConditionDataCPU>;
+
+    std::vector<BoundaryInfoCPU> boundaryInfo_;
 
     bool separateSparseSourceTerms_ = false;
 
diff --git a/opm/models/discretization/common/tpfalinearizergpukernels.hh b/opm/models/discretization/common/tpfalinearizergpukernels.hh
new file mode 100644
index 00000000000..c661c9cc0b9
--- /dev/null
+++ b/opm/models/discretization/common/tpfalinearizergpukernels.hh
@@ -0,0 +1,93 @@
+/*
+  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 TPFA_LINEARIZER_GPUKERNELS_HH
+#define TPFA_LINEARIZER_GPUKERNELS_HH
+
+/*
+    This file contains the GPU kernels that handles the GPU parallelization of the linearization.
+    This is effectively just the GPU alternative to what is a for-loop on the CPU
+*/
+
+namespace Opm
+{
+
+#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+template <class TpfaLinearizerType,
+          class ModelClass,
+          class LocalResidualT,
+          class DiagPtrType,
+          class ScalarType,
+          class DomainType,
+          class NeighborSparseTable,
+          class ResidualType,
+          class ProblemT>
+__global__ __launch_bounds__(256) void kernel_linearize(const unsigned int numCells,
+                                                        const DomainType domain,
+                                                        const NeighborSparseTable neighborInfo,
+                                                        DiagPtrType diagMatAddress,
+                                                        ResidualType residual,
+                                                        ModelClass model,
+                                                        ScalarType invLocDT,
+                                                        bool onFullDomain,
+                                                        ProblemT problem)
+{
+    const unsigned int ii = blockIdx.x * blockDim.x + threadIdx.x;
+
+    if (ii < numCells) {
+        // velocityInfo is unused on GPU (guarded by if constexpr (!useGPU)),
+        // but the parameter is T& so we need an lvalue to bind to.
+        std::nullptr_t dummyVelocityInfo = nullptr;
+        TpfaLinearizerType::template linearize_cell<true, LocalResidualT>(ii,
+                                                                          domain,
+                                                                          neighborInfo,
+                                                                          diagMatAddress,
+                                                                          residual,
+                                                                          model,
+                                                                          dummyVelocityInfo,
+                                                                          invLocDT,
+                                                                          problem,
+                                                                          onFullDomain);
+    }
+}
+
+template <class TpfaLinearizerType,
+          class IntensiveQuantities,
+          class ModelT,
+          class LocalResidualT,
+          class DiagPtrType,
+          class ResidualT,
+          class BoundaryInfo,
+          class ProblemT>
+__global__ void
+linearize_bc_threadsafe(DiagPtrType diagMatAddress,
+                        ResidualT residual,
+                        const BoundaryInfo boundaryInfo,
+                        ModelT model,
+                        ProblemT gpuProblem)
+{
+    const unsigned int ii = blockIdx.x * blockDim.x + threadIdx.x;
+    if (ii < boundaryInfo.size()) {
+        TpfaLinearizerType::template linearize_bc_threadsafe_single_cell<IntensiveQuantities,
+                                                                         LocalResidualT>(
+            diagMatAddress, residual, boundaryInfo[ii], model, gpuProblem);
+    }
+}
+
+#endif
+
+} // namespace Opm
+
+#endif // TPFA_LINEARIZER_GPUKERNELS_HH
diff --git a/opm/models/discretization/common/tpfalinearizergpuparams.hh b/opm/models/discretization/common/tpfalinearizergpuparams.hh
new file mode 100644
index 00000000000..07e0ada5654
--- /dev/null
+++ b/opm/models/discretization/common/tpfalinearizergpuparams.hh
@@ -0,0 +1,300 @@
+// -*- 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/>.
+*/
+/*!
+ * \file
+ *
+ * \brief GPU parameter setup class for TpfaLinearizer.
+ *
+ * \note This file is designed to be included from tpfalinearizer.hh inside a
+ *       \c #if HAVE_CUDA block.  It is self-contained: all required GPU headers,
+ *       model headers, and forward declarations are included below.
+ */
+#ifndef TPFA_LINEARIZER_GPU_PARAMS_HH
+#define TPFA_LINEARIZER_GPU_PARAMS_HH
+
+#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#include <opm/common/utility/gpuistl_if_available.hpp>
+
+#if USE_HIP
+#include <opm/simulators/linalg/gpuistl_hip/GpuSparseMatrixWrapper.hpp>
+#include <opm/simulators/linalg/gpuistl_hip/MiniMatrix.hpp>
+#include <opm/simulators/linalg/gpuistl_hip/MiniVector.hpp>
+#else
+#include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
+#include <opm/simulators/linalg/gpuistl/MiniMatrix.hpp>
+#include <opm/simulators/linalg/gpuistl/MiniVector.hpp>
+#endif
+
+#include <opm/grid/utility/SparseTable.hpp>
+#include <opm/models/blackoil/blackoilintensivequantities.hh>
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
+#include <opm/simulators/flow/SimpleFIBlackOilModel.hpp>
+#include <opm/simulators/flow/ThermalGasWaterFlowProblem.hpp>
+
+#include <opm/models/discretization/common/tpfalinearizerstructs.hh>
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Manages GPU buffer allocation and view setup for one linearization pass.
+ *
+ * Constructed once per call to TpfaLinearizer::linearize_() when the GPU assembly
+ * path is active.  The constructor performs all CPU-to-GPU copies and stores the
+ * resulting GPU buffers as members.  Callers retrieve lightweight GpuView objects
+ * via the accessor methods and pass them directly to the GPU kernels.  After the
+ * kernels complete, copyResidualToHost() and copyJacobianToHost() transfer the
+ * computed results back to the CPU.
+ */
+template <class TypeTag>
+class TpfaLinearizerGpuParams
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
+    using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
+    using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Model = GetPropType<TypeTag, Properties::Model>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+
+    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
+
+    using MatrixBlockCPU = typename SparseMatrixAdapter::MatrixBlock;
+    using VectorBlockCPU = Dune::FieldVector<Scalar, numEq>;
+    using MatrixBlockGPU = gpuistl::MiniMatrix<Scalar, numEq>;
+    using VectorBlockGPU = gpuistl::MiniVector<Scalar, numEq>;
+
+    using ResidualNBInfo = typename LocalResidual::ResidualNBInfo;
+    using NeighborInfoCPU = NeighborInfoStruct<ResidualNBInfo, MatrixBlockCPU>;
+    using NeighborInfoGPU = NeighborInfoStruct<ResidualNBInfo, MatrixBlockGPU>;
+
+    using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState;
+    using BoundaryConditionDataCPU = BoundaryConditionData<VectorBlockCPU, ScalarFluidState>;
+    using BoundaryInfoCPU = BoundaryInfo<BoundaryConditionDataCPU>;
+
+    // Fluid system buffer/view/ptr types derived without spelling out the fluid-system
+    // template parameters explicitly.
+    using NonStaticFluidSystem = std::decay_t<decltype(FluidSystem::getNonStaticInstance())>;
+    using GpuFluidSystemBuffer = std::decay_t<decltype(::Opm::gpuistl::copy_to_gpu(
+        std::declval<NonStaticFluidSystem&>()))>;
+    using GpuFluidSystemView
+        = std::decay_t<decltype(::Opm::gpuistl::make_view(std::declval<GpuFluidSystemBuffer&>()))>;
+    using GpuFluidSystemPtr = std::shared_ptr<GpuFluidSystemView>;
+
+public:
+    // Type aliases required by callers to instantiate the GPU kernel templates.
+    using CorrectTypeTagView =
+        typename ::Opm::Properties::TTag::to_gpu_type_t<TypeTag, gpuistl::GpuView>;
+
+    // GPUBOIQ and LocalResidualGPU are not taken from typetag as we need the
+    // outer class to template correctly here
+    using GPUBOIQ = BlackOilIntensiveQuantities<CorrectTypeTagView>;
+    using LocalResidualGPU = BlackOilLocalResidualTPFA<CorrectTypeTagView>;
+
+private:
+    using GpuScalarFluidState = typename GPUBOIQ::ScalarFluidState;
+    using BoundaryConditionDataGPU = BoundaryConditionData<VectorBlockGPU, GpuScalarFluidState>;
+    using BoundaryInfoGPU = BoundaryInfo<BoundaryConditionDataGPU>;
+
+    using GpuModelBufferType = SimpleFIBlackOilModel<CorrectTypeTagView, gpuistl::GpuBuffer>;
+    using GpuFlowProblemBufferType = ThermalGasWaterFlowProblem<Scalar, gpuistl::GpuBuffer>;
+
+    using GpuModel = GetPropType<TypeTag, Properties::GpuFIBlackOilModel>;
+
+    FullDomain<gpuistl::GpuBuffer<int>> domainBuffer_;
+    SparseTable<NeighborInfoGPU, gpuistl::GpuBuffer> neighborInfoBuffer_;
+    // diagMatAddressView_ is non-owning: the underlying buffer lives in TpfaLinearizer.
+    gpuistl::GpuView<MatrixBlockGPU*> diagMatAddressView_;
+    gpuistl::GpuBuffer<VectorBlockGPU> residualBuffer_;
+    gpuistl::GpuView<VectorBlockGPU>
+        residualView_; // stored as lvalue because mutable ref must be provided
+    // dynamicGpuFluidSystemBuffer_ must be declared before dynamicGpuFluidSystemPtr_
+    // because the ptr holds a GpuView into the buffer's GPU memory.
+    GpuFluidSystemBuffer dynamicGpuFluidSystemBuffer_;
+    // The fluid-system ptr is kept alive because gpuModelBuffer_ and
+    // boundaryInfoBuffer_ store raw GPU pointers into it.
+    GpuFluidSystemPtr dynamicGpuFluidSystemPtr_;
+    GpuModelBufferType gpuModelBuffer_;
+    GpuFlowProblemBufferType gpuFlowProblemBuffer_;
+    gpuistl::GpuBuffer<BoundaryInfoGPU> boundaryInfoBuffer_;
+
+public:
+    /*!
+     * \brief Construct and populate all GPU buffers needed for one linearization pass.
+        To avoid optionals we use initializerlists, although some are a bit ugly
+        with the lambdas
+     *
+     * \param domain                  The full domain (CPU).
+     * \param neighborInfo            CPU neighbor-info sparse table.
+     * \param gpuJacobian             Persistent GPU sparse-matrix wrapper (owned by
+     * TpfaLinearizer).
+     * \param gpuBufferDiagMatAddress Persistent GPU diagonal-pointer buffer (owned by
+     * TpfaLinearizer).
+     * \param cpuJacobian             CPU ISTL matrix (used to compute GPU pointer offsets).
+     * \param residual                CPU residual vector (read to initialise the GPU residual
+     * buffer).
+     * \param boundaryInfo            CPU boundary-condition info.
+     * \param model                   CPU model (provides intensive quantities and volumes).
+     * \param problem                 CPU problem (provides thermal transmissibilities and module
+     * params).
+     * \param numCells                Number of cells in the domain.
+     */
+    TpfaLinearizerGpuParams(const FullDomain<>& domain,
+                            const SparseTable<NeighborInfoCPU>& neighborInfo,
+                            gpuistl::GpuSparseMatrixWrapper<Scalar>& gpuJacobian,
+                            gpuistl::GpuBuffer<MatrixBlockGPU*>& gpuBufferDiagMatAddress,
+                            typename SparseMatrixAdapter::IstlMatrix& cpuJacobian,
+                            GlobalEqVector& residual,
+                            const std::vector<BoundaryInfoCPU>& boundaryInfo,
+                            Model& model,
+                            const Problem& problem,
+                            unsigned int numCells)
+        : domainBuffer_(copy_to_gpu(domain))
+        , neighborInfoBuffer_(
+              gpuistl::copy_to_gpu<MatrixBlockGPU>(neighborInfo, gpuJacobian, cpuJacobian))
+        , diagMatAddressView_(gpuistl::make_view(gpuBufferDiagMatAddress))
+        , residualBuffer_(gpuistl::copy_to_gpu_residual<GlobalEqVector, VectorBlockGPU>(residual))
+        , residualView_(gpuistl::make_view(residualBuffer_))
+        , dynamicGpuFluidSystemBuffer_(
+              ::Opm::gpuistl::copy_to_gpu(FluidSystem::getNonStaticInstance()))
+        , dynamicGpuFluidSystemPtr_(
+              gpuistl::make_gpu_shared_ptr(::Opm::gpuistl::make_view(dynamicGpuFluidSystemBuffer_)))
+        , gpuModelBuffer_([&]() -> GpuModelBufferType {
+            std::vector<Scalar> volumes(numCells);
+            for (unsigned i = 0; i < numCells; ++i) {
+                volumes[domain.cells[i]] = model.dofTotalVolume(domain.cells[i]);
+            }
+            GpuModel gpuModel(
+                model.intensiveQuantityCache()[0], model.intensiveQuantityCache()[1], volumes);
+            return gpuistl::copy_to_gpu(gpuModel, *dynamicGpuFluidSystemPtr_.get());
+        }())
+        , gpuFlowProblemBuffer_([&]() -> GpuFlowProblemBufferType {
+            std::vector<Scalar> alpha0(numCells);
+            std::vector<Scalar> alpha1(numCells);
+            std::vector<Scalar> alpha2(numCells);
+            const auto& thermalHalfTransBoundary
+                = problem.eclTransmissibilities().getThermalHalfTransBoundary();
+            for (const auto& [key, value] : thermalHalfTransBoundary) {
+                const unsigned cell = key.first;
+                const unsigned dir = key.second;
+                if (dir == 0) {
+                    alpha0[cell] = value;
+                } else if (dir == 1) {
+                    alpha1[cell] = value;
+                } else if (dir == 2) {
+                    alpha2[cell] = value;
+                } else {
+                    OPM_THROW(std::logic_error,
+                              "Invalid direction for thermal half transmissibility: "
+                                  + std::to_string(dir));
+                }
+            }
+            ThermalGasWaterFlowProblem<Scalar> gpuFlowProblem(
+                alpha0, alpha1, alpha2, problem.moduleParams());
+            return gpuistl::copy_to_gpu(gpuFlowProblem);
+        }())
+        , boundaryInfoBuffer_(
+              gpuistl::copy_to_gpu<VectorBlockGPU, GpuScalarFluidState, BoundaryInfoGPU>(
+                  boundaryInfo, *dynamicGpuFluidSystemPtr_.get()))
+    {
+    }
+
+    auto domainView()
+    {
+        return make_view(domainBuffer_);
+    }
+
+    auto neighborInfoView()
+    {
+        return gpuistl::make_view(neighborInfoBuffer_);
+    }
+
+    auto& diagMatAddressView()
+    {
+        return diagMatAddressView_;
+    }
+
+    auto& residualView()
+    {
+        return residualView_;
+    }
+
+    auto modelView()
+    {
+        return gpuistl::make_view(gpuModelBuffer_);
+    }
+
+    auto flowProblemView()
+    {
+        return gpuistl::make_view(gpuFlowProblemBuffer_);
+    }
+
+    auto boundaryInfoView()
+    {
+        return gpuistl::make_view(boundaryInfoBuffer_);
+    }
+
+    std::size_t boundaryInfoSize() const
+    {
+        return boundaryInfoBuffer_.size();
+    }
+
+    /*!
+     * \brief Copy the GPU residual buffer back to the CPU residual vector.
+     *
+     * Must be called after all GPU kernels have completed (the synchronous cudaMemcpy
+     * inside guarantees the GPU has finished).
+     */
+    void copyResidualToHost(GlobalEqVector& residual, unsigned numCells)
+    {
+        auto cpuResidualFromGpu = residualBuffer_.asStdVector();
+        std::memcpy(residual.data(), cpuResidualFromGpu.data(), numCells * numEq * sizeof(Scalar));
+    }
+
+    /*!
+     * \brief Copy the GPU Jacobian non-zero values back to the CPU ISTL matrix.
+     *
+     * Must be called after all GPU kernels have completed.
+     */
+    void copyJacobianToHost(SparseMatrixAdapter& jacobian,
+                            gpuistl::GpuSparseMatrixWrapper<Scalar>& gpuJacobian)
+    {
+        auto gpuJacobianNonZeroes = gpuJacobian.getNonZeroValues().asStdVector();
+        auto& cpuJacobian = jacobian.istlMatrix();
+        const std::size_t totalMatrixSize = cpuJacobian.nonzeroes() * numEq * numEq;
+        std::memcpy(&(cpuJacobian[0][0][0][0]),
+                    gpuJacobianNonZeroes.data(),
+                    totalMatrixSize * sizeof(Scalar));
+    }
+};
+
+} // namespace Opm
+
+#endif // HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+
+#endif // TPFA_LINEARIZER_GPU_PARAMS_HH
diff --git a/opm/models/discretization/common/tpfalinearizerstructs.hh b/opm/models/discretization/common/tpfalinearizerstructs.hh
new file mode 100644
index 00000000000..eadefba997a
--- /dev/null
+++ b/opm/models/discretization/common/tpfalinearizerstructs.hh
@@ -0,0 +1,206 @@
+/*
+  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 TPFA_LINEARIZER_STRUCTS_HH
+#define TPFA_LINEARIZER_STRUCTS_HH
+
+/*!
+ * Structs needed for tpfalinearizer and its gpuparams struct extracted
+ * to be defined in one place that both can include.
+ */
+
+#include <cassert>
+#include <vector>
+
+#include <opm/common/Exceptions.hpp>
+#include <opm/grid/utility/SparseTable.hpp>
+#include <opm/input/eclipse/Schedule/BCProp.hpp>
+
+#include <opm/common/utility/gpuistl_if_available.hpp>
+#include <opm/common/utility/pointerArithmetic.hpp>
+
+namespace Opm {
+
+template <class Storage = std::vector<int>>
+struct FullDomain
+{
+    Storage cells;
+};
+
+#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+inline FullDomain<gpuistl::GpuBuffer<int>> copy_to_gpu(FullDomain<> CPUDomain)
+{
+    if (CPUDomain.cells.size() == 0) {
+        OPM_THROW(std::runtime_error, "Cannot copy empty full domain to GPU.");
+    }
+    return FullDomain<gpuistl::GpuBuffer<int>>{
+        gpuistl::GpuBuffer<int>(CPUDomain.cells)
+    };
+};
+
+inline FullDomain<gpuistl::GpuView<int>> make_view(FullDomain<gpuistl::GpuBuffer<int>>& buffer)
+{
+    if (buffer.cells.size() == 0) {
+        OPM_THROW(std::runtime_error, "Cannot make view of empty full domain buffer.");
+    }
+    return FullDomain<gpuistl::GpuView<int>>{
+        gpuistl::make_view(buffer.cells)
+    };
+};
+#endif
+
+template <class ResidualNBInfoType, class BlockType>
+struct NeighborInfoStruct
+{
+    unsigned int neighbor;
+    ResidualNBInfoType res_nbinfo;
+    BlockType* matBlockAddress;
+
+    template <class PtrType>
+    NeighborInfoStruct(unsigned int n, const ResidualNBInfoType& r, PtrType ptr)
+        : neighbor(n)
+        , res_nbinfo(r)
+        , matBlockAddress(static_cast<BlockType*>(ptr))
+    {
+    }
+
+    // Add a default constructor
+    NeighborInfoStruct()
+        : neighbor(0)
+        , res_nbinfo()
+        , matBlockAddress(nullptr)
+    {
+    }
+};
+
+template <class VectorBlock, class ScalarFluidState>
+struct BoundaryConditionData {
+    BCType type;
+    VectorBlock massRate;
+    unsigned pvtRegionIdx;
+    unsigned boundaryFaceIndex;
+    double faceArea;
+    double faceZCoord;
+    ScalarFluidState exFluidState;
+};
+
+template <class BoundaryConditionData>
+struct BoundaryInfo {
+    unsigned int cell;
+    int dir;
+    unsigned int bfIndex;
+    BoundaryConditionData bcdata;
+};
+
+#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+namespace gpuistl {
+
+    template <class MiniMatrixType,
+              class GpuMatrixType,
+              class CpuMatrixType,
+              class MatrixBlockType,
+              class ResidualNBInfoType>
+    auto copy_to_gpu(const SparseTable<NeighborInfoStruct<ResidualNBInfoType, MatrixBlockType>>&
+                         cpuNeighborInfoTable,
+                     GpuMatrixType& gpuJacobian,
+                     CpuMatrixType& cpuJacobian)
+    {
+        // Convert the DUNE FieldVectors to MiniMatrix types
+        using StructWithMinimatrix = NeighborInfoStruct<ResidualNBInfoType, MiniMatrixType>;
+        using Scalar = typename GpuMatrixType::field_type;
+        std::vector<StructWithMinimatrix> minimatrices(cpuNeighborInfoTable.dataSize());
+        Scalar* gpuBufStart = gpuJacobian.getNonZeroValues().data();
+        Scalar* cpuBufStart = &(cpuJacobian[0][0][0][0]);
+
+        // To compute the length of the buffer of the cpuJacobian we here assume we have a blocked
+        // BCRS matrix with square blocks and that the blocks are stored as Dune::FieldMatrix
+        using CpuBlockType = typename CpuMatrixType::block_type::BaseType;
+
+        const size_t gpuBufferSize
+            = gpuJacobian.nonzeroes() * gpuJacobian.blockSize() * gpuJacobian.blockSize();
+        const size_t cpuBufferSize
+            = cpuJacobian.nonzeroes() * CpuBlockType::rows * CpuBlockType::cols;
+        assert(gpuBufferSize == cpuBufferSize);
+
+        const auto& dataStorage = cpuNeighborInfoTable.dataStorage();
+        const size_t dataSize = cpuNeighborInfoTable.dataSize();
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+        for (size_t idx = 0; idx < dataSize; ++idx) {
+            const auto& e = dataStorage[idx];
+            Scalar* cpuPtr = &((*e.matBlockAddress)[0][0]);
+
+            // convert the pointer from CPU to GPU pointer based on offset in CPU jacobian
+            Scalar* gpuPtr = ComputePtrBasedOnOffsetInOtherBuffer(
+                gpuBufStart, gpuBufferSize, cpuBufStart, cpuBufferSize, cpuPtr);
+
+            minimatrices[idx] = StructWithMinimatrix(
+                e.neighbor, e.res_nbinfo, reinterpret_cast<MiniMatrixType*>(gpuPtr));
+        }
+
+        return SparseTable<StructWithMinimatrix, gpuistl::GpuBuffer>(
+            gpuistl::GpuBuffer<StructWithMinimatrix>(minimatrices),
+            gpuistl::GpuBuffer<int>(cpuNeighborInfoTable.rowStarts())
+        );
+    }
+
+    // Handle the BoundaryInfo structs
+    template <class GpuVecBlock,
+              class GpuFluidState,
+              class BoundaryInfoTypeGPU,
+              class BoundaryInfoTypeCPU,
+              typename GpuFluidSystemPtr>
+    auto copy_to_gpu(const std::vector<BoundaryInfoTypeCPU>& cpu_boundary_info,
+                     const GpuFluidSystemPtr& ref)
+    {
+        std::vector<BoundaryInfoTypeGPU> gpu_boundary_info;
+        for (const auto& info : cpu_boundary_info) {
+            gpu_boundary_info.push_back(BoundaryInfoTypeGPU {
+                info.cell,
+                info.dir,
+                info.bfIndex,
+                BoundaryConditionData<GpuVecBlock, GpuFluidState> {
+                    info.bcdata.type,
+                    GpuVecBlock(info.bcdata.massRate),
+                    info.bcdata.pvtRegionIdx,
+                    info.bcdata.boundaryFaceIndex,
+                    info.bcdata.faceArea,
+                    info.bcdata.faceZCoord,
+                    info.bcdata.exFluidState.withOtherFluidSystem(ref)}});
+        }
+
+        return gpuistl::GpuBuffer<BoundaryInfoTypeGPU>(gpu_boundary_info);
+    }
+
+    // Implemented for residual_, which is a vector of FieldVectors.
+    // We then make a GpuBuffer of MiniVectors.
+    template <class CPUResidualType, class GpuMiniVector>
+    auto copy_to_gpu_residual(CPUResidualType& residual)
+    {
+        std::vector<GpuMiniVector> vectorOfMiniVectors;
+        for (const auto& minivec : residual) {
+            vectorOfMiniVectors.emplace_back(minivec);
+        }
+
+        return gpuistl::GpuBuffer<GpuMiniVector>(vectorOfMiniVectors);
+    }
+
+} // namespace gpuistl
+#endif
+
+} // namespace Opm
+
+#endif // TPFA_LINEARIZER_STRUCTS_HH
diff --git a/opm/simulators/flow/FIBlackoilModel.hpp b/opm/simulators/flow/FIBlackoilModel.hpp
index 5a56ed016f8..ff70053ae36 100644
--- a/opm/simulators/flow/FIBlackoilModel.hpp
+++ b/opm/simulators/flow/FIBlackoilModel.hpp
@@ -42,6 +42,15 @@
 
 #include <opm/material/fluidmatrixinteractions/EclMultiplexerMaterialParams.hpp>
 
+#if HAVE_CUDA
+#include <opm/simulators/flow/FlowProblemParameters.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuBlackoilIntensiveQuantitiesDispatcher.hpp>
+#include <memory>
+#include <variant>
+#include <vector>
+#endif
+
+#include <chrono>
 #include <cstddef>
 #include <stdexcept>
 #include <type_traits>
@@ -79,6 +88,10 @@ class FIBlackOilModel : public BlackOilModel<TypeTag>
                           Dune::Partitions::all,
                           ThreadManager::maxThreads())
     {
+#if HAVE_CUDA
+        useGpuIntensiveQuantitiesDispatcher_ =
+            Parameters::Get<Parameters::ExperimentalComputePropertiesOnGpu>();
+#endif
     }
 
     void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx) const
@@ -100,6 +113,22 @@ class FIBlackOilModel : public BlackOilModel<TypeTag>
                     elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
                 }
             }
+
+#if HAVE_CUDA
+            // After all cells are CPU-updated and written into the cache,
+            // overlay the GPU-computed BlackOil fields in one batched call.
+            if constexpr (Opm::gpuistl::GpuBlackoilIntensiveQuantitiesDispatcherSupport<TypeTag>::value)
+            {
+                if (useGpuIntensiveQuantitiesDispatcher_) {
+                    runGpuIntensiveQuantitiesDispatcher_(timeIdx);
+
+                    const std::size_t numCells = this->intensiveQuantityCache_[timeIdx].size();
+                    for (std::size_t i = 0; i < numCells; ++i) {
+                        this->setIntensiveQuantitiesCacheEntryValidity(i, timeIdx, true);
+                    }
+                }
+            }
+#endif
             OPM_END_PARALLEL_TRY_CATCH("invalidateAndUpdateIntensiveQuantities: state error",
                                        this->simulator_.vanguard().grid().comm());
         } else {
@@ -262,6 +291,21 @@ class FIBlackOilModel : public BlackOilModel<TypeTag>
                 this->template updateSingleCachedIntQuantUnchecked<Args...>(elementMapper.index(elem), timeIdx);
             }
         }
+
+#if HAVE_CUDA
+        // After the CPU per-cell update has populated all fields, optionally
+        // overlay the BlackOil intensive-quantities fields with their GPU
+        // counterparts via the experimental dispatcher. The dispatcher only
+        // overwrites the subset of fields covered by
+        // BlackOilIntensiveQuantities::overlayBlackOilFieldsFrom; everything
+        // else (mobility, energy, ...) keeps the CPU-computed value.
+        if constexpr (Opm::gpuistl::GpuBlackoilIntensiveQuantitiesDispatcherSupport<TypeTag>::value)
+        {
+            if (useGpuIntensiveQuantitiesDispatcher_) {
+                runGpuIntensiveQuantitiesDispatcher_(timeIdx);
+            }
+        }
+#endif
     }
 
     template <class ...Args>
@@ -276,6 +320,40 @@ class FIBlackOilModel : public BlackOilModel<TypeTag>
     }
 
     ElementChunks<GridView, Dune::Partitions::All> element_chunks_;
+
+#if HAVE_CUDA
+    bool useGpuIntensiveQuantitiesDispatcher_{false};
+    using GpuDispatcherStorage = std::conditional_t<
+        Opm::gpuistl::GpuBlackoilIntensiveQuantitiesDispatcherSupport<TypeTag>::value,
+        std::unique_ptr<Opm::gpuistl::GpuBlackoilIntensiveQuantitiesDispatcher<TypeTag>>,
+        std::monostate>;
+    mutable GpuDispatcherStorage gpuIntensiveQuantitiesDispatcher_{};
+
+    void runGpuIntensiveQuantitiesDispatcher_(const unsigned timeIdx) const
+    {
+        if constexpr (Opm::gpuistl::GpuBlackoilIntensiveQuantitiesDispatcherSupport<TypeTag>::value) {
+            if (!gpuIntensiveQuantitiesDispatcher_) {
+                gpuIntensiveQuantitiesDispatcher_ =
+                    std::make_unique<
+                        Opm::gpuistl::GpuBlackoilIntensiveQuantitiesDispatcher<TypeTag>>();
+            }
+            using PV = GetPropType<TypeTag, Properties::PrimaryVariables>;
+            const auto& sol = this->solution(timeIdx);
+            const std::size_t numCells = this->intensiveQuantityCache_[timeIdx].size();
+            std::vector<const PV*> priVarsPtrs(numCells);
+            std::vector<IntensiveQuantities*> outIqPtrs(numCells);
+            for (std::size_t i = 0; i < numCells; ++i) {
+                priVarsPtrs[i] = &sol[i];
+                outIqPtrs[i] = &this->intensiveQuantityCache_[timeIdx][i];
+            }
+            gpuIntensiveQuantitiesDispatcher_->update(
+                this->simulator_.problem(),
+                priVarsPtrs.data(),
+                outIqPtrs.data(),
+                numCells);
+        }
+    }
+#endif
 };
 
 } // namespace Opm
diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp
index 3f24a99d915..c88dbd9e9e9 100644
--- a/opm/simulators/flow/FlowProblem.hpp
+++ b/opm/simulators/flow/FlowProblem.hpp
@@ -833,6 +833,9 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     std::shared_ptr<const EclMaterialLawManager> materialLawManager() const
     { return materialLawManager_; }
 
+    std::shared_ptr<const EclThermalLawManager> thermalLawManager() const
+    { return thermalLawManager_; }
+
     template <class FluidState, class ...Args>
     void updateRelperms(
         std::array<Evaluation,numPhases> &mobility,
diff --git a/opm/simulators/flow/FlowProblemParameters.cpp b/opm/simulators/flow/FlowProblemParameters.cpp
index 2f5f29668da..61d8375a8d4 100644
--- a/opm/simulators/flow/FlowProblemParameters.cpp
+++ b/opm/simulators/flow/FlowProblemParameters.cpp
@@ -81,6 +81,13 @@ void registerFlowProblemParameters()
         ("Conserve inner energy and not enthalpy "
          "even if THERMAL is used.");
 
+#if HAVE_CUDA
+    Parameters::Register<Parameters::ExperimentalComputePropertiesOnGpu>
+        ("Experimental: compute BlackOilIntensiveQuantities on the GPU "
+         "via the GpuBlackoilIntensiveQuantitiesDispatcher. Only takes "
+         "effect for CO2STORE-compatible TypeTags; ignored otherwise.");
+#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..aad64841cfb 100644
--- a/opm/simulators/flow/FlowProblemParameters.hpp
+++ b/opm/simulators/flow/FlowProblemParameters.hpp
@@ -69,6 +69,14 @@ 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; };
 
+#if HAVE_CUDA
+// 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); ignored
+// otherwise.
+struct ExperimentalComputePropertiesOnGpu { static constexpr bool value = true; };
+#endif
+
 } // namespace Opm::Parameters
 
 namespace Opm {
diff --git a/opm/simulators/flow/GpuFlowProblem.hpp b/opm/simulators/flow/GpuFlowProblem.hpp
new file mode 100644
index 00000000000..77130dcc783
--- /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 <http://www.gnu.org/licenses/>.
+*/
+/*!
+ * \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 <opm/common/utility/VectorWithDefaultAllocator.hpp>
+#include <opm/common/utility/gpuDecorators.hpp>
+
+#include <opm/models/discretization/common/linearizationtype.hh>
+
+#include <opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp>
+#include <opm/material/thermal/GpuEclThermalLawManager.hpp>
+
+#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 {
+
+/*!
+ * \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 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;
+
+    /*! \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<typename EclThermalLawManager::FluidSystem, void>;
+
+    /*! \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<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)
+        : 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 <class CpuProblem,
+              class StS = Storage<Scalar>,
+              std::enable_if_t<
+                  std::is_same_v<StS, ::Opm::VectorWithDefaultAllocator<Scalar>>,
+                  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<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));
+            }
+        }
+    }
+
+    /*!
+     * \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<Scalar>, then uploads each as a single
+     * \c GpuBuffer<Scalar>. 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 <class CpuProblem,
+              class StS = Storage<Scalar>,
+              std::enable_if_t<
+                  std::is_same_v<StS, ::Opm::gpuistl::GpuBuffer<Scalar>>,
+                  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<int>(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<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];
+    }
+
+    /*! \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<unsigned>(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<unsigned>(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<unsigned>(elemIdx));
+    }
+
+    /*! \brief Default rock-pore-volume multiplier (1 if no compressibility). */
+    template <class Evaluation>
+    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 <class Evaluation>
+    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 <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);
+    }
+
+    /*! \name Direct accessors for serialization to the GPU. */
+    //!\{
+    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);
+        }
+    }
+
+    /*! \brief Build the per-cell rockFraction storage; empty when the
+     *         GpuFlowProblem instantiation has no thermal support. */
+    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);
+        }
+    }
+
+    /*! \brief Build the per-cell pvtRegionIndex storage; empty when the
+     *         GpuFlowProblem instantiation has no thermal support. */
+    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);
+            }
+        }
+    }
+
+    /*! \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 <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 {
+
+/*!
+ * \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 <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()));
+}
+
+/*!
+ * \brief Make a non-owning GpuView based GpuFlowProblem from an owning
+ *        GpuBuffer based GpuFlowProblem.
+ */
+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()));
+}
+
+} // namespace Opm::gpuistl
+
+#endif // OPM_GPU_FLOW_PROBLEM_HPP
diff --git a/opm/simulators/flow/NewTranFluxModule.hpp b/opm/simulators/flow/NewTranFluxModule.hpp
index e4105e6c20a..8c9f2bd9e20 100644
--- a/opm/simulators/flow/NewTranFluxModule.hpp
+++ b/opm/simulators/flow/NewTranFluxModule.hpp
@@ -36,6 +36,7 @@
 
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/utility/gpuDecorators.hpp>
 
 #include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
 
@@ -48,7 +49,6 @@
 #include <opm/models/blackoil/blackoilmoduleparams.hh>
 #include <opm/models/blackoil/blackoilconvectivemixingmodule.hh>
 
-#include <opm/common/utility/gpuDecorators.hpp>
 
 #include <array>
 
@@ -227,12 +227,12 @@ class NewTranExtensiveQuantities
 public:
 
     OPM_HOST_DEVICE static void volumeAndPhasePressureDifferences(std::array<short, numPhases>& upIdx,
-                                                  std::array<short, numPhases>& dnIdx,
-                                                  Evaluation (&volumeFlux)[numPhases],
-                                                  Evaluation (&pressureDifferences)[numPhases],
-                                                  const ElementContext& elemCtx,
-                                                  unsigned scvfIdx,
-                                                  unsigned timeIdx)
+                                                                  std::array<short, numPhases>& dnIdx,
+                                                                  Evaluation (&volumeFlux)[numPhases],
+                                                                  Evaluation (&pressureDifferences)[numPhases],
+                                                                  const ElementContext& elemCtx,
+                                                                  unsigned scvfIdx,
+                                                                  unsigned timeIdx)
     {
         const auto& problem = elemCtx.problem();
         const auto& stencil = elemCtx.stencil(timeIdx);
@@ -312,23 +312,23 @@ class NewTranExtensiveQuantities
         }
     }
 
-    template<class EvalType>
+    template<class EvalType, class ModuleParamsT = ModuleParams>
     OPM_HOST_DEVICE static void calculatePhasePressureDiff_(short& upIdx,
-                                            short& dnIdx,
-                                            EvalType& pressureDifference,
-                                            const IntensiveQuantities& intQuantsIn,
-                                            const IntensiveQuantities& intQuantsEx,
-                                            const unsigned phaseIdx,
-                                            const unsigned interiorDofIdx,
-                                            const unsigned exteriorDofIdx,
-                                            const Scalar Vin,
-                                            const Scalar Vex,
-                                            const unsigned globalIndexIn,
-                                            const unsigned globalIndexEx,
-                                            const Scalar distZg,
-                                            const Scalar thpresInToEx,
-                                            const Scalar thpresExToIn,
-                                            const ModuleParams& moduleParams)
+                                                            short& dnIdx,
+                                                            EvalType& pressureDifference,
+                                                            const IntensiveQuantities& intQuantsIn,
+                                                            const IntensiveQuantities& intQuantsEx,
+                                                            const unsigned phaseIdx,
+                                                            const unsigned interiorDofIdx,
+                                                            const unsigned exteriorDofIdx,
+                                                            const Scalar Vin,
+                                                            const Scalar Vex,
+                                                            const unsigned globalIndexIn,
+                                                            const unsigned globalIndexEx,
+                                                            const Scalar distZg,
+                                                            const Scalar thpresInToEx,
+                                                            const Scalar thpresExToIn,
+                                                            const ModuleParamsT& moduleParams)
     {
 
         // check shortcut: if the mobility of the phase is zero in the interior as
@@ -434,9 +434,9 @@ class NewTranExtensiveQuantities
      */
     template <class FluidState>
     OPM_HOST_DEVICE void calculateBoundaryGradients_(const ElementContext& elemCtx,
-                                     unsigned scvfIdx,
-                                     unsigned timeIdx,
-                                     const FluidState& exFluidState)
+                                                     unsigned scvfIdx,
+                                                     unsigned timeIdx,
+                                                     const FluidState& exFluidState)
     {
         const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(scvfIdx);
         const Scalar faceArea = scvf.area();
@@ -478,16 +478,16 @@ class NewTranExtensiveQuantities
      */
     template <class Problem, class FluidState, class EvaluationContainer>
     OPM_HOST_DEVICE static void calculateBoundaryGradients_(const Problem& problem,
-                                            const unsigned globalSpaceIdx,
-                                            const IntensiveQuantities& intQuantsIn,
-                                            const unsigned bfIdx,
-                                            const double faceArea,
-                                            const double zEx,
-                                            const FluidState& exFluidState,
-                                            std::array<short, numPhases>& upIdx,
-                                            std::array<short, numPhases>& dnIdx,
-                                            EvaluationContainer& volumeFlux,
-                                            EvaluationContainer& pressureDifference)
+                                                            const unsigned globalSpaceIdx,
+                                                            const IntensiveQuantities& intQuantsIn,
+                                                            const unsigned bfIdx,
+                                                            const double faceArea,
+                                                            const double zEx,
+                                                            const FluidState& exFluidState,
+                                                            std::array<short, numPhases>& upIdx,
+                                                            std::array<short, numPhases>& dnIdx,
+                                                            EvaluationContainer& volumeFlux,
+                                                            EvaluationContainer& pressureDifference)
     {
 
         bool enableBoundaryMassFlux = problem.nonTrivialBoundaryConditions();
diff --git a/opm/simulators/flow/SimpleFIBlackOilModel.hpp b/opm/simulators/flow/SimpleFIBlackOilModel.hpp
new file mode 100644
index 00000000000..74aac2a2669
--- /dev/null
+++ b/opm/simulators/flow/SimpleFIBlackOilModel.hpp
@@ -0,0 +1,203 @@
+// -*- 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 <http://www.gnu.org/licenses/>.
+
+  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
+ *
+ * \copydoc Opm::FIBlackOilModel
+ */
+#ifndef SIMPLE_FI_BLACK_OIL_MODEL_HPP
+#define SIMPLE_FI_BLACK_OIL_MODEL_HPP
+
+#include <opm/models/blackoil/blackoilconvectivemixingmodule.hh>
+#include <opm/models/blackoil/blackoilmodel.hh>
+#include <opm/models/blackoil/blackoilmoduleparams.hh>
+#include <opm/models/parallel/threadmanager.hpp>
+#include <opm/models/utils/propertysystem.hh>
+
+#include <opm/common/ErrorMacros.hpp>
+
+#include <opm/grid/CpGrid.hpp>
+#include <opm/grid/utility/ElementChunks.hpp>
+
+#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
+
+#include <opm/material/fluidmatrixinteractions/EclMultiplexerMaterialParams.hpp>
+
+#include <cassert>
+#include <cstddef>
+#include <stdexcept>
+#include <type_traits>
+#include <vector>
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+/*
+    This file implements a simplified version of the FIBlackOilModel
+    which can be used in a GPU environment.
+    This is developed while parallelizing the tpfa matrix assembly, so
+    only functionality needed there is implemented.
+    To start with this means just being able to access intensiveQuantities.
+*/
+
+namespace Opm
+{
+
+template <typename TypeTag, template <class> class Storage = Opm::VectorWithDefaultAllocator>
+class SimpleFIBlackOilModel
+{
+public:
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using TypeTagPublic = TypeTag;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+
+    SimpleFIBlackOilModel(
+        const Storage<BlackOilIntensiveQuantities<TypeTag>>& cachedIntensiveQuantities0,
+        const Storage<BlackOilIntensiveQuantities<TypeTag>>& cachedIntensiveQuantities1,
+        const Storage<Scalar>& volumes)
+        : cachedIntensiveQuantities0_(cachedIntensiveQuantities0)
+        , cachedIntensiveQuantities1_(cachedIntensiveQuantities1)
+        , volumes_(volumes)
+    {
+    }
+
+    // Special ctor because fvbasediscretization stores BOIQs with a custom allocator.
+    // OtherIQContainer can be any container (e.g. std::vector with a non-default allocator).
+    // Could this be using refs instead of making the copy? Then the copy-to-gpu would have to be
+    // rewritten a bit
+    template <typename OtherIQContainer>
+    SimpleFIBlackOilModel(const OtherIQContainer& cachedIntensiveQuantities0,
+                          const OtherIQContainer& cachedIntensiveQuantities1,
+                          const Storage<Scalar>& volumes)
+        : cachedIntensiveQuantities0_(cachedIntensiveQuantities0.begin(),
+                                      cachedIntensiveQuantities0.end())
+        , cachedIntensiveQuantities1_(cachedIntensiveQuantities1.begin(),
+                                      cachedIntensiveQuantities1.end())
+        , volumes_(volumes)
+    {
+    }
+
+    SimpleFIBlackOilModel() = delete;
+
+    OPM_HOST_DEVICE Scalar dofTotalVolume(unsigned int globalIdx) const
+    {
+        return volumes_[globalIdx];
+    }
+
+    OPM_HOST_DEVICE const auto& intensiveQuantities(unsigned int globalIdx,
+                                                    unsigned int timeIdx) const
+    {
+        assert(timeIdx == 0
+               || timeIdx
+                   == 1); // TODO: do not use assert for this because it can be run in GPU kernel...
+        if (timeIdx == 0) {
+            return cachedIntensiveQuantities0_[globalIdx];
+        } else {
+            return cachedIntensiveQuantities1_[globalIdx];
+        }
+    }
+
+public:
+    // assumes --enable-storage-cache=false so that we need not worry about updating cpu-caches from
+    // gpu kernels
+    Storage<BlackOilIntensiveQuantities<TypeTag>> cachedIntensiveQuantities0_;
+    Storage<BlackOilIntensiveQuantities<TypeTag>> cachedIntensiveQuantities1_;
+    Storage<Scalar> volumes_;
+};
+
+#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
+namespace gpuistl
+{
+    // For now we make a copy of the simplified CPU model because we need to set the fluid system
+    // pointer without breaking copy semantics
+    template <class TypeTag, typename GpuFluidSystem>
+    auto copy_to_gpu(const SimpleFIBlackOilModel<TypeTag>& cpuModel, const GpuFluidSystem& fsys)
+    {
+        using Scalar = typename SimpleFIBlackOilModel<TypeTag>::Scalar;
+        // In this copy_to_gpu we need to take care of two things:
+        // 1) Copy the cachedIntensiveQuantities_ to GPU (go from vector to GpuBuffer to allocate
+        // things on GPU) 2) Ensure that the FluidSystem pointer inside each IntensiveQuantities
+        // points to a GPU FluidSystem
+        //    The pointer should be set before we move the entire thing to the GPU.
+
+        // set pointers
+        using CorrectTypeTagView =
+            typename ::Opm::Properties::TTag::to_gpu_type_t<TypeTag, GpuView>;
+        using CorrectBOIQ = BlackOilIntensiveQuantities<CorrectTypeTagView>;
+        using SimplifiedGpuBufferModel = SimpleFIBlackOilModel<CorrectTypeTagView, GpuBuffer>;
+
+        assert(cpuModel.cachedIntensiveQuantities0_.size() == cpuModel.cachedIntensiveQuantities1_.size());
+        const std::size_t nCells = cpuModel.cachedIntensiveQuantities0_.size();
+        assert(nCells > 0);
+        // CorrectBOIQ is not default-constructible (non-static FluidSystem requires a pointer at
+        // construction time). To enable parallel fill, we create a prototype from the first element
+        // and use it to pre-size the vector with valid objects, then overwrite in parallel.
+        // Assumption: all cells share the same fsys pointer, so the prototype object is a valid
+        // stand-in until it is overwritten.
+        std::vector<CorrectBOIQ> cpuIntQuantsWithGpuPtr0;
+        std::vector<CorrectBOIQ> cpuIntQuantsWithGpuPtr1;
+        CorrectBOIQ proto0 = cpuModel.cachedIntensiveQuantities0_[0]
+                                    .template withOtherFluidSystem<CorrectTypeTagView>(fsys);
+        CorrectBOIQ proto1 = cpuModel.cachedIntensiveQuantities1_[0]
+                                    .template withOtherFluidSystem<CorrectTypeTagView>(fsys);
+        cpuIntQuantsWithGpuPtr0.resize(nCells, proto0);
+        cpuIntQuantsWithGpuPtr1.resize(nCells, proto1);
+
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static)
+#endif
+        for (std::size_t i = 1; i < nCells; ++i) {
+            cpuIntQuantsWithGpuPtr0[i]
+                = cpuModel.cachedIntensiveQuantities0_[i]
+                        .template withOtherFluidSystem<CorrectTypeTagView>(fsys);
+            cpuIntQuantsWithGpuPtr1[i]
+                = cpuModel.cachedIntensiveQuantities1_[i]
+                        .template withOtherFluidSystem<CorrectTypeTagView>(fsys);
+        }
+
+        GpuBuffer<CorrectBOIQ> gpuBuf0(cpuIntQuantsWithGpuPtr0);
+        GpuBuffer<CorrectBOIQ> gpuBuf1(cpuIntQuantsWithGpuPtr1);
+        GpuBuffer<Scalar> gpuVolumes(cpuModel.volumes_);
+
+        // create new blackoil model providing a gpubuffer allocated version of the intQuants
+        auto result = SimplifiedGpuBufferModel(
+            std::move(gpuBuf0), std::move(gpuBuf1), std::move(gpuVolumes));
+
+        return result;
+    }
+
+    template <typename LocalTypeTag>
+    auto make_view(SimpleFIBlackOilModel<LocalTypeTag, GpuBuffer>& gpuSimpleFIBlackOilModel)
+    {
+        return SimpleFIBlackOilModel<LocalTypeTag, gpuistl::GpuView>(
+            make_view(gpuSimpleFIBlackOilModel.cachedIntensiveQuantities0_),
+            make_view(gpuSimpleFIBlackOilModel.cachedIntensiveQuantities1_),
+            make_view(gpuSimpleFIBlackOilModel.volumes_));
+    }
+} // namespace gpuistl
+#endif
+
+} // namespace Opm
+
+#endif // SIMPLE_FI_BLACK_OIL_MODEL_HPP
diff --git a/opm/simulators/flow/ThermalGasWaterFlowProblem.hpp b/opm/simulators/flow/ThermalGasWaterFlowProblem.hpp
new file mode 100644
index 00000000000..571b14b5478
--- /dev/null
+++ b/opm/simulators/flow/ThermalGasWaterFlowProblem.hpp
@@ -0,0 +1,125 @@
+/*
+  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 THERMAL_GASWATER_FLOW_PROBLEM_HPP
+#define THERMAL_GASWATER_FLOW_PROBLEM_HPP
+
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/common/utility/VectorWithDefaultAllocator.hpp>
+#include <opm/common/utility/gpuDecorators.hpp>
+#include <opm/common/utility/gpuistl_if_available.hpp>
+
+/*
+    A simple gas-water flow problem with thermal effects that is GPU compatible
+*/
+
+namespace Opm
+{
+
+template <class Scalar, template <class> class Storage = Opm::VectorWithDefaultAllocator>
+class ThermalGasWaterFlowProblem
+{
+public:
+    using ModuleParams = BlackoilModuleParams<ConvectiveMixingModuleParam<Scalar, Storage>>;
+
+    ThermalGasWaterFlowProblem() = default;
+
+    ThermalGasWaterFlowProblem(Storage<Scalar> alpha0,
+                               Storage<Scalar> alpha1,
+                               Storage<Scalar> alpha2,
+                               ModuleParams moduleParams)
+        : alpha0_(alpha0)
+        , alpha1_(alpha1)
+        , alpha2_(alpha2)
+        , moduleParams_(moduleParams)
+    {
+    }
+
+    OPM_HOST_DEVICE Scalar getAlpha(unsigned globalIndex, unsigned boundaryFaceIndex) const
+    {
+        if (boundaryFaceIndex == 0) {
+            return alpha0_[globalIndex];
+        } else if (boundaryFaceIndex == 1) {
+            return alpha1_[globalIndex];
+        } else if (boundaryFaceIndex == 2) {
+            return alpha2_[globalIndex];
+        } else {
+            OPM_THROW(std::logic_error, "Invalid boundary face index");
+        }
+    }
+
+    OPM_HOST_DEVICE const ModuleParams& moduleParams() const
+    {
+        return moduleParams_;
+    }
+    OPM_HOST_DEVICE ModuleParams& moduleParams()
+    {
+        return moduleParams_;
+    }
+
+    Storage<Scalar>& alpha0()
+    {
+        return alpha0_;
+    }
+
+    Storage<Scalar>& alpha1()
+    {
+        return alpha1_;
+    }
+
+    Storage<Scalar>& alpha2()
+    {
+        return alpha2_;
+    }
+
+private:
+    Storage<Scalar> alpha0_;
+    Storage<Scalar> alpha1_;
+    Storage<Scalar> alpha2_;
+    ModuleParams moduleParams_;
+};
+
+namespace gpuistl
+{
+
+    template <class Scalar>
+    ThermalGasWaterFlowProblem<Scalar, GpuBuffer>
+    copy_to_gpu(ThermalGasWaterFlowProblem<Scalar, Opm::VectorWithDefaultAllocator>& cpuProblem)
+    {
+        using ModuleParams = BlackoilModuleParams<ConvectiveMixingModuleParam<Scalar, GpuBuffer>>;
+        return ThermalGasWaterFlowProblem<Scalar, GpuBuffer>(
+            GpuBuffer<Scalar>(cpuProblem.alpha0()),
+            GpuBuffer<Scalar>(cpuProblem.alpha1()),
+            GpuBuffer<Scalar>(cpuProblem.alpha2()),
+            ModuleParams {copy_to_gpu(cpuProblem.moduleParams().convectiveMixingModuleParam)});
+    }
+
+    template <class Scalar>
+    ThermalGasWaterFlowProblem<Scalar, GpuView>
+    make_view(ThermalGasWaterFlowProblem<Scalar, GpuBuffer>& buffer)
+    {
+        using ModuleParams = BlackoilModuleParams<ConvectiveMixingModuleParam<Scalar, GpuView>>;
+        return ThermalGasWaterFlowProblem<Scalar, GpuView>(
+            make_view(buffer.alpha0()),
+            make_view(buffer.alpha1()),
+            make_view(buffer.alpha2()),
+            ModuleParams {make_view(buffer.moduleParams().convectiveMixingModuleParam)});
+    }
+
+} // namespace gpuistl
+
+} // namespace Opm
+
+#endif // THERMAL_GASWATER_FLOW_PROBLEM_HPP
diff --git a/opm/simulators/flow/Transmissibility.hpp b/opm/simulators/flow/Transmissibility.hpp
index e5e18c11a3e..9bc3d6c62bc 100644
--- a/opm/simulators/flow/Transmissibility.hpp
+++ b/opm/simulators/flow/Transmissibility.hpp
@@ -102,6 +102,8 @@ class Transmissibility {
 
     Scalar thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const;
 
+    const std::map<std::pair<unsigned, unsigned>, Scalar>& getThermalHalfTransBoundary() const;
+
     /*!
      * \brief Return the diffusivity for the intersection between two elements.
      */
diff --git a/opm/simulators/flow/Transmissibility_impl.hpp b/opm/simulators/flow/Transmissibility_impl.hpp
index 4b78798723c..c268d2eed91 100644
--- a/opm/simulators/flow/Transmissibility_impl.hpp
+++ b/opm/simulators/flow/Transmissibility_impl.hpp
@@ -145,6 +145,13 @@ thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const
     return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx));
 }
 
+template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
+const std::map<std::pair<unsigned, unsigned>, Scalar>& Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
+getThermalHalfTransBoundary() const
+{
+    return thermalHalfTransBoundary_;
+}
+
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
 Scalar Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
 diffusivity(unsigned elemIdx1, unsigned elemIdx2) const
diff --git a/opm/simulators/linalg/gpuistl/GpuBlackoilIntensiveQuantitiesDispatcher.cu b/opm/simulators/linalg/gpuistl/GpuBlackoilIntensiveQuantitiesDispatcher.cu
new file mode 100644
index 00000000000..0b676cae588
--- /dev/null
+++ b/opm/simulators/linalg/gpuistl/GpuBlackoilIntensiveQuantitiesDispatcher.cu
@@ -0,0 +1,299 @@
+/*
+  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/>.
+*/
+#include <config.h>
+
+// Tricks for older versions of ROCm that do not support alugrid or dune-fem
+// (mirrors test_blackoilintensivequantities_gpu.cu).
+#ifdef HAVE_DUNE_ALUGRID
+#undef HAVE_DUNE_ALUGRID
+#endif
+#define HAVE_DUNE_ALUGRID 0
+#ifdef HAVE_DUNE_FEM
+#undef HAVE_DUNE_FEM
+#endif
+#define HAVE_DUNE_FEM 0
+
+#include <cuda_runtime.h>
+
+#include <opm/common/utility/gpuDecorators.hpp>
+#include <opm/material/common/ResetLocale.hpp>
+#include <opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp>
+
+#include <opm/material/fluidmatrixinteractions/EclMaterialLawTwoPhaseTypes.hpp>
+#include <opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp>
+#include <opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp>
+#include <opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp>
+
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystemNonStatic.hpp>
+
+#include <opm/models/blackoil/blackoilintensivequantities.hh>
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
+#include <opm/models/blackoil/blackoilmodel.hh>
+#include <opm/models/blackoil/blackoilprimaryvariables.hh>
+#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
+#include <opm/models/discretization/common/tpfalinearizer.hh>
+
+#include <opm/simulators/flow/FlowProblemBlackoil.hpp>
+#include <opm/simulators/flow/FlowProblemBlackoilProperties.hpp>
+#include <opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp>
+#include <opm/material/thermal/GpuEclThermalLawManager.hpp>
+#include <opm/simulators/flow/GpuFlowProblem.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuBuffer.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuView.hpp>
+#include <opm/simulators/linalg/gpuistl/detail/gpu_safe_call.hpp>
+
+#include <opm/simulators/linalg/gpuistl/GpuBlackoilIntensiveQuantitiesDispatcher.hpp>
+
+#include <opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuFlowGasWaterEnergyTypeTags.hpp>
+
+#include <opm/common/OpmLog/OpmLog.hpp>
+
+#include <format>
+#include <optional>
+#include <vector>
+
+namespace Opm::gpuistl {
+
+namespace {
+
+using DispatcherCpuTag = Opm::Properties::TTag::FlowGasWaterEnergyCpuKernelBase;
+using DispatcherGpuTag = Opm::Properties::TTag::FlowGasWaterEnergyDummyProblemGPU;
+
+using DispatcherScalar = Opm::GetPropType<DispatcherCpuTag, Opm::Properties::Scalar>;
+
+using DispatcherCpuFluidSystem = Opm::BlackOilFluidSystem<DispatcherScalar>;
+using DispatcherFluidSystemView
+    = Opm::BlackOilFluidSystemNonStatic<DispatcherScalar,
+                                        Opm::BlackOilDefaultFluidSystemIndices,
+                                        Opm::gpuistl::GpuView>;
+
+using DispatcherGpuPrimaryVariables
+    = Opm::BlackOilPrimaryVariables<DispatcherGpuTag, Opm::gpuistl::MiniVector>;
+using DispatcherGpuIntensiveQuantities = Opm::BlackOilIntensiveQuantities<DispatcherGpuTag>;
+
+using DispatcherCpuMaterialLawManager =
+    typename Opm::GetProp<DispatcherCpuTag, Opm::Properties::MaterialLaw>::EclMaterialLawManager;
+using DispatcherTraits = typename DispatcherCpuMaterialLawManager::MaterialLaw::Traits;
+using DispatcherTwoPhaseTraits =
+    Opm::TwoPhaseMaterialTraits<DispatcherScalar,
+                                DispatcherTraits::wettingPhaseIdx,
+                                DispatcherTraits::nonWettingPhaseIdx>;
+using DispatcherGpuPiecewiseLinearParamsBuf =
+    Opm::PiecewiseLinearTwoPhaseMaterialParams<DispatcherTwoPhaseTraits,
+                                               Opm::gpuistl::GpuView<const DispatcherScalar>>;
+using DispatcherGpuPiecewiseLinearLawBuf =
+    Opm::PiecewiseLinearTwoPhaseMaterial<DispatcherTwoPhaseTraits,
+                                         DispatcherGpuPiecewiseLinearParamsBuf>;
+using DispatcherGpuMaterialLawParamsBuf =
+    Opm::EclTwoPhaseMaterialParams<DispatcherTraits,
+                                   DispatcherGpuPiecewiseLinearParamsBuf,
+                                   DispatcherGpuPiecewiseLinearParamsBuf,
+                                   DispatcherGpuPiecewiseLinearParamsBuf,
+                                   Opm::gpuistl::ValueAsPointer>;
+using DispatcherGpuMaterialLawBuf =
+    Opm::EclTwoPhaseMaterial<DispatcherTraits,
+                             DispatcherGpuPiecewiseLinearLawBuf,
+                             DispatcherGpuPiecewiseLinearLawBuf,
+                             DispatcherGpuPiecewiseLinearLawBuf,
+                             DispatcherGpuMaterialLawParamsBuf>;
+using DispatcherGpuManagerBuf =
+    Opm::EclMaterialLaw::GpuManager<DispatcherTraits,
+                                    DispatcherGpuPiecewiseLinearLawBuf,
+                                    DispatcherGpuPiecewiseLinearLawBuf,
+                                    Opm::gpuistl::GpuBuffer,
+                                    DispatcherGpuMaterialLawBuf>;
+using DispatcherGpuThermalManagerBuf =
+    Opm::EclThermalLaw::GpuManager<DispatcherScalar,
+                                   DispatcherFluidSystemView,
+                                   Opm::gpuistl::GpuBuffer,
+                                   Opm::gpuistl::GpuView>;
+using DispatcherGpuFlowProblemBuf =
+    Opm::GpuFlowProblem<DispatcherScalar,
+                        DispatcherGpuManagerBuf,
+                        Opm::gpuistl::GpuBuffer,
+                        DispatcherGpuThermalManagerBuf>;
+using DispatcherGpuFlowProblemView =
+    decltype(Opm::gpuistl::make_view(std::declval<DispatcherGpuFlowProblemBuf&>()));
+
+template <class GpuProblem, class PrimaryVariablesT, class IntensiveQuantitiesT>
+__global__ void
+dispatcherUpdateAllCellsKernel(GpuProblem problem,
+                               Opm::gpuistl::GpuView<const PrimaryVariablesT> primaryVariables,
+                               Opm::gpuistl::GpuView<IntensiveQuantitiesT> outIntensiveQuantities,
+                               std::size_t numCells)
+{
+    const std::size_t i = static_cast<std::size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
+    if (i >= numCells) {
+        return;
+    }
+    IntensiveQuantitiesT& iq = outIntensiveQuantities[i];
+    iq.updateSaturations(primaryVariables[i], 0, Opm::LinearizationType{});
+    iq.update(problem, primaryVariables[i], static_cast<unsigned>(i), 0);
+    iq.updateEnergyQuantities_(problem, static_cast<unsigned>(i), 0u);
+}
+
+} // namespace
+
+// =============================================================================
+// Impl
+// =============================================================================
+template <class CpuTypeTag>
+struct GpuBlackoilIntensiveQuantitiesDispatcher<CpuTypeTag>::Impl {
+    bool initialized = false;
+    std::unique_ptr<DispatcherGpuFlowProblemBuf> problemBuf;
+    DispatcherGpuFlowProblemView problemView{};
+    DispatcherFluidSystemView* managedFluidSystemView = nullptr;
+    std::optional<DispatcherGpuIntensiveQuantities> prototype;
+
+    ~Impl()
+    {
+        if (managedFluidSystemView != nullptr) {
+            managedFluidSystemView->~DispatcherFluidSystemView();
+            (void)cudaFree(managedFluidSystemView);
+            managedFluidSystemView = nullptr;
+        }
+    }
+};
+
+template <class CpuTypeTag>
+GpuBlackoilIntensiveQuantitiesDispatcher<CpuTypeTag>::GpuBlackoilIntensiveQuantitiesDispatcher()
+    : impl_(std::make_unique<Impl>())
+{
+}
+
+template <class CpuTypeTag>
+GpuBlackoilIntensiveQuantitiesDispatcher<CpuTypeTag>::~GpuBlackoilIntensiveQuantitiesDispatcher()
+    = default;
+
+template <class CpuTypeTag>
+void GpuBlackoilIntensiveQuantitiesDispatcher<CpuTypeTag>::update(
+    const Problem& cpuProblem,
+    const PrimaryVariables* const* cpuPriVars,
+    IntensiveQuantities* const* outIQ,
+    std::size_t numDof)
+{
+    if (numDof == 0u) {
+        return;
+    }
+
+    if (!impl_->initialized) {
+        // Build the GPU FlowProblem from the CPU FlowProblem (one-time setup).
+        impl_->problemBuf = std::make_unique<DispatcherGpuFlowProblemBuf>(cpuProblem);
+        impl_->problemView = Opm::gpuistl::make_view(*impl_->problemBuf);
+
+        // Place the FluidSystemView in unified memory so its device pointer
+        // dereferences are valid both on host and device (mirrors the test).
+        auto& dynamicCpuFluidSystem = DispatcherCpuFluidSystem::getNonStaticInstance();
+        // The fluid-system buffer is owning device memory; leak it intentionally
+        // for the lifetime of the process via a function-local static.
+        using FsBufferType = decltype(Opm::gpuistl::copy_to_gpu(dynamicCpuFluidSystem));
+        static FsBufferType s_fsBuffer = Opm::gpuistl::copy_to_gpu(dynamicCpuFluidSystem);
+        auto fsView = Opm::gpuistl::make_view(s_fsBuffer);
+
+        DispatcherFluidSystemView* managed = nullptr;
+        OPM_GPU_SAFE_CALL(cudaMallocManaged(&managed, sizeof(DispatcherFluidSystemView)));
+        new (managed) DispatcherFluidSystemView(fsView);
+        impl_->managedFluidSystemView = managed;
+
+        // Build a default IntensiveQuantities prototype for the GPU side.
+        Opm::BlackOilIntensiveQuantities<DispatcherCpuTag> cpuPrototype;
+        impl_->prototype = cpuPrototype.template withOtherFluidSystem<DispatcherGpuTag>(
+            *impl_->managedFluidSystemView);
+
+        impl_->initialized = true;
+        Opm::OpmLog::info(std::format(
+            "[GpuBlackoilIntensiveQuantitiesDispatcher] initialized for {} cells",
+            cpuProblem.model().numGridDof()));
+    }
+
+    // -------------------------------------------------------------------
+    // 1. Convert CPU primary variables to the dispatcher's GPU primary
+    //    variables (host-side).
+    // -------------------------------------------------------------------
+    std::vector<DispatcherGpuPrimaryVariables> hostPriVars;
+    hostPriVars.reserve(numDof);
+    for (std::size_t i = 0; i < numDof; ++i) {
+        // The BlackOilPrimaryVariables converting copy ctor handles the
+        // TypeTag/storage difference.
+        hostPriVars.emplace_back(*cpuPriVars[i]);
+    }
+    std::vector<DispatcherGpuIntensiveQuantities> hostIQ(numDof, *impl_->prototype);
+
+    // -------------------------------------------------------------------
+    // 2. Upload host buffers to the device (host-to-device copies are
+    //    embedded in the GpuBuffer ctors).
+    // -------------------------------------------------------------------
+    Opm::gpuistl::GpuBuffer<DispatcherGpuPrimaryVariables> primaryVariablesBuffer(hostPriVars);
+    Opm::gpuistl::GpuBuffer<DispatcherGpuIntensiveQuantities> intensiveQuantitiesBuffer(hostIQ);
+
+    // -------------------------------------------------------------------
+    // 3. Launch the per-cell update kernel.
+    // -------------------------------------------------------------------
+    const unsigned blockSize = 64u;
+    const unsigned gridSize = static_cast<unsigned>((numDof + blockSize - 1u) / blockSize);
+
+    dispatcherUpdateAllCellsKernel<<<gridSize, blockSize>>>(
+        impl_->problemView,
+        Opm::gpuistl::GpuView<const DispatcherGpuPrimaryVariables>(
+            primaryVariablesBuffer.data(), primaryVariablesBuffer.size()),
+        Opm::gpuistl::GpuView<DispatcherGpuIntensiveQuantities>(
+            intensiveQuantitiesBuffer.data(), intensiveQuantitiesBuffer.size()),
+        numDof);
+
+    // -------------------------------------------------------------------
+    // 4. Read the GPU IntensiveQuantities back to host memory.
+    // -------------------------------------------------------------------
+    OPM_GPU_SAFE_CALL(cudaMemcpy(hostIQ.data(),
+                                 intensiveQuantitiesBuffer.data(),
+                                 numDof * sizeof(DispatcherGpuIntensiveQuantities),
+                                 cudaMemcpyDeviceToHost));
+
+    // -------------------------------------------------------------------
+    // 5. Field-by-field overlay onto the caller's CPU IntensiveQuantities.
+    //    The caller is expected to have run the CPU update first to
+    //    populate any fields the dispatcher does not overwrite (e.g.
+    //    \c mobility_, which the GPU relperm path currently leaves at 0).
+    // -------------------------------------------------------------------
+    for (std::size_t i = 0; i < numDof; ++i) {
+        outIQ[i]->overlayBlackOilFieldsFrom(hostIQ[i]);
+    }
+}
+
+// =============================================================================
+// Explicit instantiation for the user-facing CO2STORE simulation TypeTag,
+// namely \c FlowGasWaterEnergyProblem (used by the Flow simulator binary
+// when a CO2STORE deck with WATER+GAS+THERMAL is loaded). The GPU-specific
+// property bindings live in \c GpuFlowGasWaterEnergyTypeTags.hpp.
+// =============================================================================
+template class GpuBlackoilIntensiveQuantitiesDispatcher<
+    Opm::Properties::TTag::FlowGasWaterEnergyProblem>;
+
+// =============================================================================
+// Explicit instantiation for the GPU-assembly simulation TypeTag
+// \c FlowGasWaterEnergyProblemGPU (used by the flow_gpu binary that runs both
+// matrix assembly and intensive-quantities computation on the GPU).
+// This TypeTag inherits all physics from \c FlowGasWaterEnergyProblem and
+// therefore has the same CPU-side Problem / PrimaryVariables /
+// IntensiveQuantities types; the GPU kernel dispatch internally still uses
+// the \c FlowGasWaterEnergyKernelBaseGPU / \c FlowGasWaterEnergyDummyProblemGPU
+// chain.
+// =============================================================================
+template class GpuBlackoilIntensiveQuantitiesDispatcher<
+    Opm::Properties::TTag::FlowGasWaterEnergyProblemGPU>;
+
+} // namespace Opm::gpuistl
diff --git a/opm/simulators/linalg/gpuistl/GpuBlackoilIntensiveQuantitiesDispatcher.hpp b/opm/simulators/linalg/gpuistl/GpuBlackoilIntensiveQuantitiesDispatcher.hpp
new file mode 100644
index 00000000000..ef75fc14362
--- /dev/null
+++ b/opm/simulators/linalg/gpuistl/GpuBlackoilIntensiveQuantitiesDispatcher.hpp
@@ -0,0 +1,113 @@
+// -*- 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_BLACKOIL_INTENSIVE_QUANTITIES_DISPATCHER_HPP
+#define OPM_GPU_BLACKOIL_INTENSIVE_QUANTITIES_DISPATCHER_HPP
+
+#if HAVE_CUDA
+
+#include <opm/models/utils/propertysystem.hh>
+
+#include <cstddef>
+#include <memory>
+
+namespace Opm::Properties::TTag {
+    struct FlowGasWaterEnergyProblem;
+    struct FlowGasWaterEnergyProblemGPU;
+}
+
+namespace Opm::gpuistl {
+
+/// Compile-time predicate: does this CPU \c TypeTag describe the specific
+/// CO2STORE configuration (FlowGasWaterEnergyProblem) that the GPU
+/// intensive-quantities dispatcher currently supports? Restricted to a
+/// single TypeTag so that explicit instantiations of the dispatcher
+/// remain manageable for the experimental prototype.
+template <class CpuTypeTag>
+struct GpuBlackoilIntensiveQuantitiesDispatcherSupport {
+    static constexpr bool value = false;
+};
+
+template <>
+struct GpuBlackoilIntensiveQuantitiesDispatcherSupport<
+    Opm::Properties::TTag::FlowGasWaterEnergyProblem> {
+    static constexpr bool value = true;
+};
+
+/// Enable the dispatcher for the GPU-assembly simulation TypeTag
+/// (\c FlowGasWaterEnergyProblemGPU, declared in FlowGasWaterEnergyTypeTag.hpp).
+/// This tag inherits all physics from \c FlowGasWaterEnergyProblem and adds
+/// GPU-specific assembly properties in \c flow_gpu.cu.
+template <>
+struct GpuBlackoilIntensiveQuantitiesDispatcherSupport<
+    Opm::Properties::TTag::FlowGasWaterEnergyProblemGPU> {
+    static constexpr bool value = true;
+};
+
+/// Runs the BlackOil intensive-quantities update on the GPU for a given set
+/// of degrees of freedom. The dispatcher is process-wide (one singleton per
+/// CPU \c TypeTag). It lazily constructs the GPU-side \c GpuFlowProblem
+/// from the supplied CPU problem on the first call.
+///
+/// On every call, primary variables for the requested DoFs are uploaded to
+/// the device, the per-cell update kernel from
+/// \c test_blackoilintensivequantities_gpu.cu is launched (one thread per
+/// DoF), and the resulting intensive quantities are read back to host
+/// memory, as the design plan in \c PLAN.md requires.
+///
+/// The class is a template on the CPU \c TypeTag and is explicitly
+/// instantiated in the \c .cu translation unit.
+template <class CpuTypeTag>
+class GpuBlackoilIntensiveQuantitiesDispatcher
+{
+public:
+    using Problem            = Opm::GetPropType<CpuTypeTag, Opm::Properties::Problem>;
+    using PrimaryVariables   = Opm::GetPropType<CpuTypeTag, Opm::Properties::PrimaryVariables>;
+    using IntensiveQuantities = Opm::GetPropType<CpuTypeTag, Opm::Properties::IntensiveQuantities>;
+
+    GpuBlackoilIntensiveQuantitiesDispatcher();
+    ~GpuBlackoilIntensiveQuantitiesDispatcher();
+
+    GpuBlackoilIntensiveQuantitiesDispatcher(const GpuBlackoilIntensiveQuantitiesDispatcher&) = delete;
+    GpuBlackoilIntensiveQuantitiesDispatcher&
+    operator=(const GpuBlackoilIntensiveQuantitiesDispatcher&) = delete;
+
+    /// Run the per-cell intensive-quantities update kernel on \p numDof
+    /// DoFs. \p cpuPriVars[i] points at the CPU primary variables for DoF
+    /// \c i. The GPU-computed BlackOil intensive quantities are written
+    /// onto \p outIQ[i] field-by-field via
+    /// \c BlackOilIntensiveQuantities::overlayBlackOilFieldsFrom (so the
+    /// caller is expected to have run the CPU update first to fill in any
+    /// fields the dispatcher does not overwrite, e.g. \c mobility_).
+    void update(const Problem& cpuProblem,
+                const PrimaryVariables* const* cpuPriVars,
+                IntensiveQuantities* const* outIQ,
+                std::size_t numDof);
+
+private:
+    struct Impl;
+    std::unique_ptr<Impl> impl_;
+};
+
+} // namespace Opm::gpuistl
+
+#endif // HAVE_CUDA
+
+#endif // OPM_GPU_BLACKOIL_INTENSIVE_QUANTITIES_DISPATCHER_HPP
diff --git a/opm/simulators/linalg/gpuistl/GpuBuffer.hpp b/opm/simulators/linalg/gpuistl/GpuBuffer.hpp
index b89fc0fc3f0..11d48b301ff 100644
--- a/opm/simulators/linalg/gpuistl/GpuBuffer.hpp
+++ b/opm/simulators/linalg/gpuistl/GpuBuffer.hpp
@@ -87,6 +87,42 @@ class GpuBuffer
                                     cudaMemcpyDeviceToDevice));
     }
 
+    /**
+     * @brief Move constructor: transfers ownership of GPU memory without copying.
+     *
+     * The source buffer is left in a valid, empty state (size 0, null pointer) so
+     * that its destructor's cudaFree(nullptr) is a harmless no-op.
+     *
+     * @param other the GpuBuffer to move from
+     */
+    GpuBuffer(GpuBuffer<T>&& other) noexcept
+        : m_dataOnDevice(other.m_dataOnDevice)
+        , m_numberOfElements(other.m_numberOfElements)
+    {
+        other.m_dataOnDevice = nullptr;
+        other.m_numberOfElements = 0;
+    }
+
+    /**
+     * @brief Move assignment: frees current GPU memory, then takes ownership from other.
+     *
+     * The source buffer is left in a valid, empty state after the transfer.
+     *
+     * @param other the GpuBuffer to move from
+     * @return *this
+     */
+    GpuBuffer<T>& operator=(GpuBuffer<T>&& other) noexcept
+    {
+        if (this != &other) {
+            OPM_GPU_WARN_IF_ERROR(cudaFree(m_dataOnDevice));
+            m_dataOnDevice = other.m_dataOnDevice;
+            m_numberOfElements = other.m_numberOfElements;
+            other.m_dataOnDevice = nullptr;
+            other.m_numberOfElements = 0;
+        }
+        return *this;
+    }
+
     /**
      * @brief GpuBuffer allocates new GPU memory of the same size as data and copies the content of the data vector to
      * this newly allocated memory.
diff --git a/opm/simulators/linalg/gpuistl/GpuFlowGasWaterEnergyTypeTags.hpp b/opm/simulators/linalg/gpuistl/GpuFlowGasWaterEnergyTypeTags.hpp
new file mode 100644
index 00000000000..6daaa4a0599
--- /dev/null
+++ b/opm/simulators/linalg/gpuistl/GpuFlowGasWaterEnergyTypeTags.hpp
@@ -0,0 +1,270 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  Copyright TODO ADD YEAR AND NAME OF AUTHOR
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+/*!
+ * \file
+ *
+ * Shared GPU TypeTag derivations on top of the simulation-side
+ * \c FlowGasWaterEnergyProblem TypeTag (CO2STORE configuration).
+ *
+ * The intent is that GPU code (the dispatcher and the
+ * \c test_blackoilintensivequantities_gpu test) reuses the *exact* same
+ * TypeTag as the simulator and only overrides the small set of properties
+ * that have to change for the GPU build (storage container of the primary
+ * variables, the material law manager, the fluid system, the problem and
+ * the element context).
+ */
+#ifndef OPM_GPU_FLOW_GASWATER_ENERGY_TYPETAGS_HPP
+#define OPM_GPU_FLOW_GASWATER_ENERGY_TYPETAGS_HPP
+
+#if HAVE_CUDA
+
+#include <opm/material/fluidmatrixinteractions/EclMaterialLawTwoPhaseTypes.hpp>
+#include <opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp>
+#include <opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp>
+#include <opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp>
+
+#include <opm/material/fluidsystems/BlackOilFluidSystemNonStatic.hpp>
+
+#include <opm/models/blackoil/blackoilprimaryvariables.hh>
+
+#include <opm/material/thermal/EclSpecrockLaw.hpp>
+#include <opm/material/thermal/EclSpecrockLawParams.hpp>
+#include <opm/material/thermal/EclThconrLaw.hpp>
+#include <opm/material/thermal/EclThconrLawParams.hpp>
+
+#include <opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp>
+#include <opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp>
+#include <opm/material/thermal/GpuEclThermalLawManager.hpp>
+#include <opm/simulators/flow/GpuFlowProblem.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuBuffer.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuView.hpp>
+#include <opm/simulators/linalg/gpuistl/MiniVector.hpp>
+#include <opm/simulators/linalg/gpuistl/gpu_smart_pointer.hpp>
+
+#include <tuple>
+
+namespace Opm::Properties {
+
+namespace TTag {
+
+/// Internal kernel-base TypeTag for the GPU intensive-quantities dispatcher.
+/// Inherits all property bindings from \c FlowGasWaterEnergyProblem and
+/// overrides the primary variables, material law, fluid system and thermal
+/// laws to their GPU-compatible (GpuView / MiniVector) equivalents.
+/// This tag is an implementation detail of the dispatcher; user code should
+/// use \c FlowGasWaterEnergyProblemGPU (declared in FlowGasWaterEnergyTypeTag.hpp)
+/// as the outer simulation TypeTag.
+struct FlowGasWaterEnergyKernelBaseGPU {
+    using InheritsFrom = std::tuple<FlowGasWaterEnergyProblem>;
+};
+
+/// A second derived TypeTag used for the per-cell intensive-quantities
+/// kernel: it additionally swaps in \c GpuFlowProblem so that the GPU device side has trivially
+/// copyable \c Problem / \c ElementContext types.
+struct FlowGasWaterEnergyDummyProblemGPU {
+    using InheritsFrom = std::tuple<FlowGasWaterEnergyKernelBaseGPU>;
+};
+
+/// CPU-side TypeTag used as the source when constructing a GPU
+/// intensive-quantities prototype via the cross-TypeTag converting constructor.
+/// Inherits all CPU-side property bindings from \c FlowGasWaterEnergyProblem
+/// but disables diffusion and dispersion to match the GPU kernel's capabilities,
+/// ensuring the converting constructor's static_assert on matching settings is
+/// satisfied.
+struct FlowGasWaterEnergyCpuKernelBase {
+    using InheritsFrom = std::tuple<FlowGasWaterEnergyProblem>;
+};
+
+} // namespace TTag
+
+// -----------------------------------------------------------------------
+// Diffusion / dispersion: not yet implemented in the GPU IntensiveQuantities
+// update path; force them off on the GPU TypeTag.
+// -----------------------------------------------------------------------
+template <class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::FlowGasWaterEnergyKernelBaseGPU>
+{
+    static constexpr bool value = false;
+};
+
+template <class TypeTag>
+struct EnableDispersion<TypeTag, TTag::FlowGasWaterEnergyKernelBaseGPU>
+{
+    static constexpr bool value = false;
+};
+
+// Also disable diffusion and dispersion for the CPU prototype TypeTag so that
+// the BlackOilIntensiveQuantities converting constructor's static_assert
+// (which requires matching diffusion/dispersion settings) is satisfied when
+// building the GPU IQ prototype from a CPU IQ.
+template <class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::FlowGasWaterEnergyCpuKernelBase>
+{
+    static constexpr bool value = false;
+};
+
+template <class TypeTag>
+struct EnableDispersion<TypeTag, TTag::FlowGasWaterEnergyCpuKernelBase>
+{
+    static constexpr bool value = false;
+};
+
+// -----------------------------------------------------------------------
+// MaterialLaw: GPU-friendly EclTwoPhaseMaterial backed by GpuView storage.
+// -----------------------------------------------------------------------
+template <class TypeTag>
+struct MaterialLaw<TypeTag, TTag::FlowGasWaterEnergyKernelBaseGPU>
+{
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TTag::FlowProblem, Properties::FluidSystem>;
+    using Traits = ThreePhaseMaterialTraits<Scalar,
+                                            /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
+                                            /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
+                                            /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx,
+                                            getPropValue<TypeTag, Properties::EnableHysteresis>(),
+                                            getPropValue<TypeTag, Properties::EnableEndpointScaling>()>;
+    using TwoPhaseTraits = TwoPhaseMaterialTraits<Scalar,
+                                                  /*wettingPhaseIdx=*/Traits::wettingPhaseIdx,
+                                                  /*nonWettingPhaseIdx=*/Traits::nonWettingPhaseIdx>;
+    using TwoPhaseParams = ::Opm::PiecewiseLinearTwoPhaseMaterialParams<
+        TwoPhaseTraits, ::Opm::gpuistl::GpuView<const Scalar>>;
+    using TwoPhaseLaw = ::Opm::PiecewiseLinearTwoPhaseMaterial<TwoPhaseTraits, TwoPhaseParams>;
+    using GpuMaterialLawParams = ::Opm::EclTwoPhaseMaterialParams<
+        Traits, TwoPhaseParams, TwoPhaseParams, TwoPhaseParams, ::Opm::gpuistl::ValueAsPointer>;
+    using GpuMaterialLaw = ::Opm::EclTwoPhaseMaterial<
+        Traits, TwoPhaseLaw, TwoPhaseLaw, TwoPhaseLaw, GpuMaterialLawParams>;
+
+public:
+    using EclMaterialLawManager = ::Opm::EclMaterialLaw::GpuManager<
+        Traits, TwoPhaseLaw, TwoPhaseLaw, ::Opm::VectorWithDefaultAllocator, GpuMaterialLaw>;
+    using type = typename EclMaterialLawManager::MaterialLaw;
+};
+
+// -----------------------------------------------------------------------
+// PrimaryVariables: GPU build uses MiniVector storage so the type is
+// trivially copyable to the device.
+// -----------------------------------------------------------------------
+template <class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::FlowGasWaterEnergyKernelBaseGPU>
+{
+    using type = ::Opm::BlackOilPrimaryVariables<TypeTag, ::Opm::gpuistl::MiniVector>;
+};
+
+// -----------------------------------------------------------------------
+// FluidSystem: non-static, GpuView-backed instantiation.
+// -----------------------------------------------------------------------
+template <class TypeTag>
+struct FluidSystem<TypeTag, TTag::FlowGasWaterEnergyKernelBaseGPU>
+{
+    using type = ::Opm::BlackOilFluidSystemNonStatic<
+        GetPropType<TTag::FlowProblem, Properties::Scalar>,
+        ::Opm::BlackOilDefaultFluidSystemIndices,
+        ::Opm::gpuistl::GpuView>;
+};
+
+template <class TypeTag>
+struct FluidSystem<TypeTag, TTag::FlowGasWaterEnergyDummyProblemGPU>
+{
+    using type = ::Opm::BlackOilFluidSystemNonStatic<
+        GetPropType<TTag::FlowProblem, Properties::Scalar>,
+        ::Opm::BlackOilDefaultFluidSystemIndices,
+        ::Opm::gpuistl::GpuView>;
+};
+
+// -----------------------------------------------------------------------
+// SolidEnergyLaw / ThermalConductionLaw: bind to the GPU-portable laws
+// (EclSpecrockLaw / EclThconrLaw). Both expose an
+// \c EclThermalLawManager nested typedef matching the convention used by
+// FlowBaseProblemProperties so existing GetProp<...>::EclThermalLawManager
+// chains keep working.
+// -----------------------------------------------------------------------
+template <class TypeTag>
+struct SolidEnergyLaw<TypeTag, TTag::FlowGasWaterEnergyKernelBaseGPU>
+{
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+public:
+    using EclThermalLawManager = ::Opm::EclThermalLaw::GpuManager<
+        Scalar, FluidSystem, ::Opm::gpuistl::GpuView, ::Opm::gpuistl::GpuView>;
+
+    using type = ::Opm::EclSpecrockLaw<
+        Scalar,
+        ::Opm::EclSpecrockLawParams<Scalar, ::Opm::gpuistl::GpuView>>;
+};
+
+template <class TypeTag>
+struct ThermalConductionLaw<TypeTag, TTag::FlowGasWaterEnergyKernelBaseGPU>
+{
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+public:
+    using EclThermalLawManager = ::Opm::EclThermalLaw::GpuManager<
+        Scalar, FluidSystem, ::Opm::gpuistl::GpuView, ::Opm::gpuistl::GpuView>;
+
+    using type = ::Opm::EclThconrLaw<Scalar, FluidSystem>;
+};
+
+// -----------------------------------------------------------------------
+// Problem: GpuFlowProblem (GpuView-backed) instead of FlowProblemBlackoil.
+// -----------------------------------------------------------------------
+template <class TypeTag>
+struct Problem<TypeTag, TTag::FlowGasWaterEnergyDummyProblemGPU>
+{
+private:
+    using ScalarT = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using CpuMaterialLawManager =
+        typename Opm::GetProp<TypeTag, Opm::Properties::MaterialLaw>::EclMaterialLawManager;
+    using GpuViewMaterialLawManager =
+        ::Opm::EclMaterialLaw::GpuManager<typename CpuMaterialLawManager::Traits,
+                                          typename CpuMaterialLawManager::GasOilLaw,
+                                          typename CpuMaterialLawManager::OilWaterLaw,
+                                          ::Opm::gpuistl::GpuView,
+                                          typename CpuMaterialLawManager::MaterialLaw>;
+    using GpuViewThermalLawManager = ::Opm::EclThermalLaw::GpuManager<
+        ScalarT, FluidSystem, ::Opm::gpuistl::GpuView, ::Opm::gpuistl::GpuView>;
+
+public:
+    using type = ::Opm::GpuFlowProblem<ScalarT,
+                                       GpuViewMaterialLawManager,
+                                       ::Opm::gpuistl::GpuView,
+                                       GpuViewThermalLawManager>;
+};
+
+// -----------------------------------------------------------------------
+// ElementContext: trivially copyable GPU element context.
+// -----------------------------------------------------------------------
+template <class TypeTag>
+struct ElementContext<TypeTag, TTag::FlowGasWaterEnergyDummyProblemGPU>
+{
+    using type = std::nullptr_t;
+};
+
+} // namespace Opm::Properties
+
+#endif // HAVE_CUDA
+
+#endif // OPM_GPU_FLOW_GASWATER_ENERGY_TYPETAGS_HPP
diff --git a/opm/simulators/linalg/gpuistl/GpuView.hpp b/opm/simulators/linalg/gpuistl/GpuView.hpp
index 79940acb5a1..9de0152fc9f 100644
--- a/opm/simulators/linalg/gpuistl/GpuView.hpp
+++ b/opm/simulators/linalg/gpuistl/GpuView.hpp
@@ -80,7 +80,7 @@ class GpuView
      *
      * @param idx The index of the element
      */
-    __host__ __device__ T operator[](size_t idx) const {
+    __host__ __device__ const T& operator[](size_t idx) const {
 #ifndef NDEBUG
         assertInRange(idx);
 #endif
diff --git a/opm/simulators/linalg/gpuistl/MiniVector.hpp b/opm/simulators/linalg/gpuistl/MiniVector.hpp
index 921f3a9271a..a23d474edd6 100644
--- a/opm/simulators/linalg/gpuistl/MiniVector.hpp
+++ b/opm/simulators/linalg/gpuistl/MiniVector.hpp
@@ -181,6 +181,12 @@ class MiniVector
         return Dimension;
     }
 
+    /** @return Always \c false; a fixed‑size vector is never empty. */
+    OPM_HOST_DEVICE static constexpr bool empty() noexcept
+    {
+        return Dimension == 0;
+    }
+
 
     /**
      * @brief Fill every component with the supplied value.
diff --git a/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.cu b/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.cu
index ce20dbbcd74..3fdcc2fc818 100644
--- a/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.cu
+++ b/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.cu
@@ -21,6 +21,7 @@
 #include <opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp>
 #include <opm/simulators/linalg/gpuistl/detail/deviceBlockOperations.hpp>
 #include <opm/simulators/linalg/gpuistl/detail/gpuThreadUtils.hpp>
+#include <opm/simulators/linalg/gpuistl/MiniMatrix.hpp>
 
 namespace Opm::gpuistl::detail
 {
@@ -97,13 +98,13 @@ namespace
         }
     }
 
-    template <class T>
-    __global__ void cuComputeDiagPtrs(const int* rowIndices,
-                                      const int* colIndices,
-                                      int rows,
-                                      T** diagPtrs,
-                                      T* values,
-                                      size_t blocksize)
+    template <class MatrixBlockType, class T>
+    __global__ void cuComputeDiagPtrsTyped(const int* rowIndices,
+                                           const int* colIndices,
+                                           int rows,
+                                           MatrixBlockType** diagPtrs,
+                                           T* values,
+                                           size_t blocksize)
     {
         const auto rawIdx = blockDim.x * blockIdx.x + threadIdx.x;
         if (rawIdx < rows) {
@@ -114,8 +115,7 @@ namespace
                 ++diagIdx;
             }
 
-            // Pointer arithmetic based on first element is sufficient here
-            diagPtrs[rawIdx] = values + diagIdx*blocksize*blocksize;
+            diagPtrs[rawIdx] = reinterpret_cast<MatrixBlockType*>(values + diagIdx * blocksize * blocksize);
         }
     }
 } // namespace
@@ -166,23 +166,6 @@ copyMatDataToReorderedSplit(const T* srcMatrix,
                                                                                  numberOfRows);
 }
 
-template <class T>
-GpuBuffer<T*>
-getDiagPtrs(GpuSparseMatrixWrapper<T>& matrix)
-{
-    GpuBuffer<T*> diagPtrs(matrix.N());
-
-    // Have cuda/hip automatically pick a reasonable block size for this kernel based on static analysis
-    int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuComputeDiagPtrs<T>);
-    // Calculate the number of blocks needed
-    int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(matrix.N(), threadBlockSize);
-
-    cuComputeDiagPtrs<T>
-        <<<nThreadBlocks, threadBlockSize>>>(
-            matrix->getRowIndices().data(), matrix->getColumnIndices().data(), matrix.N(), diagPtrs.data(), matrix->getNonZeroValues().data(), matrix.blockSize());
-    return diagPtrs;
-}
-
 #define INSTANTIATE_KERNEL_WRAPPERS_WITH_SCALAR_AND_BLOCKSIZE(T, blocksize)                                                                      \
     template void copyMatDataToReordered<T, blocksize>(const T*, const int*, T*, int*, int*, size_t, int);                         \
     template void copyMatDataToReorderedSplit<T, blocksize>(const T*, const int*, const int*, T*, int*, T*, int*, T*, int*, size_t, int);
@@ -202,11 +185,46 @@ INSTANTIATE_KERNEL_WRAPPERS_WITH_SCALAR_AND_BLOCKSIZE(double, 5);
 INSTANTIATE_KERNEL_WRAPPERS_WITH_SCALAR_AND_BLOCKSIZE(double, 6);
 INSTANTIATE_KERNEL_WRAPPERS_WITH_SCALAR_AND_BLOCKSIZE(double, 7);
 
-#define INSTANTIATE_KERNEL_WRAPPER_WITH_SCALAR(T) \
-    template Opm::gpuistl::GpuBuffer<T*> Opm::gpuistl::detail::getDiagPtrs<T>(                           \
+template <class MatrixBlockType, class T>
+GpuBuffer<MatrixBlockType*>
+getDiagPtrsTyped(GpuSparseMatrixWrapper<T>& matrix)
+{
+    GpuBuffer<MatrixBlockType*> diagPtrs(matrix.N());
+
+    int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuComputeDiagPtrsTyped<MatrixBlockType, T>);
+    int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(matrix.N(), threadBlockSize);
+
+    cuComputeDiagPtrsTyped<MatrixBlockType, T>
+        <<<nThreadBlocks, threadBlockSize>>>(
+            matrix->getRowIndices().data(), matrix->getColumnIndices().data(), matrix.N(),
+            diagPtrs.data(), matrix->getNonZeroValues().data(), matrix.blockSize());
+    return diagPtrs;
+}
+
+#define INSTANTIATE_GET_DIAG_PTRS_TYPED(T, blocksize) \
+    template Opm::gpuistl::GpuBuffer<Opm::gpuistl::MiniMatrix<T, blocksize>*>                \
+    Opm::gpuistl::detail::getDiagPtrsTyped<Opm::gpuistl::MiniMatrix<T, blocksize>, T>(        \
         Opm::gpuistl::GpuSparseMatrixWrapper<T>&);
 
-INSTANTIATE_KERNEL_WRAPPER_WITH_SCALAR(float)
-INSTANTIATE_KERNEL_WRAPPER_WITH_SCALAR(double)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(float, 1)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(float, 2)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(float, 3)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(float, 4)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(float, 5)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(float, 6)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(float, 7)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(double, 1)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(double, 2)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(double, 3)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(double, 4)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(double, 5)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(double, 6)
+INSTANTIATE_GET_DIAG_PTRS_TYPED(double, 7)
+
+#define INSTANTIATE_GET_DIAG_PTRS(T) \
+    template Opm::gpuistl::GpuBuffer<T*> Opm::gpuistl::detail::getDiagPtrsTyped<T, T>(GpuSparseMatrixWrapper<T>&);
+
+INSTANTIATE_GET_DIAG_PTRS(float)
+INSTANTIATE_GET_DIAG_PTRS(double)
 
 } // namespace Opm::gpuistl::detail
diff --git a/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp b/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp
index b3c018b70d7..6faa9532619 100644
--- a/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp
+++ b/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp
@@ -78,12 +78,18 @@ void copyMatDataToReorderedSplit(const T* srcMatrix,
                                  int threadBlockSize);
 
 /**
- * @brief Return the pointers to diagonal elements in a GpuBuffer
+ * @brief Return typed pointers to diagonal block elements in a GpuBuffer
+ *
+ * Like getDiagPtrs, but returns pointers typed as MatrixBlockType* instead of T*,
+ * avoiding the need for reinterpret_cast at the call sites.
+ *
+ * @tparam MatrixBlockType The block matrix type (e.g. MiniMatrix<T, N>) whose size
+ *         must match T[blockSize^2] in the underlying storage.
  * @param matrix The matrix to extract diagonal pointers from
- * @return GpuBuffer containing pointers to the diagonal elements
+ * @return GpuBuffer containing typed pointers to the diagonal block elements
  */
-template <class T>
-GpuBuffer<T*> getDiagPtrs(GpuSparseMatrixWrapper<T>& matrix);
+template <class MatrixBlockType, class T>
+GpuBuffer<MatrixBlockType*> getDiagPtrsTyped(GpuSparseMatrixWrapper<T>& matrix);
 
 } // namespace Opm::gpuistl::detail
 
diff --git a/tests/gpuistl/test_GpuSparseMatrix.cu b/tests/gpuistl/test_GpuSparseMatrix.cu
index f38df68909b..936b7556b57 100644
--- a/tests/gpuistl/test_GpuSparseMatrix.cu
+++ b/tests/gpuistl/test_GpuSparseMatrix.cu
@@ -177,7 +177,7 @@ BOOST_AUTO_TEST_CASE(TestGetDiagPtrs)
     BOOST_CHECK_EQUAL(cpuValues[4], 5.0);
     BOOST_CHECK_EQUAL(cpuValues[5], 6.0);
 
-    auto diagPtrs = Opm::gpuistl::detail::getDiagPtrs(gpuSparseMatrix);
+    auto diagPtrs = Opm::gpuistl::detail::getDiagPtrsTyped<double>(gpuSparseMatrix);
     auto diagPtrsView = Opm::gpuistl::GpuView<double*>(diagPtrs.data(), diagPtrs.size());
 
     auto d_result = Opm::gpuistl::make_gpu_unique_ptr<std::array<double, 3>>({-1});
@@ -199,12 +199,12 @@ __global__ void checkDiagBlocksKernel (double** diagPtrs, std::array<double, 12>
         (*values)[1] = diagPtrs[0][1];  // B[0][0][0][1]
         (*values)[2] = diagPtrs[0][2];  // B[0][0][1][0]
         (*values)[3] = diagPtrs[0][3];  // B[0][0][1][1]
-        
+
         (*values)[4] = diagPtrs[1][0];  // B[1][1][0][0]
         (*values)[5] = diagPtrs[1][1];  // B[1][1][0][1]
         (*values)[6] = diagPtrs[1][2];  // B[1][1][1][0]
         (*values)[7] = diagPtrs[1][3];  // B[1][1][1][1]
-        
+
         (*values)[8] = diagPtrs[2][0];  // B[2][2][0][0]
         (*values)[9] = diagPtrs[2][1];  // B[2][2][0][1]
         (*values)[10] = diagPtrs[2][2]; // B[2][2][1][0]
@@ -251,27 +251,27 @@ BOOST_AUTO_TEST_CASE(TestGetDiagPtrsBlockedNonScalar)
     B[0][0][0][1] = 1.1;
     B[0][0][1][0] = 1.2;
     B[0][0][1][1] = 1.3;
-    
+
     B[0][1][0][0] = 2.0;
     B[0][1][0][1] = 2.1;
     B[0][1][1][0] = 2.2;
     B[0][1][1][1] = 2.3;
-    
+
     B[0][2][0][0] = 3.0;
     B[0][2][0][1] = 3.1;
     B[0][2][1][0] = 3.2;
     B[0][2][1][1] = 3.3;
-    
+
     B[1][1][0][0] = 4.0;
     B[1][1][0][1] = 4.1;
     B[1][1][1][0] = 4.2;
     B[1][1][1][1] = 4.3;
-    
+
     B[2][0][0][0] = 5.0;
     B[2][0][0][1] = 5.1;
     B[2][0][1][0] = 5.2;
     B[2][0][1][1] = 5.3;
-    
+
     B[2][2][0][0] = 6.0;
     B[2][2][0][1] = 6.1;
     B[2][2][1][0] = 6.2;
@@ -279,8 +279,8 @@ BOOST_AUTO_TEST_CASE(TestGetDiagPtrsBlockedNonScalar)
 
     auto gpuSparseMatrix = Opm::gpuistl::GpuSparseMatrixWrapper<double>::fromMatrix(B);
 
-    auto diagPtrs = Opm::gpuistl::detail::getDiagPtrs(gpuSparseMatrix);
-    
+    auto diagPtrs = Opm::gpuistl::detail::getDiagPtrsTyped<double>(gpuSparseMatrix);
+
     // Check that we have the correct number of diagonal pointers
     BOOST_CHECK_EQUAL(diagPtrs.size(), N);
 
@@ -291,18 +291,18 @@ BOOST_AUTO_TEST_CASE(TestGetDiagPtrsBlockedNonScalar)
     checkDiagBlocksKernel<<<1, 1>>>(diagPtrsView.data(), d_result.get());
 
     std::array<double, 12> h_result = Opm::gpuistl::copyFromGPU(d_result);
-    
+
     // Verify diagonal block elements
     BOOST_CHECK_EQUAL(h_result[0], 1.0);
     BOOST_CHECK_EQUAL(h_result[1], 1.1);
     BOOST_CHECK_EQUAL(h_result[2], 1.2);
     BOOST_CHECK_EQUAL(h_result[3], 1.3);
-    
+
     BOOST_CHECK_EQUAL(h_result[4], 4.0);
     BOOST_CHECK_EQUAL(h_result[5], 4.1);
     BOOST_CHECK_EQUAL(h_result[6], 4.2);
     BOOST_CHECK_EQUAL(h_result[7], 4.3);
-    
+
     BOOST_CHECK_EQUAL(h_result[8], 6.0);
     BOOST_CHECK_EQUAL(h_result[9], 6.1);
     BOOST_CHECK_EQUAL(h_result[10], 6.2);
diff --git a/tests/gpuistl/test_GpuSparseTable.cu b/tests/gpuistl/test_GpuSparseTable.cu
index a58020d1ed5..b7b8c2ff8f8 100644
--- a/tests/gpuistl/test_GpuSparseTable.cu
+++ b/tests/gpuistl/test_GpuSparseTable.cu
@@ -42,7 +42,7 @@
 
 
 __global__ void fetchSparseTableValuesInKernel(Opm::SparseTable<int, Opm::gpuistl::GpuView> tableView,
-                                               std::array<int, 6>* output)
+                                               std::array<int, 7>* output)
 {
     int idx = blockIdx.x * blockDim.x + threadIdx.x;
     if (idx == 0) {
@@ -53,6 +53,7 @@ __global__ void fetchSparseTableValuesInKernel(Opm::SparseTable<int, Opm::gpuist
                 (*output)[offset++] = e;
             }
         }
+        (*output)[offset++] = tableView.dataSize();
     }
 }
 
@@ -83,11 +84,11 @@ BOOST_AUTO_TEST_CASE(TestUsingSparseTableInKernel)
     auto gpuBufferTable = Opm::gpuistl::copy_to_gpu<int>(cpuTable);
     auto gpuTableView = Opm::gpuistl::make_view(gpuBufferTable);
 
-    auto valuesFetchedFromGpu = Opm::gpuistl::make_gpu_unique_ptr<std::array<int, 6>>({-1});
+    auto valuesFetchedFromGpu = Opm::gpuistl::make_gpu_unique_ptr<std::array<int, 7>>({-1});
 
     fetchSparseTableValuesInKernel<<<1, 1>>>(gpuTableView, valuesFetchedFromGpu.get());
 
-    auto hostValues = Opm::gpuistl::copyFromGPU<std::array<int, 6>>(valuesFetchedFromGpu);
+    auto hostValues = Opm::gpuistl::copyFromGPU<std::array<int, 7>>(valuesFetchedFromGpu);
 
     BOOST_CHECK_EQUAL(hostValues[0], 1);
     BOOST_CHECK_EQUAL(hostValues[1], 2);
@@ -95,4 +96,5 @@ BOOST_AUTO_TEST_CASE(TestUsingSparseTableInKernel)
     BOOST_CHECK_EQUAL(hostValues[3], 4);
     BOOST_CHECK_EQUAL(hostValues[4], 5);
     BOOST_CHECK_EQUAL(hostValues[5], 6);
+    BOOST_CHECK_EQUAL(hostValues[6], 6);
 }
diff --git a/tests/gpuistl/test_blackoilintensivequantities_gpu.cu b/tests/gpuistl/test_blackoilintensivequantities_gpu.cu
new file mode 100644
index 00000000000..018ac4aaa04
--- /dev/null
+++ b/tests/gpuistl/test_blackoilintensivequantities_gpu.cu
@@ -0,0 +1,1098 @@
+/*
+  Copyright 2025 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/>.
+*/
+#include <config.h>
+
+// Tricks for older versions of ROCm that do not support alugrid
+// or dune fem. Note that these changes are only valid for this test file.
+#ifdef HAVE_DUNE_ALUGRID
+#undef HAVE_DUNE_ALUGRID
+#endif
+#define HAVE_DUNE_ALUGRID 0
+#ifdef HAVE_DUNE_FEM
+#undef HAVE_DUNE_FEM
+#endif
+#define HAVE_DUNE_FEM 0
+
+#include <opm/common/utility/gpuDecorators.hpp>
+#include <opm/material/common/ResetLocale.hpp>
+#include <opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp>
+#define HAVE_ECL_INPUT 1
+
+
+#define BOOST_TEST_MODULE TestBlackOilIntensiveQuantitiesGPU
+
+#include <boost/test/unit_test.hpp>
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/models/blackoil/blackoilmodel.hh>
+#include <opm/models/blackoil/blackoilprimaryvariables.hh>
+#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
+#include <opm/models/discretization/common/tpfalinearizer.hh>
+#include <opm/models/io/dgfvanguard.hh>
+#include <opm/models/utils/start.hh>
+#include <opm/simulators/flow/Main.hpp>
+#include <opm/simulators/linalg/gpuistl/detail/gpu_safe_call.hpp>
+#include <opm/simulators/linalg/gpuistl/gpu_smart_pointer.hpp>
+#include <opm/input/eclipse/Deck/Deck.hpp>
+#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/input/eclipse/Parser/Parser.hpp>
+#include <opm/input/eclipse/Python/Python.hpp>
+#include <opm/input/eclipse/Schedule/Schedule.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystemNonStatic.hpp>
+#include <opm/models/blackoil/blackoilintensivequantities.hh>
+#include <opm/models/nonlinear/newtonmethodparams.hpp>
+#include <opm/simulators/flow/equil/EquilibrationHelpers.hpp>
+#include <opm/simulators/flow/equil/InitStateEquil.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuBuffer.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuView.hpp>
+#include <opm/material/densead/Evaluation.hpp>
+#include <opm/models/utils/simulator.hh>
+#include <opm/simulators/flow/FlowProblemBlackoil.hpp>
+#include <opm/simulators/flow/FlowProblemBlackoilProperties.hpp>
+#include <opm/simulators/utils/moduleVersion.hpp>
+#include <opm/simulators/flow/BlackoilModelParameters.hpp>
+#include <opm/simulators/flow/FlowGenericVanguard.hpp>
+#include <opm/simulators/linalg/parallelbicgstabbackend.hh>
+#include <opm/simulators/wells/BlackoilWellModel.hpp>
+
+#include <opm/material/fluidmatrixinteractions/EclMaterialLawTwoPhaseTypes.hpp>
+#include <opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp>
+#include <opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp>
+#include <opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp>
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
+
+#include <opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp>
+#include <opm/material/thermal/GpuEclThermalLawManager.hpp>
+#include <opm/simulators/flow/GpuFlowProblem.hpp>
+
+#include <opm/simulators/flow/FlowGasWaterEnergyTypeTag.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuFlowGasWaterEnergyTypeTags.hpp>
+
+#include <chrono>
+#include <filesystem>
+#include <format>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+
+/// Build a CO2STORE WATER+GAS deck with the requested DIMENS. The grid is
+/// a uniform unit-cube grid; PERMX/PORO are constant. EQUIL gives the
+/// pressure profile so different cells (different depths) get different
+/// fluid states, exercising the per-cell update on both CPU and GPU.
+static std::string makeDeckString(int nx, int ny, int nz)
+{
+    const long total = static_cast<long>(nx) * ny * nz;
+    const long areal = static_cast<long>(nx) * ny;
+    std::ostringstream o;
+    o << "RUNSPEC\n"
+      << "DIMENS\n" << nx << ' ' << ny << ' ' << nz << " /\n"
+      << "EQLDIMS\n/\n"
+      << "TABDIMS\n/\n"
+      << "WATER\nGAS\nCO2STORE\nTHERMAL\nMETRIC\n"
+      << "GRID\n"
+      << "GRIDFILE\n0 0 /\n"
+      << "DX\n" << total << "*1 /\n"
+      << "DY\n" << total << "*1 /\n"
+      << "DZ\n" << total << "*1 /\n"
+      << "TOPS\n" << areal << "*0 /\n"
+      << "PERMX\n" << total << "*1013.25 /\n"
+      << "PORO\n" << total << "*0.25 /\n"
+      << "COPY\nPERMX PERMY /\nPERMX PERMZ /\n/\n"
+      << "THCONR\n" << total << "*1.5 /\n"
+      << "PROPS\n"
+      << "SGWFN\n"
+      << "0.000000E+00 0.000000E+00 1.000000E+00 3.060000E-02\n"
+      << "1.000000E+00 1.000000E+00 0.000000E+00 3.060000E-01 /\n"
+      << "SPECROCK\n20 2125.0  100 2125.0 /\n"
+      << "SOLUTION\n"
+      << "RPTRST\n'BASIC=0' /\n"
+      << "EQUIL\n0 300 100 0 0 0 1 1 0 /\n"
+      << "RTEMPVD\n0 36.12\n1000 70.0 /\n"
+      << "SCHEDULE\n"
+      << "RPTRST\n'BASIC=0' /\n"
+      << "TSTEP\n1 /";
+    return o.str();
+}
+
+
+static constexpr const char* deckString1 = "-- =============== RUNSPEC\n"
+                                           "RUNSPEC\n"
+                                           "DIMENS\n"
+                                           "3 3 3 /\n"
+                                           "EQLDIMS\n"
+                                           "/\n"
+                                           "TABDIMS\n"
+                                           "/\n"
+                                           "WATER\n"
+                                           "GAS\n"
+                                           "CO2STORE\n"
+                                           "THERMAL\n"
+                                           "METRIC\n"
+                                           "-- =============== GRID\n"
+                                           "GRID\n"
+                                           "GRIDFILE\n"
+                                           "0 0 /\n"
+                                           "DX\n"
+                                           "27*1 /\n"
+                                           "DY\n"
+                                           "27*1 /\n"
+                                           "DZ\n"
+                                           "27*1 /\n"
+                                           "TOPS\n"
+                                           "9*0 /\n"
+                                           "PERMX\n"
+                                           "27*1013.25 /\n"
+                                           "PORO\n"
+                                           "27*0.25 /\n"
+                                           "COPY\n"
+                                           "PERMX PERMY /\n"
+                                           "PERMX PERMZ /\n"
+                                           "/\n"
+                                           "THCONR\n"
+                                           "27*1.5 /\n"
+                                           "-- =============== PROPS\n"
+                                           "PROPS\n"
+                                           "SGWFN\n"
+                                           "0.000000E+00 0.000000E+00 1.000000E+00 3.060000E-02\n"
+                                           "1.000000E+00 1.000000E+00 0.000000E+00 3.060000E-01 /\n"
+                                           "SPECROCK\n"
+                                           "20 2125.0  100 2125.0 /\n"
+                                           "-- =============== SOLUTION\n"
+                                           "SOLUTION\n"
+                                           "RPTRST\n"
+                                           "'BASIC=0' /\n"
+                                           "EQUIL\n"
+                                           "0 300 100 0 0 0 1 1 0 /\n"
+                                           "RTEMPVD\n"
+                                           "0 36.12\n"
+                                           "1000 70.0 /\n"
+                                           "-- =============== SCHEDULE\n"
+                                           "SCHEDULE\n"
+                                           "RPTRST\n"
+                                           "'BASIC=0' /\n"
+                                           "TSTEP\n"
+                                           "1 /";
+
+using ScalarToUse = Opm::GetPropType<Opm::Properties::TTag::FlowProblem, Opm::Properties::Scalar>;
+
+using GpuB = Opm::gpuistl::GpuBuffer<double>;
+using GpuV = Opm::gpuistl::GpuView<double>;
+using FluidSystem = Opm::BlackOilFluidSystem<ScalarToUse>;
+using Evaluation = Opm::DenseAd::Evaluation<double, 2>;
+using BlackOilFluidSystemView
+    = Opm::BlackOilFluidSystemNonStatic<ScalarToUse,
+                                        Opm::BlackOilDefaultFluidSystemIndices,
+                                        Opm::gpuistl::GpuView>;
+
+template <class Value>
+OPM_HOST_DEVICE double asDouble(const Value& value)
+{
+    if constexpr (requires { value.value(); }) {
+        return static_cast<double>(value.value());
+    }
+    else {
+        return static_cast<double>(value);
+    }
+}
+
+/// Compare two scalar-or-Evaluation quantities. Always compares the value
+/// part. When \p checkDerivatives is true and both operands are
+/// `DenseAd::Evaluation`s, also checks every partial derivative using
+/// \c BOOST_CHECK_CLOSE with the supplied tolerance and increments
+/// \p derivativeComparisons by the number of partial derivatives compared.
+/// If \p checkDerivatives is true but the operands are not Evaluations,
+/// the test fails with \c BOOST_FAIL so a misconfigured derivative test
+/// cannot silently pass on plain scalars.
+template <class CpuValue, class GpuValue>
+void checkValueAndDerivatives(const CpuValue& cpuValue,
+                              const GpuValue& gpuValue,
+                              double tol,
+                              bool checkDerivatives,
+                              const char* label,
+                              std::size_t& derivativeComparisons)
+{
+    BOOST_CHECK_CLOSE(asDouble(cpuValue), asDouble(gpuValue), tol);
+
+    if (!checkDerivatives) {
+        return;
+    }
+
+    if constexpr (requires {
+                      CpuValue::numVars;
+                      cpuValue.derivative(0);
+                      GpuValue::numVars;
+                      gpuValue.derivative(0);
+                  })
+    {
+        static_assert(static_cast<int>(CpuValue::numVars)
+                          == static_cast<int>(GpuValue::numVars),
+                      "CPU and GPU Evaluation types must have the same numVars");
+        static_assert(static_cast<int>(CpuValue::numVars) > 0,
+                      "Derivative comparison requires at least one partial derivative");
+        for (int derivIdx = 0; derivIdx < CpuValue::numVars; ++derivIdx) {
+            const double cpuDerivative = static_cast<double>(cpuValue.derivative(derivIdx));
+            const double gpuDerivative = static_cast<double>(gpuValue.derivative(derivIdx));
+            BOOST_CHECK_MESSAGE(std::abs(cpuDerivative - gpuDerivative)
+                                    <= tol * (1.0 + std::abs(cpuDerivative)),
+                                label << " derivative[" << derivIdx
+                                      << "] cpu=" << cpuDerivative
+                                      << " gpu=" << gpuDerivative);
+            ++derivativeComparisons;
+        }
+    }
+    else {
+        BOOST_FAIL(std::string("checkDerivatives=true but ") + label
+                   + " is not a DenseAd::Evaluation");
+    }
+}
+
+template <class TypeTag>
+struct DummyProblem {
+
+
+    OPM_HOST_DEVICE DummyProblem() {
+        // empty
+    }
+
+    OPM_HOST_DEVICE ~DummyProblem() {
+        // empty
+    }
+
+    using EclMaterialLawManager =
+        typename Opm::GetProp<TypeTag, Opm::Properties::MaterialLaw>::EclMaterialLawManager;
+    using EclThermalLawManager =
+        typename Opm::GetProp<TypeTag, Opm::Properties::SolidEnergyLaw>::EclThermalLawManager;
+    struct MaterialLawParams {};
+
+    struct {
+        struct {
+            OPM_HOST_DEVICE Opm::LinearizationType getLinearizationType() const
+            {
+                return Opm::LinearizationType();
+            }
+        } lin_;
+
+        OPM_HOST_DEVICE auto linearizer() const
+        {
+            return lin_;
+        }
+    } model_;
+
+    OPM_HOST_DEVICE auto model() const
+    {
+        return model_;
+    }
+
+    OPM_HOST_DEVICE int satnumRegionIndex(std::size_t) const
+    {
+        return 0;
+    }
+    OPM_HOST_DEVICE const MaterialLawParams& materialLawParams(std::size_t) const
+    {
+        return materialLawParams_;
+    }
+    MaterialLawParams materialLawParams_;
+    OPM_HOST_DEVICE double rockCompressibility(std::size_t) const
+    {
+        return 0.0;
+    }
+    OPM_HOST_DEVICE double rockReferencePressure(std::size_t) const
+    {
+        return 0.0;
+    }
+    OPM_HOST_DEVICE double porosity(std::size_t, unsigned int) const
+    {
+        return 0.0;
+    }
+    OPM_HOST_DEVICE double maxOilVaporizationFactor(unsigned int, std::size_t) const
+    {
+        return 0.0;
+    }
+    OPM_HOST_DEVICE double maxGasDissolutionFactor(unsigned int, std::size_t) const
+    {
+        return 0.0;
+    }
+    OPM_HOST_DEVICE double maxOilSaturation(std::size_t) const
+    {
+        return 0.0;
+    }
+
+    template <class Evaluation>
+    OPM_HOST_DEVICE Evaluation rockCompPoroMultiplier(const auto&, std::size_t) const
+    {
+        return Evaluation(0.0);
+    }
+
+    template <class FluidState, class ...Args>
+    OPM_HOST_DEVICE void updateRelperms(auto& mobility, auto& dirMob, FluidState& fluidState, unsigned globalSpaceIdx) const
+    {
+    }
+
+    template <class Evaluation>
+    OPM_HOST_DEVICE Evaluation rockCompTransMultiplier(const auto&, std::size_t) const
+    {
+        return Evaluation(0.0);
+    }
+};
+
+namespace Opm
+{
+namespace Properties
+{
+namespace TTag
+{
+
+/// CPU side mirror of the simulation TypeTag, with diffusion/dispersion
+/// disabled because the simpler \c BlackOilIntensiveQuantities::update()
+/// overload used in this test does not support them. All other property
+/// bindings (Indices, MaterialLaw, FluidSystem, etc.) are inherited as-is
+/// from \c FlowGasWaterEnergyProblem.
+struct FlowGasWaterEnergyProblemTest {
+    using InheritsFrom = std::tuple<FlowGasWaterEnergyProblem>;
+};
+
+/// CPU TypeTag used by the real-deck test: identical to the production
+/// \c FlowGasWaterEnergyProblem (does not force-disable diffusion or
+/// dispersion) so that decks enabling those keywords can construct the
+/// simulator successfully.
+struct FlowGasWaterEnergyProblemTestRealDeck {
+    using InheritsFrom = std::tuple<FlowGasWaterEnergyProblem>;
+};
+
+} // namespace TTag
+
+template <class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::FlowGasWaterEnergyProblemTest>
+{ static constexpr bool value = false; };
+
+template <class TypeTag>
+struct EnableDispersion<TypeTag, TTag::FlowGasWaterEnergyProblemTest>
+{ static constexpr bool value = false; };
+
+} // namespace Properties
+
+} // namespace Opm
+
+using TypeTag = Opm::Properties::TTag::FlowGasWaterEnergyProblemTest;
+using TypeNacht = Opm::Properties::TTag::FlowGasWaterEnergyDummyProblemGPU;
+// FlowGasWaterEnergyKernelBaseGPU is the GPU-kernel-side TypeTag with GpuView
+// FluidSystem; previously named FlowGasWaterEnergyProblemGPU in the dispatcher
+// internals.  Use it wherever a GpuView-backed ScalarFluidState is needed.
+using TypeTagGPU = Opm::Properties::TTag::FlowGasWaterEnergyKernelBaseGPU;
+
+template<class IndexTraits, class GpuProblem>
+__global__ void
+testUsingOnGPU(Opm::BlackOilFluidSystemNonStatic<double, IndexTraits, Opm::gpuistl::GpuView> fluidSystem,
+               Opm::BlackOilIntensiveQuantities<TypeNacht> intensiveQuantities,
+               Opm::BlackOilPrimaryVariables<TypeNacht, Opm::gpuistl::MiniVector> primaryVariables,
+               GpuProblem problem)
+{
+    printf("fluidSystem.phaseIsActive(0): %d\n", fluidSystem.phaseIsActive(0));
+    printf("fluidSystem.phaseIsActive(1): %d\n", fluidSystem.phaseIsActive(1));
+    printf("fluidSystem.phaseIsActive(2): %d\n", fluidSystem.phaseIsActive(2));
+    using ScalarFluidState = typename Opm::BlackOilIntensiveQuantities<TypeTagGPU>::ScalarFluidState;
+    ScalarFluidState state(fluidSystem);
+    printf("BlackOilState density before update: %f\n", asDouble(state.density(0)));
+
+    ScalarFluidState state2(fluidSystem);
+    primaryVariables.assignNaive(state2);
+    printf("BlackOilState density before update: %f\n", asDouble(state.density(0)));
+    intensiveQuantities.updateSaturations(primaryVariables, 0, Opm::LinearizationType{});
+    intensiveQuantities.update(problem, primaryVariables, 0, 0);
+
+    printf("BlackOilState density after update: %f\n", asDouble(state.density(0)));
+}
+
+__global__ void dummykernel()
+{
+    printf("Hello from dummy kernel!\n");
+}
+
+template <class GpuProblem, class PrimaryVariables, class IntensiveQuantities>
+__global__ void
+updateAllCellsKernel(GpuProblem problem,
+                     Opm::gpuistl::GpuView<const PrimaryVariables> primaryVariables,
+                     Opm::gpuistl::GpuView<IntensiveQuantities> outIntensiveQuantities,
+                     std::size_t numCells)
+{
+    const std::size_t i = static_cast<std::size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
+    if (i >= numCells) {
+        return;
+    }
+    IntensiveQuantities intensiveQuantities = outIntensiveQuantities[i];
+    intensiveQuantities.updateSaturations(primaryVariables[i], 0, Opm::LinearizationType{});
+    intensiveQuantities.update(problem, primaryVariables[i], static_cast<unsigned>(i), 0);
+    intensiveQuantities.updateEnergyQuantities_(problem, static_cast<unsigned>(i), 0u);
+    outIntensiveQuantities[i] = intensiveQuantities;
+}
+
+/// One-shot global initialization of MPI, communicator, and parameter
+/// registration. Safe to call from multiple BOOST test cases; only the
+/// first call performs the work.
+static void initSimulatorOnce()
+{
+    static bool done = false;
+    if (done) {
+        return;
+    }
+    using namespace Opm;
+    int argc1 = boost::unit_test::framework::master_test_suite().argc;
+    char** argv1 = boost::unit_test::framework::master_test_suite().argv;
+#if HAVE_DUNE_FEM
+    Dune::Fem::MPIManager::initialize(argc1, argv1);
+#else
+    Dune::MPIHelper::instance(argc1, argv1);
+#endif
+    FlowGenericVanguard::setCommunication(std::make_unique<Opm::Parallel::Communication>());
+    Opm::ThreadManager::registerParameters();
+    Opm::NewtonMethodParams<double>::registerParameters();
+    BlackoilModelParameters<ScalarToUse>::registerParameters();
+    AdaptiveTimeStepping<TypeTag>::registerParameters();
+    Parameters::Register<Parameters::EnableTerminalOutput>(
+        "Dummy added for the well model to compile.");
+    registerAllParameters_<TypeTag>(/*finalizeRegistration=*/false);
+    // Also register parameters for the RealDeck TypeTag here so that all
+    // tests in the same process share a single (still-open) registration
+    // pass. Once endRegistration() (called via setupParameters_) closes
+    // registration, no further parameters can be added globally.
+    using RealDeckTypeTag = Opm::Properties::TTag::FlowGasWaterEnergyProblemTestRealDeck;
+    AdaptiveTimeStepping<RealDeckTypeTag>::registerParameters();
+    registerAllParameters_<RealDeckTypeTag>(/*finalizeRegistration=*/false);
+    Parameters::endRegistration();
+    done = true;
+}
+
+/// Run the per-cell intensive quantities update on both GPU and CPU for the
+/// deck stored at \p deckPath, optionally measuring wall-clock time on both
+/// sides.
+///
+/// Correctness is checked by comparing CPU vs GPU intensive quantities
+/// cell-by-cell at indices \p i = 0, \p sampleStride, 2*\p sampleStride, ...
+///
+/// \param deckPath        path to a CO2STORE WATER+GAS deck written to disk
+/// \param expectedNumCells number of grid cells the deck describes
+/// \param sampleStride    stride for sampled correctness; 1 means full check
+/// \param measureTiming   if true, prints CPU and GPU wall-clock to stdout
+static void runIntensiveQuantitiesTestForDeck(const std::string& deckPath,
+                                              std::size_t expectedNumCells,
+                                              std::size_t sampleStride,
+                                              bool measureTiming,
+                                              bool checkDerivatives = false)
+{
+    Opm::resetLocale();
+    initSimulatorOnce();
+
+    using namespace Opm;
+    using Simulator = Opm::GetPropType<TypeTag, Opm::Properties::Simulator>;
+
+    const auto filenameArg = std::string{"--ecl-deck-file-name="} + deckPath;
+    const char* argv2[] = {
+        "test_gpuflowproblem",
+        filenameArg.c_str(),
+        "--check-satfunc-consistency=false",
+    };
+    Opm::setupParameters_<TypeTag>(/*argc=*/sizeof(argv2) / sizeof(argv2[0]),
+                                   argv2,
+                                   /*registerParams=*/false,
+                                   /*allowUnused=*/true,
+                                   /*handleHelp=*/false,
+                                   /*myRank=*/0);
+
+    Opm::FlowGenericVanguard::readDeck(deckPath);
+    auto sim = std::make_unique<Simulator>();
+
+    auto& cpuProblem = sim->problem();
+    auto& dynamicFluidSystem = FluidSystem::getNonStaticInstance();
+
+    auto dynamicGpuFluidSystemBuffer  = ::Opm::gpuistl::copy_to_gpu(dynamicFluidSystem);
+    auto dynamicGpuFluidSystemView = ::Opm::gpuistl::make_view(dynamicGpuFluidSystemBuffer);
+
+    // Place the FluidSystemView in unified memory so that the device-side
+    // BlackOilFluidState pointer dereference is valid both on host and on
+    // device.
+    using FluidSystemViewType = std::decay_t<decltype(dynamicGpuFluidSystemView)>;
+    FluidSystemViewType* managedFluidSystemView = nullptr;
+    OPM_GPU_SAFE_CALL(cudaMallocManaged(&managedFluidSystemView, sizeof(FluidSystemViewType)));
+    new (managedFluidSystemView) FluidSystemViewType(dynamicGpuFluidSystemView);
+
+    using CpuMaterialLawManager = typename Opm::GetProp<TypeTag, Opm::Properties::MaterialLaw>::EclMaterialLawManager;
+    using Traits = typename CpuMaterialLawManager::MaterialLaw::Traits;
+    using TwoPhaseTraits = Opm::TwoPhaseMaterialTraits<ScalarToUse,
+                                                       Traits::wettingPhaseIdx,
+                                                       Traits::nonWettingPhaseIdx>;
+    using GpuPiecewiseLinearParams = Opm::PiecewiseLinearTwoPhaseMaterialParams<
+        TwoPhaseTraits, Opm::gpuistl::GpuView<const ScalarToUse>>;
+    using GpuPiecewiseLinearLaw = Opm::PiecewiseLinearTwoPhaseMaterial<TwoPhaseTraits, GpuPiecewiseLinearParams>;
+    using GpuMaterialLawParams = Opm::EclTwoPhaseMaterialParams<
+        Traits, GpuPiecewiseLinearParams, GpuPiecewiseLinearParams, GpuPiecewiseLinearParams,
+        Opm::gpuistl::ValueAsPointer>;
+    using GpuMaterialLaw = Opm::EclTwoPhaseMaterial<
+        Traits, GpuPiecewiseLinearLaw, GpuPiecewiseLinearLaw, GpuPiecewiseLinearLaw,
+        GpuMaterialLawParams>;
+    using GpuManagerBuffer = Opm::EclMaterialLaw::GpuManager<
+        Traits, GpuPiecewiseLinearLaw, GpuPiecewiseLinearLaw,
+        Opm::gpuistl::GpuBuffer, GpuMaterialLaw>;
+    using GpuThermalManagerBuffer = Opm::EclThermalLaw::GpuManager<
+        ScalarToUse, FluidSystemViewType,
+        Opm::gpuistl::GpuBuffer, Opm::gpuistl::GpuView>;
+    using GpuProblemBuffer = Opm::GpuFlowProblem<ScalarToUse, GpuManagerBuffer, Opm::gpuistl::GpuBuffer, GpuThermalManagerBuffer>;
+
+    GpuProblemBuffer gpuProblemBuffer(cpuProblem);
+    auto gpuProblemView = Opm::gpuistl::make_view(gpuProblemBuffer);
+
+    const std::size_t numCells = cpuProblem.model().numGridDof();
+    BOOST_REQUIRE_EQUAL(numCells, expectedNumCells);
+    BOOST_REQUIRE_GT(numCells, 0u);
+
+    using PrimaryVariablesCpu = Opm::GetPropType<TypeTag, Opm::Properties::PrimaryVariables>;
+    using PrimaryVariablesGpu = Opm::BlackOilPrimaryVariables<TypeNacht, Opm::gpuistl::MiniVector>;
+    using IntensiveQuantitiesCpu = Opm::BlackOilIntensiveQuantities<TypeTag>;
+    using IntensiveQuantitiesGpu = Opm::BlackOilIntensiveQuantities<TypeNacht>;
+
+    PrimaryVariablesCpu primaryVariablesCpu;
+    primaryVariablesCpu.setPrimaryVarsMeaningPressure(Opm::BlackOil::PressureMeaning::Pg);
+    PrimaryVariablesGpu primaryVariablesGpu(primaryVariablesCpu);
+    std::vector<PrimaryVariablesGpu> hostPrimaryVariablesGpu(numCells, primaryVariablesGpu);
+    Opm::gpuistl::GpuBuffer<PrimaryVariablesGpu> primaryVariablesBuffer(hostPrimaryVariablesGpu);
+
+    IntensiveQuantitiesCpu cpuIntensiveQuantitiesPrototype;
+    IntensiveQuantitiesGpu gpuIntensiveQuantitiesPrototype = cpuIntensiveQuantitiesPrototype.template withOtherFluidSystem<TypeNacht>(*managedFluidSystemView);
+    if (measureTiming) {
+        std::cout << std::format("[timing] sizeof(IntensiveQuantitiesGpu)={}  sizeof(IntensiveQuantitiesCpu)={}  "
+                                 "sizeof(PrimaryVariablesGpu)={}  cells={}  "
+                                 "approx GPU mem for IQ buffer={:.2f} MiB\n",
+                                 sizeof(IntensiveQuantitiesGpu), sizeof(IntensiveQuantitiesCpu),
+                                 sizeof(PrimaryVariablesGpu), numCells,
+                                 (sizeof(IntensiveQuantitiesGpu) * numCells) / (1024.0 * 1024.0));
+    }
+    std::vector<IntensiveQuantitiesGpu> hostIntensiveQuantities(numCells, gpuIntensiveQuantitiesPrototype);
+    Opm::gpuistl::GpuBuffer<IntensiveQuantitiesGpu> intensiveQuantitiesBuffer(hostIntensiveQuantities);
+
+    const unsigned blockSize = 64u;
+    const unsigned gridSize = static_cast<unsigned>((numCells + blockSize - 1u) / blockSize);
+
+    // Time the GPU kernel using CUDA events (excludes the host-to-device
+    // setup and device-to-host copy already done via GpuBuffer ctors above).
+    cudaEvent_t eventStart{}, eventStop{};
+    OPM_GPU_SAFE_CALL(cudaEventCreate(&eventStart));
+    OPM_GPU_SAFE_CALL(cudaEventCreate(&eventStop));
+    OPM_GPU_SAFE_CALL(cudaEventRecord(eventStart));
+    updateAllCellsKernel<<<gridSize, blockSize>>>(
+        gpuProblemView,
+        Opm::gpuistl::GpuView<const PrimaryVariablesGpu>(primaryVariablesBuffer.data(), primaryVariablesBuffer.size()),
+        Opm::gpuistl::GpuView<IntensiveQuantitiesGpu>(intensiveQuantitiesBuffer.data(), intensiveQuantitiesBuffer.size()),
+        numCells);
+    OPM_GPU_SAFE_CALL(cudaEventRecord(eventStop));
+    OPM_GPU_SAFE_CALL(cudaEventSynchronize(eventStop));
+    OPM_GPU_SAFE_CALL(cudaGetLastError());
+    float gpuMilliseconds = 0.0f;
+    OPM_GPU_SAFE_CALL(cudaEventElapsedTime(&gpuMilliseconds, eventStart, eventStop));
+    OPM_GPU_SAFE_CALL(cudaEventDestroy(eventStart));
+    OPM_GPU_SAFE_CALL(cudaEventDestroy(eventStop));
+
+    // CPU reference: same per-cell update on the host side using the full
+    // CPU problem; time it with std::chrono.
+    std::vector<IntensiveQuantitiesCpu> cpuIntensiveQuantities(numCells);
+    const auto cpuT0 = std::chrono::steady_clock::now();
+    for (std::size_t i = 0; i < numCells; ++i) {
+        cpuIntensiveQuantities[i].updateSaturations(primaryVariablesCpu, 0, Opm::LinearizationType{});
+        cpuIntensiveQuantities[i].update(cpuProblem, primaryVariablesCpu, static_cast<unsigned>(i), 0);
+    }
+    const auto cpuT1 = std::chrono::steady_clock::now();
+    const double cpuMilliseconds =
+        std::chrono::duration<double, std::milli>(cpuT1 - cpuT0).count();
+
+    if (measureTiming) {
+        std::cout << std::format("[timing] cells={}  CPU={:.3f} ms  GPU(kernel)={:.3f} ms  "
+                                 "speedup={:.2f}x\n",
+                                 numCells, cpuMilliseconds, gpuMilliseconds,
+                                 cpuMilliseconds / gpuMilliseconds);
+    }
+
+    // Bring GPU intensive quantities back and compare cells (sampled by sampleStride).
+    std::vector<IntensiveQuantitiesGpu> gpuIntensiveQuantitiesHost(numCells, gpuIntensiveQuantitiesPrototype);
+    OPM_GPU_SAFE_CALL(cudaMemcpy(gpuIntensiveQuantitiesHost.data(),
+                                 intensiveQuantitiesBuffer.data(),
+                                 numCells * sizeof(IntensiveQuantitiesGpu),
+                                 cudaMemcpyDeviceToHost));
+    BOOST_REQUIRE_EQUAL(gpuIntensiveQuantitiesHost.size(), numCells);
+
+    constexpr double tol = 1e-6;
+    // Active phase indices for CO2STORE WATER+GAS (oil disabled).
+    constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
+    constexpr unsigned gasPhaseIdx   = FluidSystem::gasPhaseIdx;
+    const unsigned activePhases[] = {waterPhaseIdx, gasPhaseIdx};
+
+    std::size_t checked = 0;
+    std::size_t derivativeComparisons = 0;
+    for (std::size_t i = 0; i < numCells; i += sampleStride) {
+        const auto& cpuIQ = cpuIntensiveQuantities[i];
+        const auto& gpuIQ = gpuIntensiveQuantitiesHost[i];
+        const auto& cpuFluidState = cpuIQ.fluidState();
+        const auto& gpuFluidState = gpuIQ.fluidState();
+
+        // Per-phase fluid state quantities for the active phases. Restricted
+        // to fields that are stored directly on the FluidState (no call into
+        // the underlying FluidSystem), since the GPU-side FluidState's
+        // FluidSystem holds device pointers (GpuView) that are not
+        // dereferenceable on the host.
+        for (unsigned p : activePhases) {
+            checkValueAndDerivatives(cpuFluidState.saturation(p),
+                                     gpuFluidState.saturation(p),
+                                     tol, checkDerivatives, "saturation",
+                                     derivativeComparisons);
+            checkValueAndDerivatives(cpuFluidState.pressure(p),
+                                     gpuFluidState.pressure(p),
+                                     tol, checkDerivatives, "pressure",
+                                     derivativeComparisons);
+            checkValueAndDerivatives(cpuFluidState.density(p),
+                                     gpuFluidState.density(p),
+                                     tol, checkDerivatives, "density",
+                                     derivativeComparisons);
+            checkValueAndDerivatives(cpuFluidState.invB(p),
+                                     gpuFluidState.invB(p),
+                                     tol, checkDerivatives, "invB",
+                                     derivativeComparisons);
+            // Mobility now goes through the opm-common EclTwoPhaseMaterial
+            // (relativePermeabilities) on both CPU and GPU, so the values
+            // must agree to the same tolerance as the rest of the IQ.
+            checkValueAndDerivatives(cpuIQ.mobility(p),
+                                     gpuIQ.mobility(p),
+                                     tol, checkDerivatives, "mobility",
+                                     derivativeComparisons);
+        }
+
+        // Scalar / per-cell quantities. Both CPU and GPU FluidState configs
+        // have enableVapwat=true and enableDisgasInWater=true so Rsw/Rvw are
+        // stored on both sides.
+        checkValueAndDerivatives(cpuFluidState.Rsw(), gpuFluidState.Rsw(),
+                                 tol, checkDerivatives, "Rsw",
+                                 derivativeComparisons);
+        checkValueAndDerivatives(cpuFluidState.Rvw(), gpuFluidState.Rvw(),
+                                 tol, checkDerivatives, "Rvw",
+                                 derivativeComparisons);
+        BOOST_CHECK_EQUAL(static_cast<unsigned>(cpuFluidState.pvtRegionIndex()),
+                          static_cast<unsigned>(gpuFluidState.pvtRegionIndex()));
+
+        checkValueAndDerivatives(cpuIQ.porosity(), gpuIQ.porosity(),
+                                 tol, checkDerivatives, "porosity",
+                                 derivativeComparisons);
+        BOOST_CHECK_CLOSE(static_cast<double>(cpuIQ.referencePorosity()),
+                          static_cast<double>(gpuIQ.referencePorosity()), tol);
+        checkValueAndDerivatives(cpuIQ.rockCompTransMultiplier(),
+                                 gpuIQ.rockCompTransMultiplier(),
+                                 tol, checkDerivatives, "rockCompTransMultiplier",
+                                 derivativeComparisons);
+
+        ++checked;
+    }
+    BOOST_TEST_MESSAGE("Per-cell GPU vs CPU IQ comparison: checked "
+                       << checked << " / " << numCells << " cells");
+    if (checkDerivatives) {
+        BOOST_TEST_MESSAGE("Derivative comparisons performed: " << derivativeComparisons);
+        // Guard against silent regressions: if checkDerivatives is requested
+        // but the helper never actually compared a derivative (e.g. because
+        // the IQ types degraded to plain scalars), fail loudly.
+        BOOST_REQUIRE_GT(derivativeComparisons, 0u);
+    }
+}
+
+/// RAII wrapper around a temporary file. The file is deleted when the
+/// object goes out of scope.
+struct TemporaryFile {
+    std::filesystem::path path;
+    explicit TemporaryFile(std::string_view filename)
+        : path(std::filesystem::temp_directory_path() / filename) {}
+    ~TemporaryFile() { std::filesystem::remove(path); }
+    TemporaryFile(const TemporaryFile&) = delete;
+    TemporaryFile& operator=(const TemporaryFile&) = delete;
+};
+
+/// Variant of \c runIntensiveQuantitiesTestForDeck that takes the per-cell
+/// primary variables from the simulator's initial (EQUIL-driven) solution
+/// vector instead of using a single uniform prototype. This exercises the
+/// exact code path used by the dispatcher in production:
+/// \c BlackOilPrimaryVariables<TypeNacht>(cpuPrimaryVariables[i]) -> GPU
+/// kernel update -> compare CPU vs GPU IQ field-by-field. \p sampleStride
+/// controls how often we check (1 means every cell).
+static void runIntensiveQuantitiesTestFromSimulatorSolution(const std::string& deckPath,
+                                                            std::size_t sampleStride)
+{
+    Opm::resetLocale();
+
+    using namespace Opm;
+    using LocalTypeTag    = Opm::Properties::TTag::FlowGasWaterEnergyProblemTestRealDeck;
+    using LocalGpuTypeTag = Opm::Properties::TTag::FlowGasWaterEnergyDummyProblemGPU;
+    using Simulator = Opm::GetPropType<LocalTypeTag, Opm::Properties::Simulator>;
+
+    // Reuse the global init: it registers parameters for both TypeTags so
+    // this test works whether it runs first or after another test in the
+    // same process (parameter registration is closed after the first
+    // setupParameters_ call).
+    initSimulatorOnce();
+
+    const auto filenameArg = std::string{"--ecl-deck-file-name="} + deckPath;
+    const char* argv2[] = {
+        "test_gpuflowproblem",
+        filenameArg.c_str(),
+        "--check-satfunc-consistency=false",
+        "--enable-tuning=false",
+    };
+    Opm::setupParameters_<LocalTypeTag>(/*argc=*/sizeof(argv2) / sizeof(argv2[0]),
+                                        argv2,
+                                        /*registerParams=*/false,
+                                        /*allowUnused=*/true,
+                                        /*handleHelp=*/false,
+                                        /*myRank=*/0);
+
+    Opm::FlowGenericVanguard::readDeck(deckPath);
+    auto sim = std::make_unique<Simulator>();
+
+    // The Simulator constructor only calls finishInit; the solution vector
+    // is not yet EQUIL-populated. Trigger initial-solution application so
+    // model().solution(0)[i] holds the real per-cell primary variables.
+    sim->model().applyInitialSolution();
+
+    auto& cpuProblem = sim->problem();
+    auto& dynamicFluidSystem = FluidSystem::getNonStaticInstance();
+
+    BOOST_TEST_MESSAGE(std::format(
+        "FluidSystem phase activation: static[O={}, W={}, G={}]  dynamic[O={}, W={}, G={}]",
+        FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),
+        FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
+        FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx),
+        dynamicFluidSystem.phaseIsActive(FluidSystem::oilPhaseIdx),
+        dynamicFluidSystem.phaseIsActive(FluidSystem::waterPhaseIdx),
+        dynamicFluidSystem.phaseIsActive(FluidSystem::gasPhaseIdx)));
+
+    auto dynamicGpuFluidSystemBuffer  = ::Opm::gpuistl::copy_to_gpu(dynamicFluidSystem);
+    auto dynamicGpuFluidSystemView = ::Opm::gpuistl::make_view(dynamicGpuFluidSystemBuffer);
+
+    using FluidSystemViewType = std::decay_t<decltype(dynamicGpuFluidSystemView)>;
+    FluidSystemViewType* managedFluidSystemView = nullptr;
+    OPM_GPU_SAFE_CALL(cudaMallocManaged(&managedFluidSystemView, sizeof(FluidSystemViewType)));
+    new (managedFluidSystemView) FluidSystemViewType(dynamicGpuFluidSystemView);
+
+    BOOST_TEST_MESSAGE(std::format(
+        "GPU fluid-system view phase activation: O={} W={} G={}",
+        managedFluidSystemView->phaseIsActive(FluidSystem::oilPhaseIdx),
+        managedFluidSystemView->phaseIsActive(FluidSystem::waterPhaseIdx),
+        managedFluidSystemView->phaseIsActive(FluidSystem::gasPhaseIdx)));
+
+    using CpuMaterialLawManager = typename Opm::GetProp<LocalTypeTag, Opm::Properties::MaterialLaw>::EclMaterialLawManager;
+    using Traits = typename CpuMaterialLawManager::MaterialLaw::Traits;
+    using TwoPhaseTraits = Opm::TwoPhaseMaterialTraits<ScalarToUse,
+                                                       Traits::wettingPhaseIdx,
+                                                       Traits::nonWettingPhaseIdx>;
+    using GpuPiecewiseLinearParams = Opm::PiecewiseLinearTwoPhaseMaterialParams<
+        TwoPhaseTraits, Opm::gpuistl::GpuView<const ScalarToUse>>;
+    using GpuPiecewiseLinearLaw = Opm::PiecewiseLinearTwoPhaseMaterial<TwoPhaseTraits, GpuPiecewiseLinearParams>;
+    using GpuMaterialLawParams = Opm::EclTwoPhaseMaterialParams<
+        Traits, GpuPiecewiseLinearParams, GpuPiecewiseLinearParams, GpuPiecewiseLinearParams,
+        Opm::gpuistl::ValueAsPointer>;
+    using GpuMaterialLaw = Opm::EclTwoPhaseMaterial<
+        Traits, GpuPiecewiseLinearLaw, GpuPiecewiseLinearLaw, GpuPiecewiseLinearLaw,
+        GpuMaterialLawParams>;
+    using GpuManagerBuffer = Opm::EclMaterialLaw::GpuManager<
+        Traits, GpuPiecewiseLinearLaw, GpuPiecewiseLinearLaw,
+        Opm::gpuistl::GpuBuffer, GpuMaterialLaw>;
+    using GpuThermalManagerBuffer = Opm::EclThermalLaw::GpuManager<
+        ScalarToUse, FluidSystemViewType,
+        Opm::gpuistl::GpuBuffer, Opm::gpuistl::GpuView>;
+    using GpuProblemBuffer = Opm::GpuFlowProblem<ScalarToUse, GpuManagerBuffer, Opm::gpuistl::GpuBuffer, GpuThermalManagerBuffer>;
+
+    GpuProblemBuffer gpuProblemBuffer(cpuProblem);
+    auto gpuProblemView = Opm::gpuistl::make_view(gpuProblemBuffer);
+
+    const std::size_t numCells = cpuProblem.model().numGridDof();
+    BOOST_REQUIRE_GT(numCells, 0u);
+    BOOST_TEST_MESSAGE("Real-deck test: numCells=" << numCells);
+
+    using PrimaryVariablesCpu = Opm::GetPropType<LocalTypeTag, Opm::Properties::PrimaryVariables>;
+    using PrimaryVariablesGpu = Opm::BlackOilPrimaryVariables<LocalGpuTypeTag, Opm::gpuistl::MiniVector>;
+    using IntensiveQuantitiesCpu = Opm::BlackOilIntensiveQuantities<LocalTypeTag>;
+    using IntensiveQuantitiesGpu = Opm::BlackOilIntensiveQuantities<LocalGpuTypeTag>;
+
+    // Pull the EQUIL-initialized per-cell primary variables straight from
+    // the simulator's solution vector and convert them to the GPU
+    // PrimaryVariables type cell-by-cell.
+    const auto& cpuSolution = cpuProblem.model().solution(/*timeIdx=*/0);
+    BOOST_REQUIRE_EQUAL(cpuSolution.size(), numCells);
+    std::vector<PrimaryVariablesCpu> cpuPrimaryVariables(numCells);
+    std::vector<PrimaryVariablesGpu> hostPrimaryVariablesGpu;
+    hostPrimaryVariablesGpu.reserve(numCells);
+    for (std::size_t i = 0; i < numCells; ++i) {
+        cpuPrimaryVariables[i] = cpuSolution[i];
+        hostPrimaryVariablesGpu.emplace_back(cpuPrimaryVariables[i]);
+    }
+    Opm::gpuistl::GpuBuffer<PrimaryVariablesGpu> primaryVariablesBuffer(hostPrimaryVariablesGpu);
+
+    // Diagnostic: distribution of pressure-meaning values across cells.
+    {
+        std::array<std::size_t, 3> meaningCounts{0, 0, 0};
+        for (std::size_t i = 0; i < numCells; ++i) {
+            const auto m = cpuPrimaryVariables[i].primaryVarsMeaningPressure();
+            switch (m) {
+                case Opm::BlackOil::PressureMeaning::Pg: ++meaningCounts[0]; break;
+                case Opm::BlackOil::PressureMeaning::Po: ++meaningCounts[1]; break;
+                case Opm::BlackOil::PressureMeaning::Pw: ++meaningCounts[2]; break;
+                default: break;
+            }
+        }
+        BOOST_TEST_MESSAGE(std::format(
+            "PressureMeaning distribution (cells={}): Pg={} Po={} Pw={}",
+            numCells, meaningCounts[0], meaningCounts[1], meaningCounts[2]));
+    }
+    if (numCells > 0) {
+        const auto& pv0 = cpuPrimaryVariables[0];
+        BOOST_TEST_MESSAGE(std::format(
+            "Cell 0 PV: meaning(P={},W={},G={}) pvtRegion={}",
+            static_cast<int>(pv0.primaryVarsMeaningPressure()),
+            static_cast<int>(pv0.primaryVarsMeaningWater()),
+            static_cast<int>(pv0.primaryVarsMeaningGas()),
+            pv0.pvtRegionIndex()));
+        for (unsigned k = 0; k < pv0.size(); ++k) {
+            BOOST_TEST_MESSAGE(std::format("  pv[{}] = {}", k, double(pv0[k])));
+        }
+    }
+
+
+    IntensiveQuantitiesCpu cpuIntensiveQuantitiesPrototype;
+    IntensiveQuantitiesGpu gpuIntensiveQuantitiesPrototype =
+        cpuIntensiveQuantitiesPrototype.template withOtherFluidSystem<TypeNacht>(*managedFluidSystemView);
+    std::vector<IntensiveQuantitiesGpu> hostIntensiveQuantities(numCells, gpuIntensiveQuantitiesPrototype);
+    Opm::gpuistl::GpuBuffer<IntensiveQuantitiesGpu> intensiveQuantitiesBuffer(hostIntensiveQuantities);
+
+    const unsigned blockSize = 64u;
+    const unsigned gridSize = static_cast<unsigned>((numCells + blockSize - 1u) / blockSize);
+
+    updateAllCellsKernel<<<gridSize, blockSize>>>(
+        gpuProblemView,
+        Opm::gpuistl::GpuView<const PrimaryVariablesGpu>(
+            primaryVariablesBuffer.data(), primaryVariablesBuffer.size()),
+        Opm::gpuistl::GpuView<IntensiveQuantitiesGpu>(
+            intensiveQuantitiesBuffer.data(), intensiveQuantitiesBuffer.size()),
+        numCells);
+    OPM_GPU_SAFE_CALL(cudaDeviceSynchronize());
+    OPM_GPU_SAFE_CALL(cudaGetLastError());
+
+    // CPU reference: use the IntensiveQuantities the simulator itself
+    // computed during invalidateAndUpdateIntensiveQuantities(0). This is the
+    // exact CPU value that the production dispatcher path would overlay
+    // GPU-computed BlackOil fields onto, so it is a meaningful reference for
+    // the GPU comparison without requiring the (diffusion-incompatible)
+    // simpler IQ::update overload to compile here.
+    sim->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
+    std::vector<IntensiveQuantitiesCpu> cpuIntensiveQuantities(numCells);
+    for (std::size_t i = 0; i < numCells; ++i) {
+        const auto* cached = sim->model().cachedIntensiveQuantities(static_cast<unsigned>(i), 0);
+        BOOST_REQUIRE(cached != nullptr);
+        cpuIntensiveQuantities[i] = *cached;
+    }
+
+    std::vector<IntensiveQuantitiesGpu> gpuIntensiveQuantitiesHost(numCells, gpuIntensiveQuantitiesPrototype);
+    OPM_GPU_SAFE_CALL(cudaMemcpy(gpuIntensiveQuantitiesHost.data(),
+                                 intensiveQuantitiesBuffer.data(),
+                                 numCells * sizeof(IntensiveQuantitiesGpu),
+                                 cudaMemcpyDeviceToHost));
+
+    constexpr double tol = 1e-6;
+    constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
+    constexpr unsigned gasPhaseIdx   = FluidSystem::gasPhaseIdx;
+    const unsigned activePhases[] = {waterPhaseIdx, gasPhaseIdx};
+
+    std::size_t checked = 0;
+    std::size_t derivativeComparisons = 0;
+    std::size_t firstBadCell = static_cast<std::size_t>(-1);
+    for (std::size_t i = 0; i < numCells; i += sampleStride) {
+        const auto& cpuIQ = cpuIntensiveQuantities[i];
+        const auto& gpuIQ = gpuIntensiveQuantitiesHost[i];
+        const auto& cpuFs = cpuIQ.fluidState();
+        const auto& gpuFs = gpuIQ.fluidState();
+
+        // Probe non-finite values explicitly; report the first offender so
+        // we can debug the production simulator failure mode.
+        for (unsigned p : activePhases) {
+            const double cpuRho = asDouble(cpuFs.density(p));
+            const double gpuRho = asDouble(gpuFs.density(p));
+            const double cpuP   = asDouble(cpuFs.pressure(p));
+            const double gpuP   = asDouble(gpuFs.pressure(p));
+            const double cpuS   = asDouble(cpuFs.saturation(p));
+            const double gpuS   = asDouble(gpuFs.saturation(p));
+            if (!std::isfinite(gpuRho) || !std::isfinite(gpuP) || !std::isfinite(gpuS)) {
+                if (firstBadCell == static_cast<std::size_t>(-1)) {
+                    firstBadCell = i;
+                    BOOST_TEST_MESSAGE(std::format(
+                        "First non-finite GPU IQ at cell {} (phase {}): "
+                        "S(cpu={}, gpu={}) P(cpu={}, gpu={}) rho(cpu={}, gpu={})",
+                        i, p, cpuS, gpuS, cpuP, gpuP, cpuRho, gpuRho));
+                }
+            }
+        }
+
+        for (unsigned p : activePhases) {
+            checkValueAndDerivatives(cpuFs.saturation(p), gpuFs.saturation(p),
+                                     tol, /*checkDerivatives=*/false, "saturation",
+                                     derivativeComparisons);
+            checkValueAndDerivatives(cpuFs.pressure(p), gpuFs.pressure(p),
+                                     tol, /*checkDerivatives=*/false, "pressure",
+                                     derivativeComparisons);
+            checkValueAndDerivatives(cpuFs.density(p), gpuFs.density(p),
+                                     tol, /*checkDerivatives=*/false, "density",
+                                     derivativeComparisons);
+            checkValueAndDerivatives(cpuFs.invB(p), gpuFs.invB(p),
+                                     tol, /*checkDerivatives=*/false, "invB",
+                                     derivativeComparisons);
+            // mobility is intentionally NOT compared: the GPU dispatcher
+            // does not overlay it (left to the CPU-computed value in
+            // production).
+        }
+        checkValueAndDerivatives(cpuFs.Rsw(), gpuFs.Rsw(),
+                                 tol, /*checkDerivatives=*/false, "Rsw",
+                                 derivativeComparisons);
+        checkValueAndDerivatives(cpuFs.Rvw(), gpuFs.Rvw(),
+                                 tol, /*checkDerivatives=*/false, "Rvw",
+                                 derivativeComparisons);
+        BOOST_CHECK_EQUAL(static_cast<unsigned>(cpuFs.pvtRegionIndex()),
+                          static_cast<unsigned>(gpuFs.pvtRegionIndex()));
+        checkValueAndDerivatives(cpuIQ.porosity(), gpuIQ.porosity(),
+                                 tol, /*checkDerivatives=*/false, "porosity",
+                                 derivativeComparisons);
+        BOOST_CHECK_CLOSE(static_cast<double>(cpuIQ.referencePorosity()),
+                          static_cast<double>(gpuIQ.referencePorosity()), tol);
+
+        // Thermal IQ fields: temperature in the fluid state, and the three
+        // BlackOilEnergyIntensiveQuantities scalars (rockInternalEnergy,
+        // totalThermalConductivity, rockFraction). For FullyImplicitThermal
+        // the GPU dispatcher writes all of these via overlayBlackOilFieldsFrom,
+        // so they should match the CPU value field-by-field.
+        checkValueAndDerivatives(cpuFs.temperature(/*phaseIdx=*/0),
+                                 gpuFs.temperature(/*phaseIdx=*/0),
+                                 tol, /*checkDerivatives=*/false, "temperature",
+                                 derivativeComparisons);
+        checkValueAndDerivatives(cpuIQ.rockInternalEnergy(),
+                                 gpuIQ.rockInternalEnergy(),
+                                 tol, /*checkDerivatives=*/false, "rockInternalEnergy",
+                                 derivativeComparisons);
+        checkValueAndDerivatives(cpuIQ.totalThermalConductivity(),
+                                 gpuIQ.totalThermalConductivity(),
+                                 tol, /*checkDerivatives=*/false, "totalThermalConductivity",
+                                 derivativeComparisons);
+        BOOST_CHECK_CLOSE(static_cast<double>(cpuIQ.rockFraction()),
+                          static_cast<double>(gpuIQ.rockFraction()), tol);
+        for (unsigned p : activePhases) {
+            checkValueAndDerivatives(cpuFs.enthalpy(p), gpuFs.enthalpy(p),
+                                     tol, /*checkDerivatives=*/false, "enthalpy",
+                                     derivativeComparisons);
+        }
+        ++checked;
+    }
+    BOOST_TEST_MESSAGE("Real-deck per-cell GPU vs CPU IQ comparison: checked "
+                       << checked << " / " << numCells << " cells");
+}
+
+BOOST_AUTO_TEST_CASE(TestRealDeckGpuVsCpuFromSimulatorSolution)
+{
+    // Drive the per-cell IQ update with the EQUIL-initialized primary
+    // variables of the production CO2STORE deck the user reports diverges
+    // when the dispatcher is enabled. This reproduces the exact code path
+    // taken by the dispatcher in production (per-cell PrimaryVariables
+    // upload + GPU update + readback) and compares every IQ field against
+    // a CPU reference, so we can identify which field becomes non-finite.
+    //
+    // NOTE: this test must run BEFORE the synthetic-deck tests because the
+    // BlackOilFluidSystem singleton is global state that gets re-set by
+    // every readDeck call; running this after the synthetic tests in the
+    // same process leaves the static FluidSystem in a state that is
+    // inconsistent with the per-cell GPU FluidSystem copy taken here.
+    const std::string deckPath = "/workspaces/opm/thecaseiwant/deck/THECASEIWANT.DATA";
+    if (!std::filesystem::exists(deckPath)) {
+        BOOST_TEST_MESSAGE("Skipping: deck not found at " << deckPath);
+        return;
+    }
+    runIntensiveQuantitiesTestFromSimulatorSolution(deckPath, /*sampleStride=*/1u);
+}
+
+BOOST_AUTO_TEST_CASE(TestInstantiateGpuFlowProblem)
+{
+    // Persist the in-source two-phase WATER+GAS+CO2STORE deck to a file so
+    // the standard Vanguard read path can ingest it.
+    const TemporaryFile tempFile("test_blackoilintensivequantities_gpu.DATA");
+    {
+        std::ofstream f(tempFile.path);
+        f << deckString1;
+    }
+    runIntensiveQuantitiesTestForDeck(tempFile.path.string(),
+                                      /*expectedNumCells=*/27u,
+                                      /*sampleStride=*/1u,
+                                      /*measureTiming=*/false);
+}
+
+BOOST_AUTO_TEST_CASE(TestPerformanceGpuVsCpuLargeGrid)
+{
+    // 100^3 = 1'000'000 cells. Build a CO2STORE WATER+GAS deck with a uniform
+    // unit-cube grid and EQUIL initialization, then time the per-cell IQ
+    // update on both GPU and CPU. Correctness is checked at every 1000th
+    // cell (1000 sample points) so the test stays reasonably fast.
+    // 64^3 = 262'144 cells. Build a CO2STORE WATER+GAS deck with a uniform
+    // unit-cube grid and EQUIL initialization, then time the per-cell IQ
+    // update on both GPU and CPU. Correctness is checked at every 256th cell
+    // (~1024 sample points) so the test stays reasonably fast.
+    //
+    // NOTE: the upper grid size here is bounded by the per-cell device
+    // allocation pattern in GpuEclMaterialLawManager (one cudaMalloc per
+    // piecewise-linear sample array per cell, i.e. ~12 cudaMalloc per cell).
+    // With ROCm's typical 4 KiB allocation granularity, ~16 GiB of VRAM
+    // limits us to a few hundred thousand cells; a 1M-cell deck would
+    // require consolidating per-cell sample buffers (left as a future
+    // optimisation, since for a homogeneous deck all cells share the same
+    // SGWFN table).
+    constexpr int dim = 64;
+    constexpr std::size_t expected =
+        static_cast<std::size_t>(dim) * dim * dim;
+    const TemporaryFile tempFile("test_blackoilintensivequantities_gpu_1M.DATA");
+    {
+        std::ofstream f(tempFile.path);
+        f << makeDeckString(dim, dim, dim);
+    }
+    runIntensiveQuantitiesTestForDeck(tempFile.path.string(),
+                                      /*expectedNumCells=*/expected,
+                                      /*sampleStride=*/256u,
+                                      /*measureTiming=*/true);
+}
+
+BOOST_AUTO_TEST_CASE(TestGpuVsCpuIntensiveQuantitiesEvaluationDerivatives)
+{
+    // Re-runs the small (3x3x3) WATER+GAS+CO2STORE deck used by
+    // TestInstantiateGpuFlowProblem, but in addition to the value of every
+    // intensive-quantity field this test also compares every partial
+    // derivative of the underlying DenseAd::Evaluation between CPU and GPU.
+    // This guarantees that the GPU material-law / fluid-system evaluation
+    // chain is differentiated identically to the CPU one.
+    const TemporaryFile tempFile("test_blackoilintensivequantities_gpu_deriv.DATA");
+    {
+        std::ofstream f(tempFile.path);
+        f << deckString1;
+    }
+    runIntensiveQuantitiesTestForDeck(tempFile.path.string(),
+                                      /*expectedNumCells=*/27u,
+                                      /*sampleStride=*/1u,
+                                      /*measureTiming=*/false,
+                                      /*checkDerivatives=*/true);
+}
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..290a64a9edd
--- /dev/null
+++ b/tests/gpuistl/test_gpu_ecl_thermal_law_manager.cu
@@ -0,0 +1,214 @@
+/*
+  Copyright TODO ADD YEAR AND NAME OF AUTHOR
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#include <config.h>
+
+#define BOOST_TEST_MODULE TestGpuEclThermalLawManager
+
+#include <boost/test/unit_test.hpp>
+
+#include <cuda_runtime.h>
+
+#include <opm/material/common/EnsureFinalized.hpp>
+#include <opm/material/common/Tabulated1DFunction.hpp>
+#include <opm/material/densead/Evaluation.hpp>
+#include <opm/material/fluidsystems/BlackOilDefaultFluidSystemIndices.hpp>
+#include <opm/material/thermal/EclSpecrockLaw.hpp>
+#include <opm/material/thermal/EclSpecrockLawParams.hpp>
+#include <opm/material/thermal/EclThconrLaw.hpp>
+#include <opm/material/thermal/EclThconrLawParams.hpp>
+
+#include <opm/material/thermal/GpuEclThermalLawManager.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuBuffer.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuView.hpp>
+#include <opm/simulators/linalg/gpuistl/detail/gpu_safe_call.hpp>
+
+#include <array>
+#include <vector>
+
+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 <class ScalarT>
+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<Scalar, Opm::gpuistl::GpuView>;
+using GpuSpecrockParamsView = Opm::EclSpecrockLawParams<Scalar, Opm::gpuistl::GpuView>;
+using GpuSpecrockLawView = Opm::EclSpecrockLaw<Scalar, GpuSpecrockParamsView>;
+using GpuThconrLaw = Opm::EclThconrLaw<Scalar, DummyFluidSystem>;
+
+using ManagerCpu = Opm::EclThermalLaw::GpuManager<Scalar, DummyFluidSystem>;
+using ManagerBuf
+    = Opm::EclThermalLaw::GpuManager<Scalar, DummyFluidSystem, Opm::gpuistl::GpuBuffer, Opm::gpuistl::GpuView>;
+using ManagerView
+    = Opm::EclThermalLaw::GpuManager<Scalar, DummyFluidSystem, Opm::gpuistl::GpuView, Opm::gpuistl::GpuView>;
+
+// 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<std::size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
+    if (i >= numCells) {
+        return;
+    }
+    DummyFluidState<Scalar> fs;
+    fs.temperature_ = temperatures[i];
+    fs.gasSaturation_ = gasSaturations[i];
+
+    const auto solidEnergyParams = view.solidEnergyLawParams(static_cast<unsigned>(i));
+    outSolidEnergy[i] =
+        GpuSpecrockLawView::solidInternalEnergy(solidEnergyParams, fs);
+
+    const auto thermalConductionParams =
+        view.thermalConductionLawParams(static_cast<unsigned>(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<Scalar> regionTemperatures0 = { 280.0, 320.0, 360.0, 400.0 };
+    const std::vector<Scalar> regionEnergies0     = { 0.0,   4.0e7, 9.0e7, 1.5e8 };
+    const std::vector<Scalar> regionTemperatures1 = { 280.0, 350.0, 420.0 };
+    const std::vector<Scalar> regionEnergies1     = { 0.0,   3.5e7, 8.4e7 };
+
+    // Build a host-side CPU manager directly from the basic constructor.
+    Opm::EclSpecrockLawParams<Scalar> region0Params;
+    region0Params.setSamples(regionTemperatures0, regionEnergies0);
+    Opm::EclSpecrockLawParams<Scalar> region1Params;
+    region1Params.setSamples(regionTemperatures1, regionEnergies1);
+
+    std::vector<Opm::EclSpecrockLawParams<Scalar>> cpuSolidEnergyParams;
+    cpuSolidEnergyParams.push_back(region0Params);
+    cpuSolidEnergyParams.push_back(region1Params);
+
+    std::vector<int> cpuElementToRegion(numCells);
+    for (std::size_t i = 0; i < numCells; ++i) {
+        cpuElementToRegion[i] = (i < 3) ? 0 : 1;
+    }
+
+    std::vector<Opm::EclThconrLawParams<Scalar>> cpuThermalConductionParams(numCells);
+    for (std::size_t i = 0; i < numCells; ++i) {
+        cpuThermalConductionParams[i].setReferenceTotalThermalConductivity(
+            100.0 + 5.0 * static_cast<Scalar>(i));
+        cpuThermalConductionParams[i].setDTotalThermalConductivity_dSg(
+            0.20 + 0.01 * static_cast<Scalar>(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<Scalar, numCells> hostTemperatures = {
+        290.0, 310.0, 355.0, 285.0, 360.0, 410.0
+    };
+    std::array<Scalar, numCells> hostGasSaturations = {
+        0.0, 0.1, 0.25, 0.5, 0.75, 0.99
+    };
+
+    // -- Run the kernel.
+    Scalar* dT = nullptr;
+    Scalar* dSg = nullptr;
+    Scalar* dEnergy = nullptr;
+    Scalar* dConductivity = nullptr;
+    OPM_GPU_SAFE_CALL(cudaMalloc(&dT, numCells * sizeof(Scalar)));
+    OPM_GPU_SAFE_CALL(cudaMalloc(&dSg, numCells * sizeof(Scalar)));
+    OPM_GPU_SAFE_CALL(cudaMalloc(&dEnergy, numCells * sizeof(Scalar)));
+    OPM_GPU_SAFE_CALL(cudaMalloc(&dConductivity, numCells * sizeof(Scalar)));
+    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<Scalar, numCells> hostEnergyResult {};
+    std::array<Scalar, numCells> 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));
+
+    OPM_GPU_SAFE_CALL(cudaFree(dT));
+    OPM_GPU_SAFE_CALL(cudaFree(dSg));
+    OPM_GPU_SAFE_CALL(cudaFree(dEnergy));
+    OPM_GPU_SAFE_CALL(cudaFree(dConductivity));
+
+    // -- Compare against the CPU manager.
+    for (std::size_t i = 0; i < numCells; ++i) {
+        DummyFluidState<Scalar> fs;
+        fs.temperature_ = hostTemperatures[i];
+        fs.gasSaturation_ = hostGasSaturations[i];
+
+        const auto& cpuSolidParams =
+            cpuManager.solidEnergyLawParams(static_cast<unsigned>(i));
+        const Scalar cpuEnergy =
+            Opm::EclSpecrockLaw<Scalar>::solidInternalEnergy(cpuSolidParams, fs);
+        BOOST_CHECK_CLOSE(cpuEnergy, hostEnergyResult[i], 1e-9);
+
+        const auto& cpuConductionParams =
+            cpuManager.thermalConductionLawParams(static_cast<unsigned>(i));
+        const Scalar cpuConductivity =
+            GpuThconrLaw::thermalConductivity(cpuConductionParams, fs);
+        BOOST_CHECK_CLOSE(cpuConductivity, hostConductivityResult[i], 1e-9);
+    }
+}