Implementation of equilibrium reactions (#3)

* break up equilbrium reactions to impl file:

* added symmetric matrix type for the equilibrium enforcement when using reaction extents

* add script to take a system and reduce it to RREF for evaluating choice of primary and secondary species

* added printers.hpp to output arrays and matrices easily

* add log based calculation of aggregate concentrations in SpeciesUtilities file

* added aggregate based equilibrium enforcement

* add texlive-full for documentation

Author: rrsettgast

.github/workflows/continuousIntegration.yml

@@ -20,7 +20,7 @@ jobs:
       - name: Extract docker image tag hash
         id: extract_docker_image_tag_hash
         run: |
-          echo "DOCKER_IMAGE_TAG_HASH=310f70ba6065b94858af2cafeba828ea76c02a1d" >> "$GITHUB_OUTPUT"
+          echo "DOCKER_IMAGE_TAG_HASH=eaa7c33eb38fed4df436c1085c6a020ad75434ac" >> "$GITHUB_OUTPUT"
 
   # code_style:
   #   needs: [check_pull_request_is_not_a_draft]

cmake/Macros.cmake

@@ -41,24 +41,30 @@ macro(hpcReact_add_code_checks)
                         --suppress=noExplicitConstructor
                         --suppress=unusedFunction 
                         --suppress=constStatement 
-                        --suppress=unusedStructMember )
+                        --suppress=unusedStructMember
+                        -I../hpcReact/src )
 
-    blt_add_code_checks( PREFIX    ${arg_PREFIX}
-                         SOURCES   ${_sources}
-                         UNCRUSTIFY_CFG_FILE ${PROJECT_SOURCE_DIR}/src/uncrustify.cfg
-                         CPPCHECK_FLAGS ${CPPCHECK_FLAGS}
-                         )
+    if( UNCRUSTIFY_FOUND )
+        blt_add_code_checks( PREFIX    ${arg_PREFIX}
+                            SOURCES   ${_sources}
+                            UNCRUSTIFY_CFG_FILE ${PROJECT_SOURCE_DIR}/src/uncrustify.cfg
+                            CPPCHECK_FLAGS ${CPPCHECK_FLAGS}
+                            )
+        add_test( NAME testUncrustifyCheck
+                    COMMAND bash -c "${CMAKE_MAKE_PROGRAM} uncrustify_check 2> >(tee uncrustify.err) >/dev/null && exit $(cat uncrustify.err | wc -l)"
+                    WORKING_DIRECTORY ${CMAKE_BINARY_DIR} )
+    endif()
 
     if( CPPCHECK_FOUND )
         add_test( NAME testCppCheck
-                COMMAND bash -c "make cppcheck_check 2> >(tee cppcheck.err) >/dev/null && exit $(cat cppcheck.err | wc -l)"
+                COMMAND bash -c "${CMAKE_MAKE_PROGRAM} cppcheck_check 2> >(tee cppcheck.err) >/dev/null && exit $(cat cppcheck.err | wc -l)"
                 WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
                 )
     endif()
 
     if( CLANGTIDY_FOUND )
         add_test( NAME testClangTidy
-                COMMAND bash -c "make clang_tidy_check 2> >(tee tidyCheck.err) >/dev/null && exit $(cat tidyCheck.err | wc -l)"
+                COMMAND bash -c "${CMAKE_MAKE_PROGRAM} clang_tidy_check 2> >(tee tidyCheck.err) >/dev/null && exit $(cat tidyCheck.err | wc -l)"
                 WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
                 )
     endif()

hostconfigs/apple/macOS_base.cmake

@@ -26,7 +26,7 @@ set( DOXYGEN_EXECUTABLE "${HOMEBREW_DIR}/bin/doxygen" CACHE PATH "" FORCE)
 set( CPPCHECK_EXECUTABLE "${HOMEBREW_DIR}/bin/cppcheck" CACHE PATH "" FORCE)
 set( CLANGTIDY_EXECUTABLE "${HOMEBREW_DIR}/opt/llvm/bin/clang-tidy" CACHE PATH "" FORCE)
 
-set( ENABLE_COVERAGE ON CACHE BOOL "" FORCE)
+set( ENABLE_COVERAGE OFF CACHE BOOL "" FORCE)
 set( GCOV_EXECUTABLE "/usr/bin/gcov" CACHE PATH "" FORCE)
 set( LCOV_EXECUTABLE "${HOMEBREW_DIR}/bin/lcov" CACHE PATH "" FORCE)
 

