From 8bd403ed5676a93c69d082d2ae2d391f6f11953c Mon Sep 17 00:00:00 2001
From: jakobtorben <Jakob.Torben@sintef.no>
Date: Thu, 26 Jun 2025 17:41:06 +0200
Subject: [PATCH 1/8] Add storage cache validty management to
 fvbasediscretization

---
 .../common/fvbasediscretization.hh            | 74 +++++++++++++++++++
 1 file changed, 74 insertions(+)

diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh
index 37b2c0b95e..03f43c5f69 100644
--- a/opm/models/discretization/common/fvbasediscretization.hh
+++ b/opm/models/discretization/common/fvbasediscretization.hh
@@ -839,6 +839,7 @@ public:
     const EqVector& cachedStorage(unsigned globalIdx, unsigned timeIdx) const
     {
         assert(enableStorageCache_);
+        assert(storageCacheUpToDate_[timeIdx][globalIdx] != 0);
         return storageCache_[timeIdx][globalIdx];
     }
 
@@ -857,6 +858,72 @@ public:
     {
         assert(enableStorageCache_);
         storageCache_[timeIdx][globalIdx] = value;
+        storageCacheUpToDate_[timeIdx][globalIdx] = 1;
+    }
+
+    /*!
+     * \brief Returns true if the storage cache entry for a given DOF and time index is up to date.
+     *
+     * \param globalIdx The global space index for the entity
+     * \param timeIdx The index used by the time discretization
+     */
+    bool storageCacheIsUpToDate(unsigned globalIdx, unsigned timeIdx) const
+    {
+        if (!enableStorageCache_) {
+            return false;
+        }
+        return storageCacheUpToDate_[timeIdx][globalIdx] != 0;
+    }
+
+    /*!
+     * \brief Invalidate the storage cache for a given DOF and time index.
+     *
+     * \param globalIdx The global space index for the entity
+     * \param timeIdx The index used by the time discretization
+     */
+    void invalidateStorageCacheEntry(unsigned globalIdx, unsigned timeIdx) const
+    {
+        if (enableStorageCache_) {
+            storageCacheUpToDate_[timeIdx][globalIdx] = 0;
+        }
+    }
+
+    /*!
+     * \brief Invalidate the whole storage cache for a given time index.
+     *
+     * \param timeIdx The index used by the time discretization
+     */
+    void invalidateStorageCache(unsigned timeIdx) const
+    {
+        if (enableStorageCache_) {
+            std::fill(storageCacheUpToDate_[timeIdx].begin(),
+                      storageCacheUpToDate_[timeIdx].end(),
+                      /*value=*/0);
+        }
+    }
+
+    /*!
+     * \brief Shift storage cache by a given number of time step slots.
+     *
+     * This should be called at the end of each time step to move the storage cache
+     * for timeIdx=0 (current time) to timeIdx=1 (previous time) and so on.
+     *
+     * This method should only be called by the time discretization.
+     *
+     * \param numSlots The number of time step slots for which the
+     *                 storage cache should be shifted.
+     */
+    void shiftStorageCache(unsigned numSlots = 1) const
+    {
+        // If we cannot recycle first iteration storage, it does not make sense to shift the storage cache.
+        if (enableStorageCache_ && !simulator_.problem().recycleFirstIterationStorage()) {
+            for (unsigned timeIdx = 0; timeIdx < historySize - numSlots; ++timeIdx) {
+                storageCache_[timeIdx + numSlots] = storageCache_[timeIdx];
+                storageCacheUpToDate_[timeIdx + numSlots] = storageCacheUpToDate_[timeIdx];
+            }
+
+             // should we invalidate the cache for the most recent time indices? (see shiftIntensiveQuantityCache)
+        }
     }
 
     /*!
@@ -1407,6 +1474,9 @@ public:
         // make the current solution the previous one.
         solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0);
 
+        // shift the storage cache by one position in the history
+        asImp_().shiftStorageCache(/*numSlots=*/1);
+
         // shift the intensive quantities cache by one position in the
         // history
         asImp_().shiftIntensiveQuantityCache(/*numSlots=*/1);
@@ -1840,6 +1910,7 @@ protected:
             const std::size_t numDof = asImp_().numGridDof();
             for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
                 storageCache_[timeIdx].resize(numDof);
+                storageCacheUpToDate_[timeIdx].resize(numDof, /*value=*/0);
             }
         }
 
@@ -1936,6 +2007,9 @@ protected:
 
     mutable std::array<GlobalEqVector, historySize> storageCache_;
 
+    // while these are logically bools, concurrent writes to vector<bool> are not thread safe.
+    mutable std::array<std::vector<unsigned char>, historySize> storageCacheUpToDate_;
+
     bool enableGridAdaptation_;
     bool enableIntensiveQuantityCache_;
     bool enableStorageCache_;

From a9708684332d392076e0bebd47699749fafb1850 Mon Sep 17 00:00:00 2001
From: jakobtorben <Jakob.Torben@sintef.no>
Date: Thu, 26 Jun 2025 17:57:39 +0200
Subject: [PATCH 2/8] Remove the need for having an intensive quantities cache
 for time idx 1

---
 .../common/fvbasediscretization.hh            | 25 +++++++++++++------
 opm/simulators/flow/FIBlackoilModel.hpp       |  2 ++
 2 files changed, 19 insertions(+), 8 deletions(-)

diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh
index 03f43c5f69..9a1c8cbcbc 100644
--- a/opm/models/discretization/common/fvbasediscretization.hh
+++ b/opm/models/discretization/common/fvbasediscretization.hh
@@ -330,6 +330,7 @@ class FvBaseDiscretization
     enum {
         numEq = getPropValue<TypeTag, Properties::NumEq>(),
         historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>(),
+        intensiveCacheSize = 1,
     };
 
     using IntensiveQuantitiesVector = std::vector<IntensiveQuantities,
@@ -413,14 +414,16 @@ public:
 
         PrimaryVariables::init();
         const std::size_t numDof = asImp_().numGridDof();
