doxygen

Author: rrsettgast

src/common/logMath.hpp

@@ -1,124 +1,283 @@
-
 #pragma once
 
 #include "macros.hpp"
 #include <cmath>
 
-
+/**
+ * @file logMath.hpp
+ * @brief Base-selectable logarithm and exponential helpers with optional fast CUDA intrinsics.
+ *
+ * @details
+ * This header provides a small, header-only utility for computing `log` and `exp` in either:
+ * - natural base (e), or
+ * - base-10 (ten),
+ * and optionally using CUDA "fast math" intrinsics on device code.
+ *
+ * Key points:
+ * - For `float` on CUDA device (`__CUDA_ARCH__`), `MathMode::fast` selects intrinsics such as
+ *   `__logf`, `__log10f`, `__expf`, `__exp10f` when available.
+ * - For host code (and for non-fast mode), it falls back to standard `<cmath>` functions
+ *   like `::log`, `::log10`, `::exp`, `::expf`, `::logf`, `::log10f`.
+ * - Some libcs (glibc) provide `exp10/exp10f` as extensions. When unavailable, base-10
+ *   exponentiation is implemented via `exp(x * ln(10))`.
+ *
+ * @note
+ * The functions here assume their usual mathematical domains:
+ * - `log(x)` requires `x > 0`.
+ * - `exp(x)` is defined for all real `x` but may overflow for large `x`.
+ *
+ * @warning
+ * When `MathMode::fast` is used on CUDA device, results may differ slightly from the
+ * accurate host/standard-library implementations due to reduced precision/approximations.
+ */
 
 namespace hpcReact
 {
 
+/**
+ * @def HPCREACT_HAS_EXP10
+ * @brief Indicates whether `exp10/exp10f` are available as libc extensions.
+ *
+ * @details
+ * glibc commonly provides `exp10` and `exp10f` as non-standard extensions. This macro is used to
+ * select a direct call to those functions on host when available; otherwise a mathematically
+ * equivalent fallback is used:
+ * \f[
+ *   10^x = e^{x \ln(10)}
+ * \f]
+ */
 #if defined(__GLIBC__)  // glibc provides exp10/exp10f as an extension
   #define HPCREACT_HAS_EXP10 1
 #else
   #define HPCREACT_HAS_EXP10 0
 #endif
 
-enum class LogBase { e, ten };
-enum class MathMode { accurate, fast };
+/**
+ * @brief Enumerates the logarithm base used by the helper.
+ */
+enum class LogBase
+{
+  e,   /**< Natural logarithm base (ln). */
+  ten  /**< Base-10 logarithm. */
+};
+
+/**
+ * @brief Enumerates the math evaluation mode.
+ *
+ * @details
+ * `fast` is only used for `float` on CUDA device code. On host (or for `double`), the
+ * standard library functions are used.
+ */
+enum class MathMode
+{
+  accurate, /**< Prefer standard-library implementations. */
+  fast      /**< Prefer CUDA fast intrinsics for `float` on device when available. */
+};
 
+/**
+ * @brief Base-selectable logarithm and exponential functions.
+ *
+ * @tparam Base Selects the logarithm/exp base:
+ *   - `LogBase::e` for natural base
+ *   - `LogBase::ten` for base-10
+ * @tparam Mode Selects evaluation mode. See @ref MathMode.
+ *
+ * @details
+ * Provides:
+ * - `log(x)` returning \f$\ln(x)\f$ or \f$\log_{10}(x)\f$ depending on `Base`
+ * - `exp(x)` returning \f$e^x\f$ or \f$10^x\f$ depending on `Base`
+ *
+ * The overload set covers `float` and `double`.
+ */
 template< LogBase Base, MathMode Mode = MathMode::accurate >
 struct LogExp
 {
-
-  template<class T>
-  HPCREACT_HOST_DEVICE static constexpr 
-  T ln10()
+  /**
+   * @brief Compile-time constant for \f$\ln(10)\f$.
+   *
+   * @tparam T Floating-point type (typically `float` or `double`).
+   * @return \f$\ln(10)\f$ as type `T`.
+   *
+   * @details
+   * The literal is provided with long-double precision and cast to `T`.
+   */
+  template< class T >
+  HPCREACT_HOST_DEVICE
+  static constexpr T ln10()
   {
     return T(2.3025850929940456840179914546843642L);
   }
 
+  /**
+   * @brief Conversion factor for derivatives with respect to `log(x)` vs `ln(x)`.
+   *
+   * @tparam T Floating-point type (typically `float` or `double`).
+   * @return A constant factor depending on @p Base.
+   *
+   * @details
+   * This returns the scalar \f$c\f$ such that:
+   * \f[
+   *   \frac{d}{d(\log_b x)} = c \, \frac{d}{d(\ln x)}
+   * \f]
+   * where:
+   * - if `Base == e`, then \f$c = 1\f$
+   * - if `Base == ten`, then \f$c = \ln(10)\f$
+   *
+   * Equivalently, since \f$\log_{10}(x) = \ln(x) / \ln(10)\f$:
+   * \f[
+   *   \frac{d}{d(\log_{10} x)} = \ln(10) \, \frac{d}{d(\ln x)}
+   * \f]
+   */
   template< typename T >
   HPCREACT_HOST_DEVICE
-  static constexpr inline 
-  T dWrtLogConst()
+  static constexpr inline T dWrtLogConst()
   {
     if constexpr ( Base == LogBase::e ) { return T(1.0); }
     else                                { return T(ln10<T>()); }
   }
 
+  // ***** log function *********************************************************************************************
 
-
-// ***** log function *********************************************************************************************
-  HPCREACT_HOST_DEVICE 
-  static inline 
-  double log( double const x )
+  /**
+   * @brief Logarithm in the selected base for `double`.
+   *
+   * @param x Input value (must be > 0).
+   * @return `ln(x)` if `Base == LogBase::e`, otherwise `log10(x)`.
+   */
+  HPCREACT_HOST_DEVICE
+  static constexpr inline double log( double const x )
   {
-    if constexpr ( Base == LogBase::e ) { return ::log(x); }
-    else                                { return ::log10(x); }
+    if constexpr ( Base == LogBase::e ) { return ::log( x ); }
+    else                                { return ::log10( x ); }
   }
-  
-  HPCREACT_HOST_DEVICE 
-  static inline 
-  float log( float const x )
+
+  /**
+   * @brief Logarithm in the selected base for `float`.
+   *
+   * @param x Input value (must be > 0).
+   * @return `ln(x)` if `Base == LogBase::e`, otherwise `log10(x)`.
+   *
+   * @details
+   * On CUDA device code (`__CUDA_ARCH__`) and when `Mode == MathMode::fast`,
+   * uses CUDA intrinsics:
+   * - `__logf` / `__log10f`
+   * Otherwise falls back to `<cmath>`:
+   * - `::logf` / `::log10f`
+   */
+  HPCREACT_HOST_DEVICE
+  static inline float log( float const x )
   {
 #if defined(__CUDA_ARCH__)
     if constexpr ( Mode == MathMode::fast )
     {
-      if constexpr ( Base == LogBase::e) { return __logf(x); }
-      else                               { return __log10f(x); }
+      if constexpr ( Base == LogBase::e ) { return __logf( x ); }
+      else                                { return __log10f( x ); }
     }
 #endif
-    if constexpr ( Base == LogBase::e ) { return ::logf(x); }
-    else                                { return ::log10f(x); }
+    if constexpr ( Base == LogBase::e ) { return ::logf( x ); }
+    else                                { return ::log10f( x ); }
   }
 
-// ***** exp function *********************************************************************************************
-  HPCREACT_HOST_DEVICE 
-  static inline 
-  double exp( double const x )
+  // ***** exp function *********************************************************************************************
+
+  /**
+   * @brief Exponential in the selected base for `double`.
+   *
+   * @param x Exponent.
+   * @return `exp(x)` if `Base == LogBase::e`, otherwise `10^x`.
+   *
+   * @details
+   * For base-10:
+   * - If `HPCREACT_HAS_EXP10` is true, uses `::exp10(x)` (glibc extension).
+   * - Otherwise uses `::exp(x * ln(10))`.
+   */
+  HPCREACT_HOST_DEVICE
+  static inline double exp( double const x )
   {
-    if constexpr ( Base == LogBase::e) 
-    { 
-      return ::exp(x); 
+    if constexpr ( Base == LogBase::e )
+    {
+      return ::exp( x );
     }
     else
-    { 
-      #if HPCREACT_HAS_EXP10
-        return ::exp10(x);
-      #else
-        return ::exp( x * ln10<double>() ); ; 
-      #endif
+    {
+#if HPCREACT_HAS_EXP10
+      return ::exp10( x );
+#else
+      return ::exp( x * ln10<double>() );
+#endif
     }
   }
 
-
-  HPCREACT_HOST_DEVICE 
-  static inline 
-  float exp( float const x )
+  /**
+   * @brief Exponential in the selected base for `float`.
+   *
+   * @param x Exponent.
+   * @return `exp(x)` if `Base == LogBase::e`, otherwise `10^x`.
+   *
+   * @details
+   * On CUDA device code (`__CUDA_ARCH__`) and when `Mode == MathMode::fast`,
+   * uses CUDA intrinsics:
+   * - `__expf` / `__exp10f`
+   *
+   * Otherwise:
+   * - For base-e uses `::expf(x)`
+   * - For base-10 uses `::exp10f(x)` if available (glibc extension), else `::expf(x * ln(10))`.
+   */
+  HPCREACT_HOST_DEVICE
+  static inline float exp( float const x )
   {
 #if defined(__CUDA_ARCH__)
     if constexpr ( Mode == MathMode::fast )
     {
-      if constexpr ( Base == LogBase::e ) { return __expf(x); }
-      else                                { return __exp10f(x); }
+      if constexpr ( Base == LogBase::e ) { return __expf( x ); }
+      else                                { return __exp10f( x ); }
     }
 #endif
-    if constexpr ( Base == LogBase::e) 
-    { 
-      return ::expf(x); 
+    if constexpr ( Base == LogBase::e )
+    {
+      return ::expf( x );
     }
     else
-    { 
+    {
 #if HPCREACT_HAS_EXP10
-      return ::exp10f(x);
+      return ::exp10f( x );
 #else
-      return ::expf( x * ln10<float>() ); 
+      return ::expf( x * ln10<float>() );
 #endif
     }
   }
 }; // struct LogExp
 
-
+/**
+ * @def HPC_REACT_LOG_TYPE
+ * @brief Compile-time selection of the default `logmath` alias.
+ *
+ * @details
+ * Values:
+ * - `0`: natural log/exp (base-e), fast mode
+ * - `1`: base-10 log/exp, fast mode (default)
+ *
+ * This macro is intended to allow build-time selection of the logarithm base used in code
+ * that relies on the `hpcReact::logmath` alias.
+ */
 #if !defined(HPC_REACT_LOG_TYPE)
 #define HPC_REACT_LOG_TYPE 1
 #endif
 
+/**
+ * @brief Default log/exp helper type selected by `HPC_REACT_LOG_TYPE`.
+ *
+ * @details
+ * - `HPC_REACT_LOG_TYPE == 0` selects `LogBase::e`
+ * - `HPC_REACT_LOG_TYPE == 1` selects `LogBase::ten`
+ *
+ * Both selections currently use `MathMode::fast`. Note that `MathMode::fast` only changes
+ * behavior for `float` on CUDA device code.
+ */
 #if HPC_REACT_LOG_TYPE == 0
   using logmath = LogExp< LogBase::e, MathMode::fast >;
 #elif HPC_REACT_LOG_TYPE == 1
   using logmath = LogExp< LogBase::ten, MathMode::fast >;
 #endif
 
-}  // namespace hpcReact
+}  // namespace hpcReact

src/reactions/exampleSystems/unitTests/testEquilibriumReactions.cpp

@@ -27,11 +27,11 @@ TEST( testEquilibriumReactions, computeResidualAndJacobianTest )
 
   {
     std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
-    double const expectedResiduals[] = { -37.534508668465 / logmath::dWrtLogConst<double>(), 
-                                         -72.989575795250 / logmath::dWrtLogConst<double>() };
+    double const expectedResiduals[] = { -37.534508668465 / logmath::dWrtLogConst< double >(),
+                                         -72.989575795250 / logmath::dWrtLogConst< double >() };
     double const expectedJacobian[2][2] =
-    { { 1.0e16 / logmath::dWrtLogConst<double>(), -2.0 / logmath::dWrtLogConst<double>() },
-      { -2.0 / logmath::dWrtLogConst<double>(), 4.0e16 / logmath::dWrtLogConst<double>() } };
+    { { 1.0e16 / logmath::dWrtLogConst< double >(), -2.0 / logmath::dWrtLogConst< double >() },
+      { -2.0 / logmath::dWrtLogConst< double >(), 4.0e16 / logmath::dWrtLogConst< double >() } };
 
     computeResidualAndJacobianTest< double, 2 >( bulkGeneric::simpleTestRateParams,
                                                  initialSpeciesConcentration,

src/reactions/exampleSystems/unitTests/testKineticReactions.cpp

@@ -59,7 +59,8 @@ TEST( testKineticReactions, computeSpeciesRatesTest_simpleKineticTestRateParams
 // TEST( testKineticReactions, testTimeStep )
 // {
 //   double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
-//   double const expectedSpeciesConcentrations[5] = { 3.92138293924124e-01, 3.03930853037938e-01, 5.05945480771998e-01, 7.02014627734060e-01, 5.95970744531880e-01 };
+//   double const expectedSpeciesConcentrations[5] = { 3.92138293924124e-01, 3.03930853037938e-01, 5.05945480771998e-01,
+// 7.02014627734060e-01, 5.95970744531880e-01 };
 
 //   timeStepTest< double >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
 //                           2.0,

src/reactions/massActions/unitTests/testMassActions.cpp

@@ -188,7 +188,7 @@ void testcalculateAggregatePrimaryConcentrationsWrtLogCHelper()
     for( int j=0; j<numPrimarySpecies; ++j )
     {
       EXPECT_NEAR( data.dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, j ),
-                   logmath::dWrtLogConst<double>() * expected_dAggregatePrimarySpeciesConcentration_dLogPrimarySpeciesConcentrations[i][j],
+                   logmath::dWrtLogConst< double >() * expected_dAggregatePrimarySpeciesConcentration_dLogPrimarySpeciesConcentrations[i][j],
                    1.0e-8 );
     }
   }

src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp

@@ -103,7 +103,7 @@ EquilibriumReactions< REAL_TYPE,
     // compute the jacobian
     for( IndexType b=0; b<numReactions; ++b )
     {
-      jacobian( a, b ) = ( -dForwardProduct_dxi_divProduct[b] + dReverseProduct_dxi_divProduct[b] ) / logmath::dWrtLogConst<REAL_TYPE>();
+      jacobian( a, b ) = ( -dForwardProduct_dxi_divProduct[b] + dReverseProduct_dxi_divProduct[b] ) / logmath::dWrtLogConst< REAL_TYPE >();
     }
   }
 }

src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp

@@ -255,7 +255,7 @@ void timeStepTest( PARAMS_DATA const & params,
           {
             printf( "%e ", dataCopy->speciesConcentration[i] );
           }
-          printf( "\n" );                                          
+          printf( "\n" );
           dataCopy->time += dt;
         }
       } );