scripts/ci_build_and_test_in_container.sh

@@ -20,9 +20,9 @@ fi
 echo "Using package manager: $PACKAGE_MANAGER"
 
 if [ "$PACKAGE_MANAGER" == "apt" ]; then
-    apt update && apt-get install -y libblas-dev liblapack-dev
+    apt update && apt-get install -y libblas-dev liblapack-dev texlive-full
 elif [ "$PACKAGE_MANAGER" == "yum" ]; then
-    yum install -y blas lapack
+    yum update && yum install -y blas lapack 
 fi
 
 # The or_die function run the passed command line and

src/CMakeLists.txt

@@ -3,9 +3,13 @@ set( hpcReact_headers
      common/macros.hpp
      common/CArrayWrapper.hpp
      reactions/bulkGeneric/EquilibriumReactions.hpp
+     reactions/bulkGeneric/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp
+     reactions/bulkGeneric/EquilibriumReactionsReactionExtents_impl.hpp
      reactions/bulkGeneric/KineticReactions.hpp
+     reactions/bulkGeneric/KineticReactions_impl.hpp
      reactions/bulkGeneric/Parameters.hpp
      reactions/bulkGeneric/ParametersPredefined.hpp
+     reactions/bulkGeneric/SpeciesUtilities.hpp
      )
 
 set( hpcReact_sources
@@ -35,6 +39,7 @@ target_include_directories( hpcReact
                             $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src>
                             $<INSTALL_INTERFACE:include> )
 
+blt_print_target_properties( TARGET hpcReact )
 # install( FILES ${hpcReact_headers} 
 #          DESTINATION include )
 

src/common/DirectSystemSolve.hpp

@@ -1,13 +1,145 @@
 #pragma once
 
 #include "macros.hpp"
+#include "symmetricMatrix.hpp"
 #include <cmath>
 
 namespace hpcReact
 {
+
+template< typename REAL_TYPE, int N >
+bool isPositiveDefinite( REAL_TYPE const (&A)[N][N] )
+{
+  REAL_TYPE temp[N][N];
+
+  // Copy A into temp to avoid modifying original
+  for( int i = 0; i < N; i++ )
+  {
+    for( int j = 0; j < N; j++ )
+    {
+      temp[i][j] = A[i][j];
+    }
+  }
+
+  for( int k = 0; k < N; k++ )
+  {
+    // Compute determinant of k-th leading principal minor using Gaussian elimination
+    if( temp[k][k] <= 0 )
+      return false; // Must be positive
+
+    for( int i = k + 1; i < N; i++ )
+    {
+      REAL_TYPE factor = temp[i][k] / temp[k][k];
+      for( int j = k; j < N; j++ )
+      {
+        temp[i][j] -= factor * temp[k][j];
+      }
+    }
+  }
+  return true;
+}
+
+
+template< typename REAL_TYPE, int N >
+void solveNxN_Cholesky( REAL_TYPE const (&A)[N][N], REAL_TYPE const (&b)[N], REAL_TYPE (& x)[N] )
+{
+  REAL_TYPE L[N][N] = {{0}};
+
+  // **Cholesky Decomposition**
+  for( int i = 0; i < N; i++ )
+  {
+    for( int j = 0; j <= i; j++ )
+    {
+      REAL_TYPE sum = 0;
+      for( int k = 0; k < j; k++ )
+      {
+        sum += L[i][k] * L[j][k];
+      }
+      if( i == j )
+        L[i][j] = sqrt( A[i][i] - sum );
+      else
+        L[i][j] = (A[i][j] - sum) / L[j][j];
+    }
+  }
+
+  // **Forward Substitution: Solve L y = b**
+  //REAL_TYPE y[N];
+  for( int i = 0; i < N; i++ )
+  {
+    x[i] = b[i];
+    for( int j = 0; j < i; j++ )
+      x[i] -= L[i][j] * x[j];
+    x[i] /= L[i][i];
+  }
+
+  // **Backward Substitution: Solve L^T x = y**
+  for( int i = N - 1; i >= 0; i-- )
+  {
+//    x[i] = y[i];
+    for( int j = i + 1; j < N; j++ )
+      x[i] -= L[j][i] * x[j];
+    x[i] /= L[i][i];
+  }
+}
+
+
+
+template< typename REAL_TYPE, int N >
+void solveNxN_Cholesky( symmetricMatrix< REAL_TYPE, int, N > const & A,
+                        REAL_TYPE const (&b)[N],
+                        REAL_TYPE (& x)[N] )
+{
+  symmetricMatrix< REAL_TYPE, int, N > L = {0};
+
+  // **Cholesky Decomposition**
+  for( int i = 0; i < N; i++ )
+  {
+    for( int j = 0; j <= i; j++ )
+    {
+      REAL_TYPE sum = 0;
+      for( int k = 0; k < j; k++ )
+      {
+        sum += L( i, k ) * L( j, k );
+      }
+      if( i == j )
+      {
+        L( i, j ) = sqrt( A( i, i ) - sum );
+      }
+      else
+      {
+        L( i, j ) = (A( i, j ) - sum) / L[j][j];
+      }
+    }
+  }
+
+  // **Forward Substitution: Solve L y = b**
+  REAL_TYPE y[N];
+  for( int i = 0; i < N; i++ )
+  {
+    y[i] = b[i];
+    for( int j = 0; j < i; j++ )
+    {
+      y[i] -= L( i, j ) * y[j];
+    }
+    y[i] /= L( i, i );
+  }
+
+  // **Backward Substitution: Solve L^T x = y**
+  for( int i = N - 1; i >= 0; i-- )
+  {
+    x[i] = y[i];
+    for( int j = i + 1; j < N; j++ )
+    {
+      x[i] -= L( j, i ) * x[j];
+    }
+    x[i] /= L( i, i );
+  }
+}
+
+
 template< typename REAL_TYPE, int N >
 HPCREACT_HOST_DEVICE
-void solveNxN_pivoted( REAL_TYPE A[N][N], REAL_TYPE b[N], REAL_TYPE x[N] )
+void solveNxN_pivoted( REAL_TYPE (& A)[N][N], REAL_TYPE (& b)[N], REAL_TYPE (& x)[N] )
 {
 
 

src/common/macros.hpp

@@ -19,3 +19,48 @@
 /// This macro is used to ignore warnings that that a variable is
 /// unused.
 #define HPCREACT_UNUSED_VAR( ... ) (void)( __VA_ARGS__ )
+
+
+#if defined( __clang__ )
+#define HPCREACT_NO_MISSING_BRACES( ... ) \
+        _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\"")
+#define HPCREACT_NO_MISSING_BRACES_CLOSE \
+        _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")
+
+#define HPCREACT_NO_MISSING_BRACES_OPEN \
+        _Pragma("GCC diagnostic push") \
+        _Pragma("GCC diagnostic ignored \"-Wmissing-braces\"")
+#define HPCREACT_NO_MISSING_BRACES_CLOSE \
+        _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))
+
+#define HPCREACT_NO_MISSING_BRACES_OPEN \
+        __pragma(warning(push)) \
+        __pragma(warning(disable : 4351))
+#define HPCREACT_NO_MISSING_BRACES_CLOSE \
+        __pragma(warning(pop))
+#else
+#define HPCREACT_NO_MISSING_BRACES( ... ) __VA_ARGS__ // No-op for unknown compilers
+#define HPCREACT_NO_MISSING_BRACES_OPEN
+#define HPCREACT_NO_MISSING_BRACES_CLOSE
+#endif