-        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
+        for (unsigned timeIdx = 0; timeIdx < intensiveCacheSize; ++timeIdx) {
             if (storeIntensiveQuantities()) {
                 intensiveQuantityCache_[timeIdx].resize(numDof);
                 intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof, /*value=*/false);
             }
-
+        }
+        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
             if (enableStorageCache_) {
                 storageCache_[timeIdx].resize(numDof);
+                storageCacheUpToDate_[timeIdx].resize(numDof, /*value=*/false);
             }
         }
 
@@ -525,7 +528,7 @@ public:
         resizeAndResetIntensiveQuantitiesCache_();
         if (storeIntensiveQuantities()) {
             // invalidate all cached intensive quantities
-            for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
+            for (unsigned timeIdx = 0; timeIdx < intensiveCacheSize; ++timeIdx) {
                 invalidateIntensiveQuantitiesCache(timeIdx);
             }
         }
@@ -658,6 +661,8 @@ public:
             return nullptr;
         }
 
+        assert(timeIdx == 0);
+
         // With the storage cache enabled, usually only the
         // intensive quantities for the most recent time step are
         // cached. However, this may be false for some Problem
@@ -686,6 +691,8 @@ public:
             return;
         }
 
+        assert(timeIdx == 0);
+
         intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
         intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = 1;
     }
