cleanup use of logarithms

src/common/logMath.hpp

@@ -0,0 +1,124 @@
+
+#pragma once
+
+#include "macros.hpp"
+#include <cmath>
+
+
+
+namespace hpcReact
+{
+
+#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 };
+
+template< LogBase Base, MathMode Mode = MathMode::accurate >
+struct LogExp
+{
+
+  template<class T>
+  HPCREACT_HOST_DEVICE static constexpr 
+  T ln10()
+  {
+    return T(2.3025850929940456840179914546843642L);
+  }
+
+  template< typename T >
+  HPCREACT_HOST_DEVICE
+  static constexpr inline 
+  T dWrtLogConst()
+  {
+    if constexpr ( Base == LogBase::e ) { return T(1.0); }
+    else                                { return T(ln10<T>()); }
+  }
+
+
+
+// ***** log function *********************************************************************************************
+  HPCREACT_HOST_DEVICE 
+  static inline 
+  double log( double const x )
+  {
+    if constexpr ( Base == LogBase::e ) { return ::log(x); }
+    else                                { return ::log10(x); }
+  }
+
+  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); }
+    }
+#endif
+    if constexpr ( Base == LogBase::e ) { return ::logf(x); }
+    else                                { return ::log10f(x); }
+  }
+
+// ***** exp function *********************************************************************************************
+  HPCREACT_HOST_DEVICE 
+  static inline 
+  double exp( double const x )
+  {
+    if constexpr ( Base == LogBase::e) 
+    { 
+      return ::exp(x); 
+    }
+    else
+    { 
+      #if HPCREACT_HAS_EXP10
+        return ::exp10(x);
+      #else
+        return ::exp( x * ln10<double>() ); ; 
+      #endif
+    }
+  }
+
+
+  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); }
+    }
+#endif
+    if constexpr ( Base == LogBase::e) 
+    { 
+      return ::expf(x); 
+    }
+    else
+    { 
+#if HPCREACT_HAS_EXP10
+      return ::exp10f(x);
+#else
+      return ::expf( x * ln10<float>() ); 
+#endif
+    }
+  }
+}; // struct LogExp
+
+
+#if !defined(HPC_REACT_LOG_TYPE)
+#define HPC_REACT_LOG_TYPE 1
+#endif
+
+#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

src/common/macros.hpp

@@ -39,42 +39,42 @@
 
 #if defined( __clang__ )
 #define HPCREACT_NO_MISSING_BRACES( ... ) \
-        _Pragma("clang diagnostic push") \
-        _Pragma("clang diagnostic ignored \"-Wmissing-braces\"") \
-        __VA_ARGS__ \
-        _Pragma("clang diagnostic pop")
+  _Pragma("clang diagnostic push") \
+  _Pragma("clang diagnostic ignored \"-Wmissing-braces\"") \
+  __VA_ARGS__ \
+  _Pragma("clang diagnostic pop")
 
 #define HPCREACT_NO_MISSING_BRACES_OPEN \
-        _Pragma("clang diagnostic push") \
-        _Pragma("clang diagnostic ignored \"-Wmissing-braces\"")
+  _Pragma("clang diagnostic push") \
+  _Pragma("clang diagnostic ignored \"-Wmissing-braces\"")
 #define HPCREACT_NO_MISSING_BRACES_CLOSE \
-        _Pragma("clang diagnostic pop")
+  _Pragma("clang diagnostic pop")
 
 #elif defined(__GNUC__)
 #define HPCREACT_NO_MISSING_BRACES( ... ) \
-        _Pragma("GCC diagnostic push") \
-        _Pragma("GCC diagnostic ignored \"-Wmissing-braces\"") \
-        __VA_ARGS__ \
-        _Pragma("GCC diagnostic pop")
+  _Pragma("GCC diagnostic push") \
+  _Pragma("GCC diagnostic ignored \"-Wmissing-braces\"") \
+  __VA_ARGS__ \
+  _Pragma("GCC diagnostic pop")
 
 #define HPCREACT_NO_MISSING_BRACES_OPEN \
-        _Pragma("GCC diagnostic push") \
-        _Pragma("GCC diagnostic ignored \"-Wmissing-braces\"")
+  _Pragma("GCC diagnostic push") \
+  _Pragma("GCC diagnostic ignored \"-Wmissing-braces\"")
 #define HPCREACT_NO_MISSING_BRACES_CLOSE \
-        _Pragma("GCC diagnostic pop")
+  _Pragma("GCC diagnostic pop")
 
 #elif defined(_MSC_VER)
 #define HPCREACT_NO_MISSING_BRACES( ... ) \
-        __pragma(warning(push)) \
-        __pragma(warning(disable : 4351)) \
-        __VA_ARGS__ \
-        __pragma(warning(pop))
+  __pragma(warning(push)) \
+  __pragma(warning(disable : 4351)) \
+  __VA_ARGS__ \
+  __pragma(warning(pop))
 
 #define HPCREACT_NO_MISSING_BRACES_OPEN \
