Allocate the per-cell material law parameters from an arena

opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.cpp

@@ -85,6 +85,16 @@ run(const IntLookupFunction& fieldPropIntOnLeafAssigner,
     std::vector<std::vector<MaterialLawParams>*> mlpArray;
     initArrays_(satnumArray, imbnumArray, mlpArray);
     const auto num_arrays = mlpArray.size();
+    // One contiguous arena per parameter array. Without this each cell's
+    // nested parameter object is a separate heap allocation reached through
+    // its own shared_ptr control block, so the evaluation loop walks the
+    // parameter set in allocation order rather than in cell order.
+    const auto approach = this->parent_.threePhaseApproach();
+    params_.materialLawParamArenas.resize(num_arrays);
+    for (unsigned i = 0; i < num_arrays; i++) {
+        params_.materialLawParamArenas[i] =
+            MaterialLawParams::makeArena(approach, this->numCompressedElems_);
+    }
     for (unsigned i = 0; i < num_arrays; i++) {
 #ifdef _OPENMP
 #pragma omp parallel for
@@ -111,7 +121,9 @@ run(const IntLookupFunction& fieldPropIntOnLeafAssigner,
                 hystParams.setImbibitionParamsGasWater(elemIdx, imbRegionIdx, lookupIdxOnLevelZeroAssigner);
             }
             hystParams.finalize();
-            initThreePhaseParams_(hystParams, (*mlpArray[i])[elemIdx], satRegionIdx, elemIdx);
+            initThreePhaseParams_(hystParams, (*mlpArray[i])[elemIdx], satRegionIdx, elemIdx,
+                                  MaterialLawParams::arenaSlot(params_.materialLawParamArenas[i],
+                                                               approach, elemIdx));
         }
     }
 }
@@ -229,14 +241,15 @@ InitParams<Traits>::
 initThreePhaseParams_(HystParams<Traits>& hystParams,
                       MaterialLawParams& materialParams,
                       unsigned satRegionIdx,
-                      unsigned elemIdx)
+                      unsigned elemIdx,
+                      std::shared_ptr<void> pooledParams)
 {
     const auto& epsInfo = this->params_.oilWaterScaledEpsInfoDrainage[elemIdx];
 
     auto oilWaterParams = hystParams.getOilWaterParams();
     auto gasOilParams = hystParams.getGasOilParams();
     auto gasWaterParams = hystParams.getGasWaterParams();
-    materialParams.setApproach(this->parent_.threePhaseApproach());
+    materialParams.setApproach(this->parent_.threePhaseApproach(), std::move(pooledParams));
     switch (materialParams.approach()) {
         case EclMultiplexerApproach::Stone1: {
             auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::Stone1>();

opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.hpp

@@ -33,6 +33,7 @@
 
 #include <cstddef>
 #include <functional>
+#include <memory>
 #include <string>
 #include <vector>
 
@@ -99,7 +100,8 @@ class InitParams
     void initThreePhaseParams_(HystParams<Traits>& hystParams,
                                MaterialLawParams& materialParams,
                                unsigned satRegionIdx,
-                               unsigned elemIdx);
+                               unsigned elemIdx,
+                               std::shared_ptr<void> pooledParams);
 
     void readEffectiveParameters_();
 

opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp

@@ -122,6 +122,10 @@ class Manager
         std::vector<int> imbnumRegionArray{};
         std::vector<MaterialLawParams> materialLawParams{};
         DirectionalMaterialLawParamsPtr dirMaterialLawParams{};
+        // Contiguous backing storage for the nested parameter objects the
+        // entries above point at: one arena per parameter array instead of
+        // one heap allocation per cell.
+        std::vector<std::shared_ptr<void>> materialLawParamArenas{};
         bool onlyPiecewiseLinear = true;
 
         bool hasDirectionalRelperms() const

opm/material/fluidmatrixinteractions/EclMultiplexerMaterialParams.hpp

@@ -33,8 +33,11 @@
 #include "EclTwoPhaseMaterial.hpp"
 
 #include <cassert>
+#include <cstddef>
 #include <memory>
 #include <type_traits>
+#include <utility>
+#include <vector>
 
 #include <opm/material/common/EnsureFinalized.hpp>
 
@@ -120,6 +123,67 @@ class EclMultiplexerMaterialParams : public Traits, public EnsureFinalized
         return *this;
     }
 
+    /*!
+     * \brief Contiguous storage for \a n parameter objects of one approach.
+     *
+     * Hand the slots out with arenaSlot(). Keep the returned pointer alive for
+     * as long as any of them is in use.
+     */
+    static ParamPointerType makeArena(EclMultiplexerApproach approach, std::size_t n)
+    {
+        switch (approach) {
+        case EclMultiplexerApproach::Stone1:
+            return std::make_shared<std::vector<Stone1Params>>(n);
+        case EclMultiplexerApproach::Stone2:
+            return std::make_shared<std::vector<Stone2Params>>(n);
+        case EclMultiplexerApproach::Default:
+            return std::make_shared<std::vector<DefaultParams>>(n);
+        case EclMultiplexerApproach::TwoPhase:
+            return std::make_shared<std::vector<TwoPhaseParams>>(n);
+        case EclMultiplexerApproach::OnePhase:
+            break;  // no parameters
+        }
+        return {};
+    }
+
+    /*!
+     * \brief Slot \a i of an arena, as an aliasing pointer sharing the arena's
+     *        single control block.
+     */
+    static ParamPointerType arenaSlot(const ParamPointerType& arena,
+                                      EclMultiplexerApproach approach,
+                                      std::size_t i)
+    {
+        const auto at = [&arena, i](auto* tag) {
+            using V = std::vector<std::remove_pointer_t<decltype(tag)>>;
+            return ParamPointerType{arena, &(*std::static_pointer_cast<V>(arena))[i]};
+        };
+        switch (approach) {
+        case EclMultiplexerApproach::Stone1:   return at(static_cast<Stone1Params*>(nullptr));
+        case EclMultiplexerApproach::Stone2:   return at(static_cast<Stone2Params*>(nullptr));
+        case EclMultiplexerApproach::Default:  return at(static_cast<DefaultParams*>(nullptr));
+        case EclMultiplexerApproach::TwoPhase: return at(static_cast<TwoPhaseParams*>(nullptr));
+        case EclMultiplexerApproach::OnePhase: break;
+        }
+        return {};
+    }
+
+    /*!
+     * \brief Point this object at storage owned by someone else.
+     *
+     * The caller keeps ownership; \a pooled is expected to be an aliasing
+     * shared_ptr into an arena, so a whole grid's worth of parameter objects
+     * costs one allocation and one control block instead of one of each per
+     * cell. Type and lifetime are the caller's responsibility, exactly as for
+     * the self-allocating overload.
+     */
+    void setApproach(EclMultiplexerApproach newApproach, ParamPointerType pooled)
+    {
+        assert(realParams_ == nullptr);
+        approach_ = newApproach;
+        realParams_ = std::move(pooled);
+    }
+
     void setApproach(EclMultiplexerApproach newApproach)
     {
         assert(realParams_ == 0);