@@ -715,6 +722,8 @@ public:
      */
     void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
     {
+        assert(timeIdx == 0);
+
         if (storeIntensiveQuantities()) {
             std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
                       intensiveQuantityCacheUpToDate_[timeIdx].end(),
@@ -782,6 +791,7 @@ public:
      */
     void shiftIntensiveQuantityCache(unsigned numSlots = 1)
     {
+        // TODO: Consider removing this method as it i not needed when we only have one time step in the cache.
         if (!storeIntensiveQuantities()) {
             return;
         }
@@ -797,8 +807,7 @@ public:
         }
 
         assert(numSlots > 0);
-
-        for (unsigned timeIdx = 0; timeIdx < historySize - numSlots; ++timeIdx) {
+        for (unsigned timeIdx = 0; timeIdx < intensiveCacheSize - numSlots; ++timeIdx) {
             intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
             intensiveQuantityCacheUpToDate_[timeIdx + numSlots] = intensiveQuantityCacheUpToDate_[timeIdx];
         }
@@ -1917,7 +1926,7 @@ protected:
         // allocate the intensive quantities cache
         if (storeIntensiveQuantities()) {
             const std::size_t numDof = asImp_().numGridDof();
-            for(unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
+            for(unsigned timeIdx = 0; timeIdx < intensiveCacheSize; ++timeIdx) {
                 intensiveQuantityCache_[timeIdx].resize(numDof);
                 intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof);
                 invalidateIntensiveQuantitiesCache(timeIdx);
@@ -1992,10 +2001,10 @@ protected:
 
     // cur is the current iterative solution, prev the converged
     // solution of the previous time step
-    mutable std::array<IntensiveQuantitiesVector, historySize> intensiveQuantityCache_;
+    mutable std::array<IntensiveQuantitiesVector, intensiveCacheSize> intensiveQuantityCache_;
 
     // while these are logically bools, concurrent writes to vector<bool> are not thread safe.
-    mutable std::array<std::vector<unsigned char>, historySize> intensiveQuantityCacheUpToDate_;
+    mutable std::array<std::vector<unsigned char>, intensiveCacheSize> intensiveQuantityCacheUpToDate_;
 
     mutable std::array<std::unique_ptr<DiscreteFunction>, historySize> solution_;
 
diff --git a/opm/simulators/flow/FIBlackoilModel.hpp b/opm/simulators/flow/FIBlackoilModel.hpp
index 501a800d0a..2f3d8041b1 100644
--- a/opm/simulators/flow/FIBlackoilModel.hpp
+++ b/opm/simulators/flow/FIBlackoilModel.hpp
@@ -187,6 +187,8 @@ public:
                       "Run without intensive quantites not enabled: "
                       "Use --enable-intensive-quantity=true");
         }
+
+        assert(timeIdx == 0);
         const auto* intquant = this->cachedIntensiveQuantities(globalIdx, timeIdx);
         if (!intquant) {
             OPM_THROW(std::logic_error, "Intensive quantites need to be updated in code");

From a49ec4e429d15e12d83b552d2a01e661fdb4fc30 Mon Sep 17 00:00:00 2001
From: jakobtorben <Jakob.Torben@sintef.no>
Date: Fri, 27 Jun 2025 12:59:19 +0200
Subject: [PATCH 3/8] Introduces a new TypeTag to manage the history size for
 intensiveQuantityCache

---
 .../common/fvbasediscretization.hh            | 28 +++++++++++--------
 .../discretization/common/fvbaseproperties.hh |  6 ++++
 .../flow/FlowBaseProblemProperties.hpp        |  4 +++
 3 files changed, 27 insertions(+), 11 deletions(-)

diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh
index 9a1c8cbcbc..14e7cda4d7 100644
--- a/opm/models/discretization/common/fvbasediscretization.hh
+++ b/opm/models/discretization/common/fvbasediscretization.hh
@@ -251,6 +251,10 @@ template<class TypeTag>
 struct TimeDiscHistorySize<TypeTag, TTag::FvBaseDiscretization>
 { static constexpr int value = 2; };
 
+template<class TypeTag>
+struct IntensiveQuantityHistorySize<TypeTag, TTag::FvBaseDiscretization>
+{ static constexpr int value = 2; };
+
 //! Most models use extensive quantities for their storage term (so far, only the Stokes
 //! model does), so we disable this by default.
 template<class TypeTag>
@@ -330,7 +334,7 @@ class FvBaseDiscretization
     enum {
         numEq = getPropValue<TypeTag, Properties::NumEq>(),
         historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>(),
-        intensiveCacheSize = 1,
+        intensiveQuantityHistorySize = getPropValue<TypeTag, Properties::IntensiveQuantityHistorySize>(),
     };
 
     using IntensiveQuantitiesVector = std::vector<IntensiveQuantities,
@@ -414,7 +418,7 @@ public:
 
         PrimaryVariables::init();
         const std::size_t numDof = asImp_().numGridDof();
-        for (unsigned timeIdx = 0; timeIdx < intensiveCacheSize; ++timeIdx) {
+        for (unsigned timeIdx = 0; timeIdx < intensiveQuantityHistorySize; ++timeIdx) {
             if (storeIntensiveQuantities()) {
                 intensiveQuantityCache_[timeIdx].resize(numDof);
                 intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof, /*value=*/false);
@@ -528,7 +532,7 @@ public:
         resizeAndResetIntensiveQuantitiesCache_();
         if (storeIntensiveQuantities()) {
             // invalidate all cached intensive quantities
-            for (unsigned timeIdx = 0; timeIdx < intensiveCacheSize; ++timeIdx) {
+            for (unsigned timeIdx = 0; timeIdx < intensiveQuantityHistorySize; ++timeIdx) {
                 invalidateIntensiveQuantitiesCache(timeIdx);
             }
         }
@@ -661,7 +665,7 @@ public:
             return nullptr;
         }
 
-        assert(timeIdx == 0);
+        assert(timeIdx < intensiveQuantityHistorySize);
 
         // With the storage cache enabled, usually only the
         // intensive quantities for the most recent time step are
@@ -691,10 +695,10 @@ public:
             return;
         }
 
-        assert(timeIdx == 0);
+        assert(timeIdx < intensiveQuantityHistorySize);
 
         intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
-        intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = 1;
+        intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = true;
     }
 
     /*!
@@ -722,7 +726,7 @@ public:
      */
     void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
     {
-        assert(timeIdx == 0);
+        assert(timeIdx < intensiveQuantityHistorySize);
 
         if (storeIntensiveQuantities()) {
             std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
@@ -807,7 +811,7 @@ public:
         }
 
         assert(numSlots > 0);
-        for (unsigned timeIdx = 0; timeIdx < intensiveCacheSize - numSlots; ++timeIdx) {
+        for (unsigned timeIdx = 0; timeIdx < intensiveQuantityHistorySize - numSlots; ++timeIdx) {
             intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
             intensiveQuantityCacheUpToDate_[timeIdx + numSlots] = intensiveQuantityCacheUpToDate_[timeIdx];
         }
@@ -848,6 +852,7 @@ public:
     const EqVector& cachedStorage(unsigned globalIdx, unsigned timeIdx) const
     {
         assert(enableStorageCache_);
+        assert(timeIdx < historySize);
         assert(storageCacheUpToDate_[timeIdx][globalIdx] != 0);
         return storageCache_[timeIdx][globalIdx];
     }
@@ -866,6 +871,7 @@ public:
     void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector& value) const
     {
         assert(enableStorageCache_);
+        assert(timeIdx < historySize);
         storageCache_[timeIdx][globalIdx] = value;
         storageCacheUpToDate_[timeIdx][globalIdx] = 1;
     }
@@ -1926,7 +1932,7 @@ protected:
         // allocate the intensive quantities cache
         if (storeIntensiveQuantities()) {
             const std::size_t numDof = asImp_().numGridDof();
-            for(unsigned timeIdx = 0; timeIdx < intensiveCacheSize; ++timeIdx) {
+            for(unsigned timeIdx = 0; timeIdx < intensiveQuantityHistorySize; ++timeIdx) {
                 intensiveQuantityCache_[timeIdx].resize(numDof);
                 intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof);
                 invalidateIntensiveQuantitiesCache(timeIdx);
@@ -2001,10 +2007,10 @@ protected:
 
     // cur is the current iterative solution, prev the converged
     // solution of the previous time step
-    mutable std::array<IntensiveQuantitiesVector, intensiveCacheSize> intensiveQuantityCache_;
+    mutable std::array<IntensiveQuantitiesVector, intensiveQuantityHistorySize> intensiveQuantityCache_;
 
     // while these are logically bools, concurrent writes to vector<bool> are not thread safe.
-    mutable std::array<std::vector<unsigned char>, intensiveCacheSize> intensiveQuantityCacheUpToDate_;
+    mutable std::array<std::vector<unsigned char>, intensiveQuantityHistorySize> intensiveQuantityCacheUpToDate_;
 
     mutable std::array<std::unique_ptr<DiscreteFunction>, historySize> solution_;
 
diff --git a/opm/models/discretization/common/fvbaseproperties.hh b/opm/models/discretization/common/fvbaseproperties.hh
index bcd773a3cd..2187ffadf0 100644
--- a/opm/models/discretization/common/fvbaseproperties.hh
+++ b/opm/models/discretization/common/fvbaseproperties.hh
@@ -224,6 +224,12 @@ struct DofMapper { using type = UndefinedProperty; };
 template<class TypeTag, class MyTypeTag>
 struct TimeDiscHistorySize { using type = UndefinedProperty; };
 
+/*!
+ * \brief The history size required by the intensive quantities
+ */
+template<class TypeTag, class MyTypeTag>
+struct IntensiveQuantityHistorySize { using type = UndefinedProperty; };
+
 /*!
  * \brief Specify whether the storage terms use extensive quantities or not.
  *
diff --git a/opm/simulators/flow/FlowBaseProblemProperties.hpp b/opm/simulators/flow/FlowBaseProblemProperties.hpp
index 4b7dd4066a..352c2260d6 100644
--- a/opm/simulators/flow/FlowBaseProblemProperties.hpp
+++ b/opm/simulators/flow/FlowBaseProblemProperties.hpp
@@ -241,6 +241,10 @@ template<class TypeTag>
 struct EnableDebuggingChecks<TypeTag, TTag::FlowBaseProblem>
 { static constexpr bool value = true; };
 
+template<class TypeTag>
+struct IntensiveQuantityHistorySize<TypeTag, TTag::FlowBaseProblem>
+{ static constexpr int value = 1; };
+
 } // namespace Opm::Properties
 
 #endif // OPM_BASE_FLOW_PROBLEM_PROPERTIES_HPP

From 32cb7887a5a7d6779263f7d3b19e1043d4ddeced Mon Sep 17 00:00:00 2001
From: jakobtorben <Jakob.Torben@sintef.no>
Date: Wed, 2 Jul 2025 16:00:21 +0200
Subject: [PATCH 4/8] Make intensiveQuantityHistorySize a runtime controlled
 check

---
 .../common/fvbasediscretization.hh            | 44 +++++++------------
 .../discretization/common/fvbaseproblem.hh    | 17 ++++++-
 .../discretization/common/fvbaseproperties.hh |  5 ---
 opm/simulators/flow/FIBlackoilModel.hpp       |  2 +-
 .../flow/FlowBaseProblemProperties.hpp        |  3 --
 opm/simulators/flow/FlowProblemBlackoil.hpp   |  2 +-
 6 files changed, 35 insertions(+), 38 deletions(-)

diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh
index 14e7cda4d7..573889e542 100644
--- a/opm/models/discretization/common/fvbasediscretization.hh
+++ b/opm/models/discretization/common/fvbasediscretization.hh
@@ -251,9 +251,6 @@ template<class TypeTag>
 struct TimeDiscHistorySize<TypeTag, TTag::FvBaseDiscretization>
 { static constexpr int value = 2; };
 
-template<class TypeTag>
-struct IntensiveQuantityHistorySize<TypeTag, TTag::FvBaseDiscretization>
-{ static constexpr int value = 2; };
 
 //! Most models use extensive quantities for their storage term (so far, only the Stokes
 //! model does), so we disable this by default.
@@ -334,7 +331,6 @@ class FvBaseDiscretization
     enum {
         numEq = getPropValue<TypeTag, Properties::NumEq>(),
         historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>(),
-        intensiveQuantityHistorySize = getPropValue<TypeTag, Properties::IntensiveQuantityHistorySize>(),
     };
 
     using IntensiveQuantitiesVector = std::vector<IntensiveQuantities,
@@ -417,21 +413,7 @@ public:
         }
 
         PrimaryVariables::init();
-        const std::size_t numDof = asImp_().numGridDof();
-        for (unsigned timeIdx = 0; timeIdx < intensiveQuantityHistorySize; ++timeIdx) {
-            if (storeIntensiveQuantities()) {
-                intensiveQuantityCache_[timeIdx].resize(numDof);
-                intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof, /*value=*/false);
-            }
-        }
-        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
-            if (enableStorageCache_) {
-                storageCache_[timeIdx].resize(numDof);
-                storageCacheUpToDate_[timeIdx].resize(numDof, /*value=*/false);
-            }
-        }
-
-        resizeAndResetIntensiveQuantitiesCache_();
+        // Setting up the intensive quantities cache and storage cache is done in finishInit()
         asImp_().registerOutputModules_();
     }
 
@@ -532,7 +514,8 @@ public:
         resizeAndResetIntensiveQuantitiesCache_();
         if (storeIntensiveQuantities()) {
             // invalidate all cached intensive quantities
-            for (unsigned timeIdx = 0; timeIdx < intensiveQuantityHistorySize; ++timeIdx) {
+            const unsigned intensiveHistorySize = simulator_.problem().intensiveQuantityHistorySize();
+            for (unsigned timeIdx = 0; timeIdx < intensiveHistorySize; ++timeIdx) {
                 invalidateIntensiveQuantitiesCache(timeIdx);
             }
         }
@@ -665,7 +648,7 @@ public:
             return nullptr;
         }
 
-        assert(timeIdx < intensiveQuantityHistorySize);
+        assert(timeIdx < simulator_.problem().intensiveQuantityHistorySize());
 
         // With the storage cache enabled, usually only the
         // intensive quantities for the most recent time step are
@@ -695,7 +678,7 @@ public:
             return;
         }
 
-        assert(timeIdx < intensiveQuantityHistorySize);
+        assert(timeIdx < simulator_.problem().intensiveQuantityHistorySize());
 
         intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
         intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = true;
@@ -726,7 +709,7 @@ public:
      */
     void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
     {
-        assert(timeIdx < intensiveQuantityHistorySize);
+        assert(timeIdx < simulator_.problem().intensiveQuantityHistorySize());
 
         if (storeIntensiveQuantities()) {
             std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
@@ -811,7 +794,8 @@ public:
         }
 
         assert(numSlots > 0);
-        for (unsigned timeIdx = 0; timeIdx < intensiveQuantityHistorySize - numSlots; ++timeIdx) {
+        const unsigned intensiveHistorySize = simulator_.problem().intensiveQuantityHistorySize();
+        for (unsigned timeIdx = 0; timeIdx < intensiveHistorySize - numSlots; ++timeIdx) {
             intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
             intensiveQuantityCacheUpToDate_[timeIdx + numSlots] = intensiveQuantityCacheUpToDate_[timeIdx];
         }
@@ -1932,7 +1916,13 @@ protected:
         // allocate the intensive quantities cache
         if (storeIntensiveQuantities()) {
             const std::size_t numDof = asImp_().numGridDof();
-            for(unsigned timeIdx = 0; timeIdx < intensiveQuantityHistorySize; ++timeIdx) {
+            const unsigned intensiveHistorySize = simulator_.problem().intensiveQuantityHistorySize();
+
+            // resize the vectors based on runtime history size
+            intensiveQuantityCache_.resize(intensiveHistorySize);
+            intensiveQuantityCacheUpToDate_.resize(intensiveHistorySize);
+
+            for(unsigned timeIdx = 0; timeIdx < intensiveHistorySize; ++timeIdx) {
                 intensiveQuantityCache_[timeIdx].resize(numDof);
                 intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof);
                 invalidateIntensiveQuantitiesCache(timeIdx);
@@ -2007,10 +1997,10 @@ protected:
 
     // cur is the current iterative solution, prev the converged
     // solution of the previous time step
-    mutable std::array<IntensiveQuantitiesVector, intensiveQuantityHistorySize> intensiveQuantityCache_;
+    mutable std::vector<IntensiveQuantitiesVector> intensiveQuantityCache_;
 
     // while these are logically bools, concurrent writes to vector<bool> are not thread safe.
-    mutable std::array<std::vector<unsigned char>, intensiveQuantityHistorySize> intensiveQuantityCacheUpToDate_;
+    mutable std::vector<std::vector<unsigned char>> intensiveQuantityCacheUpToDate_;
 
     mutable std::array<std::unique_ptr<DiscreteFunction>, historySize> solution_;
 
diff --git a/opm/models/discretization/common/fvbaseproblem.hh b/opm/models/discretization/common/fvbaseproblem.hh
index 9c999dd347..1767135c0a 100644
--- a/opm/models/discretization/common/fvbaseproblem.hh
+++ b/opm/models/discretization/common/fvbaseproblem.hh
@@ -189,9 +189,24 @@ public:
      * This is only relevant if the storage cache is enabled and is usually the case,
      * i.e., this method only needs to be overwritten in rare corner cases.
      */
-    bool recycleFirstIterationStorage() const
+    virtual bool recycleFirstIterationStorage() const
     { return true; }
 
+    /*!
+     * \brief Returns the required history size for intensive quantities cache.
+     *
+     * If recycleFirstIterationStorage() is true and storage cache is enabled, only the current time step intensive
+     * quantities need to be cached (size = 1). Otherwise, both current and previous
+     * time step intensive quantities are needed (size = 2).
+     */
+    unsigned intensiveQuantityHistorySize() const
+    {
+        if (Parameters::Get<Parameters::EnableStorageCache>() && recycleFirstIterationStorage()) {
+            return 1;
+        }
+        return 2;
+    }
+
     /*!
      * \brief Determine the directory for simulation output.
      *
diff --git a/opm/models/discretization/common/fvbaseproperties.hh b/opm/models/discretization/common/fvbaseproperties.hh
index 2187ffadf0..9075d7242f 100644
--- a/opm/models/discretization/common/fvbaseproperties.hh
+++ b/opm/models/discretization/common/fvbaseproperties.hh
@@ -224,11 +224,6 @@ struct DofMapper { using type = UndefinedProperty; };
 template<class TypeTag, class MyTypeTag>
 struct TimeDiscHistorySize { using type = UndefinedProperty; };
 
-/*!
- * \brief The history size required by the intensive quantities
- */
-template<class TypeTag, class MyTypeTag>
-struct IntensiveQuantityHistorySize { using type = UndefinedProperty; };
 
 /*!
  * \brief Specify whether the storage terms use extensive quantities or not.
diff --git a/opm/simulators/flow/FIBlackoilModel.hpp b/opm/simulators/flow/FIBlackoilModel.hpp
index 2f3d8041b1..8b18f71655 100644
--- a/opm/simulators/flow/FIBlackoilModel.hpp
+++ b/opm/simulators/flow/FIBlackoilModel.hpp
@@ -188,7 +188,7 @@ public:
                       "Use --enable-intensive-quantity=true");
         }
 
-        assert(timeIdx == 0);
+        assert(timeIdx < this->simulator_.problem().intensiveQuantityHistorySize());
         const auto* intquant = this->cachedIntensiveQuantities(globalIdx, timeIdx);
         if (!intquant) {
             OPM_THROW(std::logic_error, "Intensive quantites need to be updated in code");
diff --git a/opm/simulators/flow/FlowBaseProblemProperties.hpp b/opm/simulators/flow/FlowBaseProblemProperties.hpp
index 352c2260d6..3d88742a04 100644
--- a/opm/simulators/flow/FlowBaseProblemProperties.hpp
+++ b/opm/simulators/flow/FlowBaseProblemProperties.hpp
@@ -241,9 +241,6 @@ template<class TypeTag>
 struct EnableDebuggingChecks<TypeTag, TTag::FlowBaseProblem>
 { static constexpr bool value = true; };
 
-template<class TypeTag>
-struct IntensiveQuantityHistorySize<TypeTag, TTag::FlowBaseProblem>
-{ static constexpr int value = 1; };
 
 } // namespace Opm::Properties
 
diff --git a/opm/simulators/flow/FlowProblemBlackoil.hpp b/opm/simulators/flow/FlowProblemBlackoil.hpp
index 0052215d78..a2d8ed90f0 100644
--- a/opm/simulators/flow/FlowProblemBlackoil.hpp
+++ b/opm/simulators/flow/FlowProblemBlackoil.hpp
@@ -866,7 +866,7 @@ public:
      * or DRVDT are active. Nor if the porosity is changes between timesteps
      * using a pore volume multiplier (i.e., poreVolumeMultiplier() != 1.0)
      */
-    bool recycleFirstIterationStorage() const
+    bool recycleFirstIterationStorage() const override
     {
         int episodeIdx = this->episodeIndex();
         return !this->mixControls_.drsdtActive(episodeIdx) &&

From 5b94c000fadf5e7dc4ff851bcc04d90e0b66ef27 Mon Sep 17 00:00:00 2001
From: jakobtorben <Jakob.Torben@sintef.no>
Date: Fri, 4 Jul 2025 11:46:08 +0200
Subject: [PATCH 5/8] Move setup of intensiveQuantities after all
 problem-specific setup

---
 .../discretization/common/fvbasediscretization.hh    | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh
index 573889e542..ebe017f2f4 100644
--- a/opm/models/discretization/common/fvbasediscretization.hh
+++ b/opm/models/discretization/common/fvbasediscretization.hh
@@ -512,13 +512,6 @@ public:
         }
 
         resizeAndResetIntensiveQuantitiesCache_();
-        if (storeIntensiveQuantities()) {
-            // invalidate all cached intensive quantities
-            const unsigned intensiveHistorySize = simulator_.problem().intensiveQuantityHistorySize();
-            for (unsigned timeIdx = 0; timeIdx < intensiveHistorySize; ++timeIdx) {
-                invalidateIntensiveQuantitiesCache(timeIdx);
-            }
-        }
 
         newtonMethod_.finishInit();
     }
@@ -573,6 +566,11 @@ public:
             solution(timeIdx) = solution(/*timeIdx=*/0);
         }
 
+        // Initialize intensive quantities cache now that all problem-specific parameters are available.
+        // This ensures intensiveQuantityHistorySize is correct based on recycleFirstIterationStorage().
+        // TODO: Where this is done should perhaps be changed once finishInit() is refactored.
+        resizeAndResetIntensiveQuantitiesCache_();
+
         simulator_.problem().initialSolutionApplied();
 
 #ifndef NDEBUG

From 8f6b0fc893551702341ae6213e3ea6af9f5525a6 Mon Sep 17 00:00:00 2001
From: jakobtorben <Jakob.Torben@sintef.no>
Date: Fri, 4 Jul 2025 14:36:14 +0200
Subject: [PATCH 6/8] Make intensive quantities update logic in
 FvBaseElementContext cache size aware

---
 .../common/fvbaseelementcontext.hh            | 29 ++++++++++++-------
 1 file changed, 19 insertions(+), 10 deletions(-)

diff --git a/opm/models/discretization/common/fvbaseelementcontext.hh b/opm/models/discretization/common/fvbaseelementcontext.hh
index 3ad54b6d71..7733c30be5 100644
--- a/opm/models/discretization/common/fvbaseelementcontext.hh
+++ b/opm/models/discretization/common/fvbaseelementcontext.hh
@@ -552,24 +552,33 @@ protected:
         // update the intensive quantities for the whole history
         const SolutionVector& globalSol = model().solution(timeIdx);
 
+        const unsigned intensiveHistorySize = problem().intensiveQuantityHistorySize();
+
         // update the non-gradient quantities
         for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx) {
             unsigned globalIdx = globalSpaceIndex(dofIdx, timeIdx);
             const PrimaryVariables& dofSol = globalSol[globalIdx];
             dofVars_[dofIdx].priVars[timeIdx] = &dofSol;
 
-            dofVars_[dofIdx].thermodynamicHint[timeIdx] =
-                model().thermodynamicHint(globalIdx, timeIdx);
-
-            const auto *cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx);
-            if (cachedIntQuants) {
-                dofVars_[dofIdx].intensiveQuantities[timeIdx] = *cachedIntQuants;
+            if (timeIdx >= intensiveHistorySize) {
+                // If we do not have a cache for this time index, it does not make sense to try to fetch it or
+                // update the intensive quantities cache for this time index.
+                updateSingleIntQuants_(dofSol, dofIdx, timeIdx);
             }
             else {
-                updateSingleIntQuants_(dofSol, dofIdx, timeIdx);
-                model().updateCachedIntensiveQuantities(dofVars_[dofIdx].intensiveQuantities[timeIdx],
-                                                        globalIdx,
-                                                        timeIdx);
+                // If we have a cache for this time index, we can fetch the intensive quantities from the cache.
+                // We can then also fetch the thermodynamic hint from the cache.
+                dofVars_[dofIdx].thermodynamicHint[timeIdx] = model().thermodynamicHint(globalIdx, timeIdx);
+                const auto *cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx);
+                if (cachedIntQuants) {
+                    dofVars_[dofIdx].intensiveQuantities[timeIdx] = *cachedIntQuants;
+                }
+                else {
+                    updateSingleIntQuants_(dofSol, dofIdx, timeIdx);
+                    model().updateCachedIntensiveQuantities(dofVars_[dofIdx].intensiveQuantities[timeIdx],
+                                                            globalIdx,
+                                                            timeIdx);
+                }
             }
         }
     }