-        __pragma(warning(push)) \
-        __pragma(warning(disable : 4351))
+  __pragma(warning(push)) \
+  __pragma(warning(disable : 4351))
 #define HPCREACT_NO_MISSING_BRACES_CLOSE \
-        __pragma(warning(pop))
+  __pragma(warning(pop))
 #else
 #define HPCREACT_NO_MISSING_BRACES( ... ) __VA_ARGS__ // No-op for unknown compilers
 #define HPCREACT_NO_MISSING_BRACES_OPEN

src/common/pmpl.hpp

@@ -22,17 +22,17 @@ namespace hpcReact
 {
 #if defined(HPCREACT_USE_DEVICE)
   #if defined(HPCREACT_USE_CUDA)
-    #define deviceMalloc( PTR, BYTES ) cudaMalloc( PTR, BYTES );
-    #define deviceMallocManaged( PTR, BYTES ) cudaMallocManaged( PTR, BYTES );
-    #define deviceDeviceSynchronize() cudaDeviceSynchronize();
-    #define deviceMemCpy( DST, SRC, BYTES, KIND ) cudaMemcpy( DST, SRC, BYTES, KIND );
-    #define deviceFree( PTR ) cudaFree( PTR );
+#define deviceMalloc( PTR, BYTES ) cudaMalloc( PTR, BYTES );
+#define deviceMallocManaged( PTR, BYTES ) cudaMallocManaged( PTR, BYTES );
+#define deviceDeviceSynchronize() cudaDeviceSynchronize();
+#define deviceMemCpy( DST, SRC, BYTES, KIND ) cudaMemcpy( DST, SRC, BYTES, KIND );
+#define deviceFree( PTR ) cudaFree( PTR );
   #elif defined(HPCREACT_USE_HIP)
-    #define deviceMalloc( PTR, BYTES ) hipMalloc( PTR, BYTES );
-    #define deviceMallocManaged( PTR, BYTES ) hipMallocManaged( PTR, BYTES );
-    #define deviceDeviceSynchronize() hipDeviceSynchronize();
-    #define deviceMemCpy( DST, SRC, BYTES, KIND ) hipMemcpy( DST, SRC, BYTES, KIND );
-    #define deviceFree( PTR ) hipFree( PTR );
+#define deviceMalloc( PTR, BYTES ) hipMalloc( PTR, BYTES );
+#define deviceMallocManaged( PTR, BYTES ) hipMallocManaged( PTR, BYTES );
+#define deviceDeviceSynchronize() hipDeviceSynchronize();
+#define deviceMemCpy( DST, SRC, BYTES, KIND ) hipMemcpy( DST, SRC, BYTES, KIND );
+#define deviceFree( PTR ) hipFree( PTR );
   #endif
 #endif
 
@@ -130,12 +130,18 @@ void genericKernelWrapper( int const N, DATA_TYPE * const hostData, LAMBDA && fu
   genericKernel <<< 1, 1 >>> ( std::forward< LAMBDA >( func ), deviceData );
 
   cudaError_t e = cudaGetLastError();
-  if (e != cudaSuccess) { fprintf(stderr, "launch error: %s\n", cudaGetErrorString(e)); abort(); }
+  if( e != cudaSuccess )
+  {
+    fprintf( stderr, "launch error: %s\n", cudaGetErrorString( e )); abort();
+  }
 
   deviceDeviceSynchronize();
 
   e = cudaGetLastError();
-  if (e != cudaSuccess) { fprintf(stderr, "post-sync error: %s\n", cudaGetErrorString(e)); abort(); }
+  if( e != cudaSuccess )
+  {
+    fprintf( stderr, "post-sync error: %s\n", cudaGetErrorString( e )); abort();
+  }
 
   deviceMemCpy( hostData, deviceData, N * sizeof(DATA_TYPE), cudaMemcpyDeviceToHost );
   deviceFree( deviceData );

src/reactions/exampleSystems/unitTests/testEquilibriumReactions.cpp

@@ -27,10 +27,11 @@ TEST( testEquilibriumReactions, computeResidualAndJacobianTest )
 
   {
     std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
-    double const expectedResiduals[] = { -37.534508668465, -72.989575795250 };
+    double const expectedResiduals[] = { -37.534508668465 / logmath::dWrtLogConst<double>(), 
+                                         -72.989575795250 / logmath::dWrtLogConst<double>() };
     double const expectedJacobian[2][2] =
-    { { 1.0e16, -2.0 },
-      { -2.0, 4.0e16 } };
+    { { 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/testGenericChainReactions.cpp

@@ -44,16 +44,11 @@ TEST( testChainGenericKineticReactions, computeReactionRatesTest_chainReactionPa
   { { 0.05, 0.0, 0.0 },
     { 0.0, 0.03, 0.0 },
     { 0.0, 0.0, 0.02 } };
-  computeReactionRatesTest< double, false >( serialAllKineticParams.kineticReactionsParameters(),
-                                             initialSpeciesConcentration,
-                                             surfaceArea, // No use. Just to pass something here
-                                             expectedReactionRates,
-                                             expectedReactionRatesDerivatives );
-  computeReactionRatesTest< double, true >( serialAllKineticParams.kineticReactionsParameters(),
-                                            initialSpeciesConcentration,
-                                            surfaceArea, // No use. Just to pass something here
-                                            expectedReactionRates,
-                                            expectedReactionRatesDerivatives );
+  computeReactionRatesTest< double >( serialAllKineticParams.kineticReactionsParameters(),
+                                      initialSpeciesConcentration,
+                                      surfaceArea,       // No use. Just to pass something here
+                                      expectedReactionRates,
+                                      expectedReactionRatesDerivatives );
 }
 
 int main( int argc, char * * argv )

src/reactions/exampleSystems/unitTests/testKineticReactions.cpp

@@ -29,16 +29,11 @@ TEST( testKineticReactions, computeReactionRatesTest_simpleKineticTestRateParams
   double const expectedReactionRatesDerivatives[][5] =
   { { 2.0, -0.5, 0.0, 0.0, 0.0 },
     { 0.0, 0.0, 0.5, 0.25, 0.0 } };
-  computeReactionRatesTest< double, false >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
-                                             initialSpeciesConcentration,
-                                             surfaceArea, // No use. Just to pass something here
-                                             expectedReactionRates,
-                                             expectedReactionRatesDerivatives );
-  computeReactionRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
-                                            initialSpeciesConcentration,
-                                            surfaceArea, // No use. Just to pass something here
-                                            expectedReactionRates,
-                                            expectedReactionRatesDerivatives );
+  computeReactionRatesTest< double >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
+                                      initialSpeciesConcentration,
+                                      surfaceArea,       // No use. Just to pass something here
+                                      expectedReactionRates,
+                                      expectedReactionRatesDerivatives );
 }
 
 
@@ -52,37 +47,27 @@ TEST( testKineticReactions, computeSpeciesRatesTest_simpleKineticTestRateParams
     {  0.0, 0.0, -0.5, -0.25, 0.0 },
     {  0.0, 0.0, 1.0, 0.5, 0.0 } };
 
-  computeSpeciesRatesTest< double, false >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
-                                            initialSpeciesConcentration,
-                                            expectedSpeciesRates,
-                                            expectedSpeciesRatesDerivatives );
 
-  computeSpeciesRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
-                                           initialSpeciesConcentration,
-                                           expectedSpeciesRates,
-                                           expectedSpeciesRatesDerivatives );
+  computeSpeciesRatesTest< double >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
+                                     initialSpeciesConcentration,
+                                     expectedSpeciesRates,
+                                     expectedSpeciesRatesDerivatives );
 
 }
 
 
