Skip to content

cleanup use of logarithms #19

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rrsettgast wants to merge 10 commits into
base: main
Choose a base branch
from
Open

cleanup use of logarithms #19

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
91e192b
added logMath aliases. replaced log and exp calls everywhere with ali…
rrsettgast Dec 17, 2025
9acef98
fix incorrect derivative correction
rrsettgast Dec 17, 2025
dc141c0
more derivatives
rrsettgast Dec 18, 2025
3ddc42c
uncrustify
rrsettgast Dec 18, 2025
d0ca8f4
uncrustify
rrsettgast Dec 18, 2025
b7a8739
doxygen
rrsettgast Dec 19, 2025
a4724bf
Merge branch 'main' into feature/logCleanup
rrsettgast Dec 19, 2025
dbf6434
added new scaling function for ln/log derivatives
rrsettgast Dec 19, 2025
9909192
franks fixes
rrsettgast Dec 19, 2025
0f6cc11
fix doxygen
rrsettgast Dec 20, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
124 changes: 124 additions & 0 deletions src/common/logMath.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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()); }
}



// ***** 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
7 changes: 1 addition & 6 deletions src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,12 +44,7 @@ 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(),
computeReactionRatesTest< double >( serialAllKineticParams.kineticReactionsParameters(),
initialSpeciesConcentration,
surfaceArea, // No use. Just to pass something here
expectedReactionRates,
Expand Down
21 changes: 3 additions & 18 deletions src/reactions/exampleSystems/unitTests/testKineticReactions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,12 +29,7 @@ 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(),
computeReactionRatesTest< double >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
initialSpeciesConcentration,
surfaceArea, // No use. Just to pass something here
expectedReactionRates,
Expand All @@ -52,12 +47,8 @@ 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(),
computeSpeciesRatesTest< double >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
initialSpeciesConcentration,
expectedSpeciesRates,
expectedSpeciesRatesDerivatives );
Expand All @@ -70,18 +61,12 @@ 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(),
timeStepTest< double >( 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 );
}

int main( int argc, char * * argv )
Expand Down
12 changes: 6 additions & 6 deletions src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,11 +53,11 @@ 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,
Expand All @@ -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] );
}

}
Expand Down
12 changes: 6 additions & 6 deletions src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,11 +54,11 @@ void testMoMasMediumEquilibriumHelper()

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,
Expand All @@ -79,7 +79,7 @@ void testMoMasMediumEquilibriumHelper()

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


Expand Down
16 changes: 8 additions & 8 deletions src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,13 +122,13 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 )

double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] =
{
log( initialPrimarySpeciesConcentration[0] ),
log( initialPrimarySpeciesConcentration[1] ),
log( initialPrimarySpeciesConcentration[2] ),
log( initialPrimarySpeciesConcentration[3] ),
log( initialPrimarySpeciesConcentration[4] ),
log( initialPrimarySpeciesConcentration[5] ),
log( initialPrimarySpeciesConcentration[6] )
logmath::log( initialPrimarySpeciesConcentration[0] ),
logmath::log( initialPrimarySpeciesConcentration[1] ),
logmath::log( initialPrimarySpeciesConcentration[2] ),
logmath::log( initialPrimarySpeciesConcentration[3] ),
logmath::log( initialPrimarySpeciesConcentration[4] ),
logmath::log( initialPrimarySpeciesConcentration[5] ),
logmath::log( initialPrimarySpeciesConcentration[6] )
};

double logPrimarySpeciesConcentration[numPrimarySpecies];
Expand All @@ -150,7 +150,7 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 )

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


Expand Down
32 changes: 8 additions & 24 deletions src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,12 +81,7 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic )

};

computeReactionRatesTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(),
initialSpeciesConcentration,
surfaceArea, // No use. Just to pass something here
expectedReactionRates,
expectedReactionRatesDerivatives );
computeReactionRatesTest< double, true >( carbonateSystemAllKinetic.kineticReactionsParameters(),
computeReactionRatesTest< double >( carbonateSystemAllKinetic.kineticReactionsParameters(),
initialSpeciesConcentration,
surfaceArea, // No use. Just to pass something here
expectedReactionRates,
Expand Down Expand Up @@ -124,12 +119,7 @@ TEST( testKineticReactions, computeReactionRatesQuotientTest_carbonateSystem )
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.0877659574468075e-03, -3.0877659574468075e-03, -2.9999999999999997e-02, 0, 0, 0, 0 }
};

computeReactionRatesTest< double, false >( carbonateSystem.kineticReactionsParameters(),
initialSpeciesConcentration,
surfaceArea,
expectedReactionRates,
expectedReactionRatesDerivatives );
computeReactionRatesTest< double, true >( carbonateSystem.kineticReactionsParameters(),
computeReactionRatesTest< double >( carbonateSystem.kineticReactionsParameters(),
initialSpeciesConcentration,
surfaceArea,
expectedReactionRates,
Expand Down Expand Up @@ -218,18 +208,12 @@ TEST( testKineticReactions, testTimeStep_carbonateSystemAllKinetic )
1.072307827865370e+00 // Na+1
};

timeStepTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(),
10.0,
10000,
initialSpeciesConcentration,
expectedSpeciesConcentrations );

// ln(c) as the primary variable results in a singular system.
// timeStepTest< double, true >( simpleKineticTestRateParams,
// 2.0,
// 10,
// initialSpeciesConcentration,
// expectedSpeciesConcentrations );
timeStepTest< double >( carbonateSystemAllKinetic.kineticReactionsParameters(),
10.0,
10000,
initialSpeciesConcentration,
expectedSpeciesConcentrations );

}

int main( int argc, char * * argv )
Expand Down
2 changes: 1 addition & 1 deletion src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ TEST( testMixedReactions, testTimeStep_carbonateSystem )
1.070434904554991 // Na+1
};

timeStepTest< double, true >( carbonateSystem,
timeStepTest< double >( carbonateSystem,
1.0,
10,
initialAggregateSpeciesConcentration,
Expand Down
Loading
Loading