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);