From 282f8606bc4ecb711c33f2cfbf88cd073929c287 Mon Sep 17 00:00:00 2001
From: jakobtorben <Jakob.Torben@sintef.no>
Date: Mon, 4 Aug 2025 17:04:50 +0200
Subject: [PATCH 7/8] Cache intensive quantity history size calculation

---
 .../common/fvbasediscretization.hh            | 20 ++++++++++++++-----
 .../common/fvbaseelementcontext.hh            |  4 ++--
 .../discretization/common/fvbaseproblem.hh    | 14 ++++++-------
 .../discretization/common/fvbaseproperties.hh |  1 -
 opm/simulators/flow/FIBlackoilModel.hpp       |  2 +-
 .../flow/FlowBaseProblemProperties.hpp        |  1 -
 6 files changed, 25 insertions(+), 17 deletions(-)

diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh
index ebe017f2f4..57ba14c680 100644
--- a/opm/models/discretization/common/fvbasediscretization.hh
+++ b/opm/models/discretization/common/fvbasediscretization.hh
@@ -404,6 +404,7 @@ public:
         , enableIntensiveQuantityCache_(Parameters::Get<Parameters::EnableIntensiveQuantityCache>())
         , enableStorageCache_(Parameters::Get<Parameters::EnableStorageCache>())
         , enableThermodynamicHints_(Parameters::Get<Parameters::EnableThermodynamicHints>())
