From 1cb73b756afee623e0f6e8e963c580adfc46d1a8 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 12 Aug 2026 11:35:17 +0200 Subject: [PATCH] Allocate the per-cell material law parameters from an arena EclMultiplexerMaterialParams::setApproach() allocates the nested parameter object with a bare new, so a grid costs one allocation and one shared_ptr control block per cell. On Norne that is 44927 allocations of 304 B spread over 105.7 MB of address space for 13.7 MB of data. Add makeArena()/arenaSlot() and a setApproach() overload taking externally owned storage, and have the manager hand out aliasing pointers into one vector per parameter array. The aliasing constructor shares the arena's single control block, so the per-cell overhead drops to the pointer itself and the objects land at a fixed 304 B stride. This is not a measured speedup: on Norne the relperm+pc evaluation loop is unchanged in cell order (91.5 -> 91.3 ns/cell, best of 5 alternating runs) and at best 2-3% in shuffled order, which is inside the noise. The allocator already hands out long contiguous runs when the objects are built back to back. What changes is the allocation count and the address span. All 228 opm-common tests pass. --- .../EclMaterialLawInitParams.cpp | 19 +++++- .../EclMaterialLawInitParams.hpp | 4 +- .../EclMaterialLawManager.hpp | 4 ++ .../EclMultiplexerMaterialParams.hpp | 64 +++++++++++++++++++ 4 files changed, 87 insertions(+), 4 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.cpp index ec52934e763..f575e4214eb 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.cpp @@ -85,6 +85,16 @@ run(const IntLookupFunction& fieldPropIntOnLeafAssigner, std::vector*> 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:: initThreePhaseParams_(HystParams& hystParams, MaterialLawParams& materialParams, unsigned satRegionIdx, - unsigned elemIdx) + unsigned elemIdx, + std::shared_ptr 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(); diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.hpp index 7b286dc9e0b..cf037ffae5d 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.hpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.hpp @@ -33,6 +33,7 @@ #include #include +#include #include #include @@ -99,7 +100,8 @@ class InitParams void initThreePhaseParams_(HystParams& hystParams, MaterialLawParams& materialParams, unsigned satRegionIdx, - unsigned elemIdx); + unsigned elemIdx, + std::shared_ptr pooledParams); void readEffectiveParameters_(); diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp index 920df2bd7d4..8302dc231d8 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp @@ -122,6 +122,10 @@ class Manager std::vector imbnumRegionArray{}; std::vector 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> materialLawParamArenas{}; bool onlyPiecewiseLinear = true; bool hasDirectionalRelperms() const diff --git a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterialParams.hpp b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterialParams.hpp index a2caeefcdad..a5bf6044a65 100644 --- a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterialParams.hpp +++ b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterialParams.hpp @@ -33,8 +33,11 @@ #include "EclTwoPhaseMaterial.hpp" #include +#include #include #include +#include +#include #include @@ -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>(n); + case EclMultiplexerApproach::Stone2: + return std::make_shared>(n); + case EclMultiplexerApproach::Default: + return std::make_shared>(n); + case EclMultiplexerApproach::TwoPhase: + return std::make_shared>(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>; + return ParamPointerType{arena, &(*std::static_pointer_cast(arena))[i]}; + }; + switch (approach) { + case EclMultiplexerApproach::Stone1: return at(static_cast(nullptr)); + case EclMultiplexerApproach::Stone2: return at(static_cast(nullptr)); + case EclMultiplexerApproach::Default: return at(static_cast(nullptr)); + case EclMultiplexerApproach::TwoPhase: return at(static_cast(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);