Merge pull request #470 from lanl/jmm/indexable-types

Jmm/indexable types

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR470]](https://github.com/lanl/singularity-eos/pull/470) Add the ability to access lambda elements by named types
 - [[PR459]](https://github.com/lanl/singularity-eos/pull/459) Add electron and ion tables to EOSPAC and SpinerEOS backends
 - [[PR453]](https://github.com/lanl/singularity-eos/pull/453) A PT space PTE solver
 - [[PR444]](https://github.com/lanl/singularity-eos/pull/444) Add Z split modifier and electron ideal gas EOS

doc/sphinx/src/using-closures.rst

@@ -457,8 +457,9 @@ and
 
 provide the number of real numbers (i.e., either ``float`` or
 ``double``) required for a single cell given a number of materials in
-equilibriun for either the ``RhoT`` or ``RhoU`` solver. The equivalent
-functions
+equilibriun for either the ``RhoT`` or ``RhoU`` solver.
+
+The equivalent functions
 
 .. cpp:function:: size_t PTESolverRhoTRequiredScratchInBytes(const int nmat);
 

doc/sphinx/src/using-eos.rst

@@ -659,6 +659,107 @@ minimize latency due to ``malloc`` and ``free``. Several models, such
 as ``SpinerEOS`` also use the persistency of these arrays to cache
 useful quantities for a performance boost.
 
+Named Arguments to Lambda Indexers
+-----------------------------------
+
+Lambdas support a more free-form kind of indexer. In particular, you
+may define indexers that take **names** via a type system. For
+example, you can write code such as:
+
+.. code-block:: cpp
+
+  Lambda_t lambda;
+  lambda[MeanIonizationState()] = zbar;
+  Real P = eos.PressureFromDensityTemperature(rho, T, lambda);
+
+The available types currently supported by default are:
+
+.. code-block:: cpp
+
+  struct singularity::IndexableTypes::MeanIonizationState;
+  struct singularity::IndexableTypes::LogDensity;
+  struct singularity::IndexableTypes::LogTemperature;
+  struct singularity::IndexableTypes::MeanAtomicMass;
+  struct singularity::IndexableTypes::MeanAtomicNumber;
+  struct singularity::IndexableTypes::ElectronFraction;
+
+However if you are not limited to these types. Any type will do and
+you can define your own as you like. For example:
+
+.. code-block::
+
+  namespace IndexableTypes {
+  struct MyLambdaParameter;
+  }
+
+To use an indexable type, you must define an indexer with an overload
+of the ``[]`` operator that takes your type. For example:
+
+.. code-block:: cpp
+
+  class MyLambda_t {
+  public:
+    MyLambda_t() = default;
+    PORTABLE_FORCEINLINE_FUNCTION
+    Real &operator[](const std::size_t idx) const {
+      return data_[idx];
+    }
+    PORTABLE_FORCEINLINE_FUNCTION
+    Real &operator[](const MeanIonizationState &zbar) const {
+      return data_[2];
+    }
+  private:
+    std::array<Real, 3> data_;
+  };
+
+which might be used as
+
+.. code-blocK:: cpp
+
+  MyLambda_t lambda;
+  lambda[0] = lRho;
+  lambda[1] = lT;
+  lambda[MeanIonizationState()] = Z;
+
+where ``MeanIonizationState`` is shorthand for index 2, since you
+defined that overload. Note that the ``operator[]`` must be marked
+``const``. To more easily enable mixing and matching integer-based
+indexing with type-based indexing, the function
+
+.. code-block:: cpp
+
+  template<typename Name_t, typename Indexer_t>
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real &Get(Indexer_t &&lambda, std::size_t idx = 0);
+
+will return a reference to the value at named index ``Name_t()`` if
+that overload is defined in ``Indexer_t`` and otherwise return a
+reference at index ``idx``.
+
+As a convenience tool, the struct
+
+.. code-block:: cpp
+
+  template<typename... Ts>
+  singularity::IndexableTypes::VariadicIndexer;
+
+automatically defines an indexer that accepts all named indices in the
+variadic list ``Ts...`` and also integer indexing. It's a fixed-size
+array under the hood.
+
+You can check if an equation of state is compatible with a given named
+index type by calling
+
+.. code-block:: cpp
+
+   eos.NeedsLambda(NamedIndex())
+
+which returns ``true`` if the EOS is compatible with ``NamedIndex``
+and ``false`` otherwise.
+
+All of the above functionality is available in
+the header file ``singularity-eos/base/indexable_types.hpp``.
+
 EOS Modifiers
 --------------
 

singularity-eos/CMakeLists.txt

@@ -25,6 +25,7 @@ register_headers(
 
     # Normal files
     base/fast-math/logs.hpp
+    base/indexable_types.hpp
     base/robust_utils.hpp
     base/root-finding-1d/root_finding.hpp
     base/serialization_utils.hpp

singularity-eos/base/indexable_types.hpp

@@ -0,0 +1,84 @@
+//------------------------------------------------------------------------------
+// © 2025. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#ifndef SINGULARITY_EOS_BASE_INDEXABLE_TYPES_
+#define SINGULARITY_EOS_BASE_INDEXABLE_TYPES_
+
+#include <array>
+#include <type_traits>
+#include <utility>
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-eos/base/variadic_utils.hpp>
+
+namespace singularity {
+namespace IndexerUtils {
+// Convenience function for accessing an indexer by either type or
+// natural number index depending on what is available
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION auto &Get(Indexer_t &&lambda, std::size_t idx = 0) {
+  if constexpr (variadic_utils::is_indexable_v<Indexer_t, T>) {
+    return lambda[T()];
+  } else {
+    return lambda[idx];
+  }
+}
+
+// This is a convenience struct to easily build a small indexer with
+// a set of indexable types.
+template <typename Data_t, typename... Ts>
+class VariadicIndexerBase {
+ public:
+  VariadicIndexerBase() = default;
+  PORTABLE_FORCEINLINE_FUNCTION
+  VariadicIndexerBase(const Data_t &data) : data_(data) {}
+  template <typename T,
+            typename = std::enable_if_t<variadic_utils::contains<T, Ts...>::value>>
+  PORTABLE_FORCEINLINE_FUNCTION Real &operator[](const T &t) {
+    constexpr std::size_t idx = variadic_utils::GetIndexInTL<T, Ts...>();
+    return data_[idx];
+  }
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real &operator[](const std::size_t idx) { return data_[idx]; }
+  template <typename T,
+            typename = std::enable_if_t<variadic_utils::contains<T, Ts...>::value>>
+  PORTABLE_FORCEINLINE_FUNCTION const Real &operator[](const T &t) const {
+    constexpr std::size_t idx = variadic_utils::GetIndexInTL<T, Ts...>();
+    return data_[idx];
+  }
+  PORTABLE_FORCEINLINE_FUNCTION
+  const Real &operator[](const std::size_t idx) const { return data_[idx]; }
+  static inline constexpr std::size_t size() { return sizeof...(Ts); }
+
+ private:
+  Data_t data_;
+};
+// uses an array
+template <typename... Ts>
+using VariadicIndexer = VariadicIndexerBase<std::array<Real, sizeof...(Ts)>, Ts...>;
+// uses a Real*
+template <typename... Ts>
+using VariadicPointerIndexer = VariadicIndexerBase<Real *, Ts...>;
+} // namespace IndexerUtils
+
+namespace IndexableTypes {
+struct MeanIonizationState {};
+struct LogDensity {};
+struct LogTemperature {};
+struct MeanAtomicMass {};
+struct MeanAtomicNumber {};
+struct ElectronFraction {};
+} // namespace IndexableTypes
+} // namespace singularity
+#endif // SINGULARITY_EOS_BASE_INDEXABLE_TYPES_

singularity-eos/base/variadic_utils.hpp

@@ -86,6 +86,31 @@ struct type_list {};
 template <template <typename> class... Ts>
 struct adapt_list {};
 
+// provide index of a type in a type list
+template <typename T, typename Head, typename... Ts>
+constexpr std::size_t GetIndexInTL(std::size_t current_index) {
+  if constexpr (std::is_same_v<T, Head>) {
+    return current_index;
+  } else {
+    static_assert(sizeof...(Ts) > 0, "Type T must be in type list!");
+    return GetIndexInTL<T, Ts...>(current_index + 1);
+  }
+}
+template <typename T, typename Head, typename... Ts>
+constexpr std::size_t GetIndexInTL() {
+  return GetIndexInTL<T, Head, Ts...>(0);
+}
+
+// is_indexable similar to is_invokable
+template <typename, typename, typename = void>
+struct is_indexable : std::false_type {};
+template <typename T, typename Index>
+struct is_indexable<T, Index,
+                    std::void_t<decltype(std::declval<T>()[std::declval<Index>()])>>
+    : std::true_type {};
+template <typename T, typename Index>
+constexpr bool is_indexable_v = is_indexable<T, Index>::value;
+
 // this flattens a typelist of typelists to a single typelist
 
 // first parameter - accumulator