+        , cachedIntensiveQuantityHistorySize_(-1)
     {
         const bool isEcfv = std::is_same_v<Discretization, EcfvDiscretization<TypeTag>>;
         if (enableGridAdaptation_ && !isEcfv) {
@@ -646,7 +647,7 @@ public:
             return nullptr;
         }
 
-        assert(timeIdx < simulator_.problem().intensiveQuantityHistorySize());
+        assert(timeIdx < cachedIntensiveQuantityHistorySize_);
 
         // With the storage cache enabled, usually only the
         // intensive quantities for the most recent time step are
@@ -676,7 +677,7 @@ public:
             return;
         }
 
-        assert(timeIdx < simulator_.problem().intensiveQuantityHistorySize());
+        assert(timeIdx < cachedIntensiveQuantityHistorySize_);
 
         intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
         intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = true;
@@ -700,6 +701,12 @@ public:
         intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = newValue ? 1 : 0;
     }
 
+    /*!
+     * \brief Get the cached intensive quantity history size.
+     */
+    unsigned cachedIntensiveQuantityHistorySize() const
+    { return cachedIntensiveQuantityHistorySize_; }
+
     /*!
      * \brief Invalidate the whole intensive quantity cache for time index.
      *
@@ -707,7 +714,7 @@ public:
      */
     void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
     {
-        assert(timeIdx < simulator_.problem().intensiveQuantityHistorySize());
+        assert(timeIdx < cachedIntensiveQuantityHistorySize_);
 
         if (storeIntensiveQuantities()) {
             std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
@@ -792,7 +799,7 @@ public:
         }
 
         assert(numSlots > 0);
-        const unsigned intensiveHistorySize = simulator_.problem().intensiveQuantityHistorySize();
+        const unsigned intensiveHistorySize = cachedIntensiveQuantityHistorySize_;
         for (unsigned timeIdx = 0; timeIdx < intensiveHistorySize - numSlots; ++timeIdx) {
             intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
             intensiveQuantityCacheUpToDate_[timeIdx + numSlots] = intensiveQuantityCacheUpToDate_[timeIdx];
@@ -1914,7 +1921,8 @@ protected:
         // allocate the intensive quantities cache
         if (storeIntensiveQuantities()) {
             const std::size_t numDof = asImp_().numGridDof();
-            const unsigned intensiveHistorySize = simulator_.problem().intensiveQuantityHistorySize();
+            cachedIntensiveQuantityHistorySize_ = simulator_.problem().intensiveQuantityHistorySize();
+            const unsigned intensiveHistorySize = cachedIntensiveQuantityHistorySize_;
 
             // resize the vectors based on runtime history size
             intensiveQuantityCache_.resize(intensiveHistorySize);
@@ -2017,6 +2025,8 @@ protected:
     bool enableIntensiveQuantityCache_;
     bool enableStorageCache_;
     bool enableThermodynamicHints_;
+
+    mutable unsigned cachedIntensiveQuantityHistorySize_;
 };
 
 /*!
diff --git a/opm/models/discretization/common/fvbaseelementcontext.hh b/opm/models/discretization/common/fvbaseelementcontext.hh
index 7733c30be5..140769bcc8 100644
--- a/opm/models/discretization/common/fvbaseelementcontext.hh
+++ b/opm/models/discretization/common/fvbaseelementcontext.hh
@@ -552,7 +552,7 @@ protected:
         // update the intensive quantities for the whole history
         const SolutionVector& globalSol = model().solution(timeIdx);
 
-        const unsigned intensiveHistorySize = problem().intensiveQuantityHistorySize();
+        const unsigned intensiveHistorySize = model().cachedIntensiveQuantityHistorySize();
 
         // update the non-gradient quantities
         for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx) {
@@ -569,7 +569,7 @@ protected:
                 // If we have a cache for this time index, we can fetch the intensive quantities from the cache.
                 // We can then also fetch the thermodynamic hint from the cache.
                 dofVars_[dofIdx].thermodynamicHint[timeIdx] = model().thermodynamicHint(globalIdx, timeIdx);
-                const auto *cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx);
+                const auto* cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx);
                 if (cachedIntQuants) {
                     dofVars_[dofIdx].intensiveQuantities[timeIdx] = *cachedIntQuants;
                 }
diff --git a/opm/models/discretization/common/fvbaseproblem.hh b/opm/models/discretization/common/fvbaseproblem.hh
index 1767135c0a..78f759c4cc 100644
--- a/opm/models/discretization/common/fvbaseproblem.hh
+++ b/opm/models/discretization/common/fvbaseproblem.hh
@@ -199,13 +199,13 @@ public:
      * quantities need to be cached (size = 1). Otherwise, both current and previous
      * time step intensive quantities are needed (size = 2).
      */
-    unsigned intensiveQuantityHistorySize() const
-    {
-        if (Parameters::Get<Parameters::EnableStorageCache>() && recycleFirstIterationStorage()) {
-            return 1;
-        }
-        return 2;
-    }
+     unsigned intensiveQuantityHistorySize() const
+     {
+         if (Parameters::Get<Parameters::EnableStorageCache>() && recycleFirstIterationStorage()) {
+             return 1;
+         }
+         return 2;
+     }
 
     /*!
      * \brief Determine the directory for simulation output.
diff --git a/opm/models/discretization/common/fvbaseproperties.hh b/opm/models/discretization/common/fvbaseproperties.hh
index 9075d7242f..bcd773a3cd 100644
--- a/opm/models/discretization/common/fvbaseproperties.hh
+++ b/opm/models/discretization/common/fvbaseproperties.hh
@@ -224,7 +224,6 @@ struct DofMapper { using type = UndefinedProperty; };
 template<class TypeTag, class MyTypeTag>
 struct TimeDiscHistorySize { using type = UndefinedProperty; };
 
-
 /*!
  * \brief Specify whether the storage terms use extensive quantities or not.
  *
diff --git a/opm/simulators/flow/FIBlackoilModel.hpp b/opm/simulators/flow/FIBlackoilModel.hpp
index 8b18f71655..3cd9f8c6bb 100644
--- a/opm/simulators/flow/FIBlackoilModel.hpp
+++ b/opm/simulators/flow/FIBlackoilModel.hpp
@@ -188,7 +188,7 @@ public:
                       "Use --enable-intensive-quantity=true");
         }
 
-        assert(timeIdx < this->simulator_.problem().intensiveQuantityHistorySize());
+        assert(timeIdx < this->cachedIntensiveQuantityHistorySize());
         const auto* intquant = this->cachedIntensiveQuantities(globalIdx, timeIdx);
         if (!intquant) {
             OPM_THROW(std::logic_error, "Intensive quantites need to be updated in code");
diff --git a/opm/simulators/flow/FlowBaseProblemProperties.hpp b/opm/simulators/flow/FlowBaseProblemProperties.hpp
index 3d88742a04..4b7dd4066a 100644
--- a/opm/simulators/flow/FlowBaseProblemProperties.hpp
+++ b/opm/simulators/flow/FlowBaseProblemProperties.hpp
@@ -241,7 +241,6 @@ template<class TypeTag>
 struct EnableDebuggingChecks<TypeTag, TTag::FlowBaseProblem>
 { static constexpr bool value = true; };
 
-
 } // namespace Opm::Properties
 
 #endif // OPM_BASE_FLOW_PROBLEM_PROPERTIES_HPP

From bf10320089d9d37b8c7688d45c0dac473ef92be2 Mon Sep 17 00:00:00 2001
From: jakobtorben <Jakob.Torben@sintef.no>
Date: Thu, 7 Aug 2025 12:03:40 +0200
Subject: [PATCH 8/8] Improve cache validity checks for intensive quantities
 and storage

---
 .../common/fvbasediscretization.hh            | 46 +++++++++++--------
 1 file changed, 27 insertions(+), 19 deletions(-)

diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh
index 57ba14c680..6705ec9b40 100644
--- a/opm/models/discretization/common/fvbasediscretization.hh
+++ b/opm/models/discretization/common/fvbasediscretization.hh
@@ -251,7 +251,6 @@ template<class TypeTag>
 struct TimeDiscHistorySize<TypeTag, TTag::FvBaseDiscretization>
 { static constexpr int value = 2; };
 
-
 //! Most models use extensive quantities for their storage term (so far, only the Stokes
 //! model does), so we disable this by default.
 template<class TypeTag>
@@ -415,6 +414,7 @@ public:
 
         PrimaryVariables::init();
         // Setting up the intensive quantities cache and storage cache is done in finishInit()
+        // and applyInitialSolution() to ensure that the history size is correct.
         asImp_().registerOutputModules_();
     }
 
@@ -643,12 +643,12 @@ public:
      */
     const IntensiveQuantities* cachedIntensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
     {
-        if (!enableIntensiveQuantityCache_ || !intensiveQuantityCacheUpToDate_[timeIdx][globalIdx]) {
+        if (!enableIntensiveQuantityCache_ ||
+            timeIdx >= cachedIntensiveQuantityHistorySize_ ||
+            !intensiveQuantityCacheUpToDate_[timeIdx][globalIdx]) {
             return nullptr;
         }
 
-        assert(timeIdx < cachedIntensiveQuantityHistorySize_);
-
         // With the storage cache enabled, usually only the
         // intensive quantities for the most recent time step are
         // cached. However, this may be false for some Problem
@@ -673,12 +673,10 @@ public:
                                          unsigned globalIdx,
                                          unsigned timeIdx) const
     {
-        if (!storeIntensiveQuantities()) {
+        if (!storeIntensiveQuantities() || timeIdx >= cachedIntensiveQuantityHistorySize_) {
             return;
         }
 
-        assert(timeIdx < cachedIntensiveQuantityHistorySize_);
-
         intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
         intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = true;
     }
@@ -714,7 +712,9 @@ public:
      */
     void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
     {
-        assert(timeIdx < cachedIntensiveQuantityHistorySize_);
+        if (timeIdx >= cachedIntensiveQuantityHistorySize_) {
+            return;
+        }
 
         if (storeIntensiveQuantities()) {
             std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
@@ -783,8 +783,7 @@ public:
      */
     void shiftIntensiveQuantityCache(unsigned numSlots = 1)
     {
-        // TODO: Consider removing this method as it i not needed when we only have one time step in the cache.
-        if (!storeIntensiveQuantities()) {
+        if (!storeIntensiveQuantities() || numSlots <= 0) {
             return;
         }
 
@@ -798,7 +797,6 @@ public:
             return;
         }
 
-        assert(numSlots > 0);
         const unsigned intensiveHistorySize = cachedIntensiveQuantityHistorySize_;
         for (unsigned timeIdx = 0; timeIdx < intensiveHistorySize - numSlots; ++timeIdx) {
             intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
@@ -835,14 +833,22 @@ public:
      * volume unit at a given time. The user is responsible for making sure that the
      * value of this is correct and that it can be used before this method is called.
      *
+     * \attention If the storage cache is disabled, or if the entry is not up to date,
+     *            this method will throw a std::logic_error.
+     *
      * \param globalDofIdx The index of the relevant degree of freedom in a grid-global vector
      * \param timeIdx The relevant index for the time discretization
      */
     const EqVector& cachedStorage(unsigned globalIdx, unsigned timeIdx) const
     {
-        assert(enableStorageCache_);
-        assert(timeIdx < historySize);
-        assert(storageCacheUpToDate_[timeIdx][globalIdx] != 0);
+        if (!enableStorageCache_ ||
+            timeIdx >= historySize ||
+            !storageCacheUpToDate_[timeIdx][globalIdx]) {
+            throw std::logic_error("Cached storage is not available or up to date for the requested "
+                                   "global index and time index. Make sure storage cache is enabled "
+                                   "and the entry is valid before calling this method.");
+        }
+
         return storageCache_[timeIdx][globalIdx];
     }
 
@@ -859,8 +865,10 @@ public:
      */
     void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector& value) const
     {
-        assert(enableStorageCache_);
-        assert(timeIdx < historySize);
+        if (!enableStorageCache_ || timeIdx >= historySize) {
+            return;
+        }
+
         storageCache_[timeIdx][globalIdx] = value;
         storageCacheUpToDate_[timeIdx][globalIdx] = 1;
     }
@@ -873,7 +881,7 @@ public:
      */
     bool storageCacheIsUpToDate(unsigned globalIdx, unsigned timeIdx) const
     {
-        if (!enableStorageCache_) {
+        if (!enableStorageCache_ || timeIdx >= historySize) {
             return false;
         }
         return storageCacheUpToDate_[timeIdx][globalIdx] != 0;
@@ -887,7 +895,7 @@ public:
      */
     void invalidateStorageCacheEntry(unsigned globalIdx, unsigned timeIdx) const
     {
-        if (enableStorageCache_) {
+        if (enableStorageCache_ && timeIdx < historySize) {
             storageCacheUpToDate_[timeIdx][globalIdx] = 0;
         }
     }
@@ -899,7 +907,7 @@ public:
      */
     void invalidateStorageCache(unsigned timeIdx) const
     {
-        if (enableStorageCache_) {
+        if (enableStorageCache_ && timeIdx < historySize) {
             std::fill(storageCacheUpToDate_[timeIdx].begin(),
                       storageCacheUpToDate_[timeIdx].end(),
                       /*value=*/0);