src/common/printers.hpp

@@ -0,0 +1,43 @@
+#pragma once
+
+#include <stdio.h>
+
+template< typename T, int N >
+void print( T const (&a)[N], char const * name, int const numDigits )
+{
+  char format[8];
+  snprintf( format, 8, "%%%d.%de", numDigits+7, numDigits );
+  printf( "%-9s = { ", name );
+  for( int i=0; i<N; ++i )
+  {
+    printf( format, a[i] );
+    if( i < N-1 )
+    {
+      printf( ", " );
+    }
+  }
+  printf( " }\n" );
+}
+
+template< typename T, int NROWS, int NCOLS >
+void print( CArrayWrapper< T, NROWS, NCOLS > const & matrix, char const * name, int const numDigits )
+{
+  char format[8];
+  snprintf( format, 8, "%%%d.%de", numDigits+7, numDigits );
+
+  printf( "%-9s = \n{ \n", name );
+  for( int i=0; i<NROWS; ++i )
+  {
+    printf( " { " );
+    for( int j=0; j<NCOLS; ++j )
+    {
+      printf( format, matrix( i, j ) );
+      if( j < NCOLS-1 )
+      {
+        printf( ", " );
+      }
+    }
+    printf( " },\n" );
+  }
+  printf( " }\n" );
+}