-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 };
-
-  timeStepTest< double, false >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
-                                 2.0,
-                                 10,
-                                 initialSpeciesConcentration,
-                                 expectedSpeciesConcentrations );
-
-  // ln(c) as the primary variable results in a singular system.
-  // timeStepTest< double, true >( simpleKineticTestRateParams,
-  //                               2.0,
-  //                               10,
-  //                               initialSpeciesConcentration,
-  //                               expectedSpeciesConcentrations );
-}
+// 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 };
+
+//   timeStepTest< double >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
+//                           2.0,
+//                           10,
+//                           initialSpeciesConcentration,
+//                           expectedSpeciesConcentrations );
+
+// }
 
 int main( int argc, char * * argv )
 {

src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp

@@ -53,19 +53,19 @@ void testMoMasAllEquilibriumHelper()
 
     double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] =
     {
-      log( initialPrimarySpeciesConcentration[0] ),
-      log( initialPrimarySpeciesConcentration[1] ),
-      log( initialPrimarySpeciesConcentration[2] ),
-      log( initialPrimarySpeciesConcentration[3] ),
-      log( initialPrimarySpeciesConcentration[4] )
+      logmath::log( initialPrimarySpeciesConcentration[0] ),
+      logmath::log( initialPrimarySpeciesConcentration[1] ),
+      logmath::log( initialPrimarySpeciesConcentration[2] ),
+      logmath::log( initialPrimarySpeciesConcentration[3] ),
+      logmath::log( initialPrimarySpeciesConcentration[4] )
     };
 
     EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
                                                             hpcReact::MoMasBenchmark::easyCaseParams.equilibriumReactionsParameters(),
                                                             targetAggregatePrimarySpeciesConcentration,
                                                             logInitialPrimarySpeciesConcentration,
                                                             logPrimarySpeciesConcentrationCopy );
-  });
+  } );
 
   double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] =
   {
@@ -78,7 +78,7 @@ void testMoMasAllEquilibriumHelper()
 
   for( int r=0; r<numPrimarySpecies; ++r )
   {
-    EXPECT_NEAR( exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
+    EXPECT_NEAR( logmath::exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
   }
 
 }