franks fixes

Author: rrsettgast

src/common/logMath.hpp

@@ -151,9 +151,10 @@ struct LogExp
    */
   template< typename T >
   HPCREACT_HOST_DEVICE
-  static constexpr inline void dWrtLogScale( T & x )
+  static constexpr inline T dWrtLogScale( T x )
   {
-    if constexpr ( Base == LogBase::ten ) { x = x * ln10<T>(); }
+    if constexpr ( Base == LogBase::ten ) { return x * ln10<T>(); }
+    else                                  { return x; }
   }
 
   // ***** log function *********************************************************************************************

src/reactions/exampleSystems/unitTests/testKineticReactions.cpp

@@ -41,7 +41,8 @@ TEST( testKineticReactions, computeSpeciesRatesTest_simpleKineticTestRateParams
 {
   double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
   double const expectedSpeciesRates[5] = { -2.0, 1.0, 0.75, -0.25, 0.5 };
-  double const expectedSpeciesRatesDerivatives[5][5] = { { -4.0, 1.0, 0.0, 0.0, 0.0 },
+  double const expectedSpeciesRatesDerivatives[5][5] =
+  { { -4.0, 1.0, 0.0, 0.0, 0.0 },
     {  2.0, -0.5, 0.0, 0.0, 0.0 },
     {  2.0, -0.5, -0.5, -0.25, 0.0 },
     {  0.0, 0.0, -0.5, -0.25, 0.0 },

src/reactions/reactionsSystems/KineticReactions_impl.hpp

@@ -62,7 +62,7 @@ KineticReactions< REAL_TYPE,
     // set reaction rate to zero
     reactionRates[r] = 0.0;
     // get/calculate the forward and reverse rate constants for this reaction
-    RealType const forwardRateConstant = params.rateConstantForward( r ); //* logmath::exp( -params.m_activationEnergy[r] / ( constants::R *
+    RealType const forwardRateConstant = params.rateConstantForward( r ); // exp( -params.m_activationEnergy[r] / ( constants::R *
                                                                           // temperature ) );
     RealType const reverseRateConstant = params.rateConstantReverse( r );
 
@@ -90,16 +90,18 @@ KineticReactions< REAL_TYPE,
 
     if constexpr( CALCULATE_DERIVATIVES )
     {
+      RealType const dFactorForward = logmath::dWrtLogScale( forwardRateConstant * logmath::exp( productConcForward ) );
+      RealType const dFactorReverse = logmath::dWrtLogScale( -reverseRateConstant * logmath::exp( productConcReverse ) );
       for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
       {
         RealType const s_ri = params.stoichiometricMatrix( r, i );
         if( s_ri < 0.0 )
         {
-          reactionRatesDerivatives( r, i ) = forwardRateConstant * logmath::exp( productConcForward ) * (-s_ri);
+          reactionRatesDerivatives( r, i ) = dFactorForward * (-s_ri);
         }
         else if( s_ri > 0.0 )
         {
-          reactionRatesDerivatives( r, i ) = -reverseRateConstant * logmath::exp( productConcReverse ) * s_ri;
+          reactionRatesDerivatives( r, i ) = dFactorReverse * s_ri;
         }
         else
         {
@@ -150,8 +152,8 @@ KineticReactions< REAL_TYPE,
     }
 
     // get/calculate the forward and reverse rate constants for this reaction
-    RealType const rateConstant = params.rateConstantForward( r ); //* logmath::exp( -params.m_activationEnergy[r] / ( constants::R *
-    // temperature ) );
+    RealType const rateConstant = params.rateConstantForward( r ); //* exp( -params.m_activationEnergy[r] / ( constants::R * temperature )
+                                                                   // );
     RealType const equilibriumConstant = params.equilibriumConstant( r );
 
 
@@ -257,7 +259,7 @@ KineticReactions< REAL_TYPE,
 
 
   REAL_TYPE residualNorm = 0.0;
-  for( int k=0; k<20; ++k ) // newton loop
+  for( int k=0; k<20; ++k )  // newton loop
   {
     printf( "iteration %2d: \n", k );
 

src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp

@@ -56,7 +56,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE,
                                                              ARRAY_1D_PRIMARY & aggregateSpeciesRates,
                                                              ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations )
 {
-  if constexpr ( PARAMS_DATA::numEquilibriumReactions() > 0 )
+  if constexpr( PARAMS_DATA::numEquilibriumReactions() > 0 )
   {
     // 1. Compute new aggregate species from primary species
     massActions::calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE,
@@ -168,7 +168,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE,
   {
     for( IntType j = 0; j < numPrimarySpecies; ++j )
     {
-      dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) = logmath::dWrtLogConst< REAL_TYPE >() * reactionRatesDerivatives( i, j + numSecondarySpecies );
+      dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) = reactionRatesDerivatives( i, j + numSecondarySpecies );
 
       for( IntType k = 0; k < numSecondarySpecies; ++k )
       {

src/reactions/reactionsSystems/Parameters.hpp

@@ -231,7 +231,7 @@ struct MixedReactionsParameters
         RealType const absDiff = fabs( K - ( kf / kr ) );
         RealType const effectiveMagnitude = max( fabs( K ), fabs( kf/kr ));
         RealType const tolerance = effectiveMagnitude * pow( 10, -num_digits );
-        if( absDiff > tolerance ) // Tolerance for floating point precision
+        if( absDiff > tolerance )  // Tolerance for floating point precision
         {
           throw std::runtime_error( "Error: Inconsistent equilibrium relation for reaction " + std::to_string( i ));
         }

src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp

@@ -113,7 +113,7 @@ void computeReactionRatesTest( PARAMS_DATA const & params,
     for( int i = 0; i < numSpecies; ++i )
     {
 
-      data.reactionRatesDerivatives( r, i ) = data.reactionRatesDerivatives( r, i ) * logmath::exp( -data.speciesConcentration[i] );
+      data.reactionRatesDerivatives( r, i ) = data.reactionRatesDerivatives( r, i ) * logmath::exp( -data.speciesConcentration[i] ) / logmath::dWrtLogConst< double >();
       EXPECT_NEAR( data.reactionRatesDerivatives( r, i ), expectedReactionRatesDerivatives[r][i], std::max( magScale, fabs( expectedReactionRatesDerivatives[r][i] ) ) * 1.0e-8 );
     }
   }
@@ -180,7 +180,7 @@ void computeSpeciesRatesTest( PARAMS_DATA const & params,
   {
     for( int j = 0; j < numSpecies; ++j )
     {
-      data.speciesRatesDerivatives( i, j ) = data.speciesRatesDerivatives( i, j ) * logmath::exp( -data.speciesConcentration[j] );
+      data.speciesRatesDerivatives( i, j ) = data.speciesRatesDerivatives( i, j ) * logmath::exp( -data.speciesConcentration[j] ) / logmath::dWrtLogConst< double >();
       EXPECT_NEAR( data.speciesRatesDerivatives( i, j ), expectedSpeciesRatesDerivatives[i][j], 1.0e-8 );
     }
   }

src/uncrustify.cfg

@@ -210,7 +210,7 @@ sp_angle_shift                  = force      # ignore/add/remove/force
 sp_permit_cpp11_shift           = false    # false/true
 
 # Add or remove space before '(' of 'if', 'for', 'switch', 'while', etc.
-sp_before_sparen                = force   # ignore/add/remove/force
+sp_before_sparen                = remove   # ignore/add/remove/force
 
 # Add or remove space inside if-condition '(' and ')'.
 sp_inside_sparen                = force  # ignore/add/remove/force