add initial constitutive relations api

Author: rrsettgast

src/CMakeLists.txt

@@ -2,6 +2,7 @@
 set( hpcReact_headers
      common/macros.hpp
      common/CArrayWrapper.hpp
+     constitutive/activity/Bdot.hpp
      reactions/exampleSystems/BulkGeneric.hpp
      reactions/geochemistry/Carbonate.hpp
      reactions/geochemistry/Forge.hpp
@@ -67,10 +68,11 @@ message(STATUS "HPCReact/src CMAKE_CURRENT_SOURCE_DIR: ${CMAKE_CURRENT_SOURCE_DI
 # hpcReact_add_code_checks( PREFIX hpcReact
 #                        EXCLUDES "blt/*" )
 
+add_subdirectory( common/unitTests )
+add_subdirectory( constitutive/unitTests)
 add_subdirectory( reactions/exampleSystems/unitTests )
 add_subdirectory( reactions/geochemistry/unitTests )
 add_subdirectory( reactions/massActions/unitTests )
-add_subdirectory( common/unitTests )
 add_subdirectory( docs )
 
 hpcReact_add_code_checks( PREFIX hpcReact

src/constitutive/activity/Bdot.hpp

@@ -0,0 +1,63 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: (BSD-3-Clause)
+ *
+ * Copyright (c) 2025- Lawrence Livermore National Security LLC
+ * All rights reserved
+ *
+ * See top level LICENSE files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+#pragma once
+
+#include "common/macros.hpp"
+
+namespace hpcReact
+{
+
+template< typename REAL_TYPE,
+          typename INDEX_TYPE,
+          typename IONIC_STRENGTH_TYPE >
+class Bdot
+{
+public:
+  using RealType = REAL_TYPE;
+  using IndexType = INDEX_TYPE;
+
+struct Params : public IONIC_STRENGTH_TYPE::PARAMS
+{
+
+};
+
+
+template< typename ARRAY_1D_TO_CONST,
+          typename ARRAY_1D >
+static inline HPCREACT_HOST_DEVICE 
+void
+calculateActivities( IONIC_STRENGTH_TYPE::PARAMS const & params,
+                     ARRAY_1D_TO_CONST const & speciesConcentrations,
+                     ARRAY_1D & activities )
+{
+
+  RealType const ionicStrength = IONIC_STRENGTH_TYPE::calculate( params, speciesConcentrations );
+  RealType const sqrtI = sqrt(ionicStrength);
+  RealType const A_gamma = 2;
+  RealType const B_gamma = 1.6;
+  auto const & speciesCharge = params.m_speciesCharge;
+  auto const & a = params.m_ionSizeParameter;
+  auto const & b = params.m_bdotParameter;
+
+
+
+  constexpr IndexType numSpecies = params.numSpecies();
+  for( IndexType i=0; i<numSpecies; ++i )
+  {
+    RealType const gamma_coeff = -A_gamma * sqrtI / ( 1 + a[i] * B_gamma * sqrtI );
+    activities[i] = gamma_coeff * speciesCharge[i] * speciesCharge[i] + b[i] * ionicStrength;
+  }
+}
+
+};
+
+
+} // namespace hpcReact

src/constitutive/ionicStrength/SpeciatedIonicStrength.hpp

@@ -0,0 +1,59 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: (BSD-3-Clause)
+ *
+ * Copyright (c) 2025- Lawrence Livermore National Security LLC
+ * All rights reserved
+ *
+ * See top level LICENSE files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+#pragma once
+
+#include "common/CArrayWrapper.hpp"
+#include "common/macros.hpp"
+
+namespace hpcReact
+{
+
+template< typename REAL_TYPE,
+          typename INDEX_TYPE,
+          int NUM_SPECIES >
+class SpeciatedIonicStrength
+{
+public:
+  using RealType = REAL_TYPE;
+  using IndexType = INDEX_TYPE;
+
+struct Params
+{
+  
+  HPCREACT_HOST_DEVICE static constexpr IndexType numSpecies() { return NUM_SPECIES; }
+  HPCREACT_HOST_DEVICE constexpr CArrayWrapper< RealType, NUM_SPECIES > & speciesCharge() { return m_speciesCharge; }
+
+  CArrayWrapper< RealType, NUM_SPECIES > m_speciesCharge;
+};
+
+
+template< typename ARRAY_1D_TO_CONST >
+static inline HPCREACT_HOST_DEVICE 
+REAL_TYPE
+calculate( Params const & params,
+           ARRAY_1D_TO_CONST const & speciesConcentration )
+{
+  REAL_TYPE I = 0.0;
+  auto const &  numSpecies = params.numSpecies();
+  auto const & speciesCharge = params.m_speciesCharge;
+  for( int i=0; i<numSpecies; ++i )
+  {
+    I += speciesConcentration[i] * speciesCharge[i] * speciesCharge[i];
+  }
+  I *= 0.5;
+  return I;
+}
+
+
+};
+
+
+} // namespace hpcReact

src/constitutive/ionicStrength/StoichiometricIonicStrength.hpp

@@ -0,0 +1,42 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: (BSD-3-Clause)
+ *
+ * Copyright (c) 2025- Lawrence Livermore National Security LLC
+ * All rights reserved
+ *
+ * See top level LICENSE files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+#pragma once
+
+#include "common/macros.hpp"
+
+namespace hpcReact
+{
+
+template< typename REAL_TYPE >
+class SpeciatedIonicStrength
+{
+public:
+
+template< typename ARRAY_1D_TO_CONST >
+static inline HPCREACT_HOST_DEVICE 
+REAL_TYPE
+calculate( ARRAY_1D_TO_CONST const & speciesConcentration,
+           ARRAY_1D_TO_CONST const & speciesCharge,
+           int const numSpecies )
+{
+  REAL_TYPE I = 0.0;
+  for( int i=0; i<numSpecies; ++i )
+  {
+    I += speciesConcentration[i] * speciesCharge[i] * speciesCharge[i];
+  }
+  I *= 0.5;
+  return I;
+}
+
+};
+
+
+} // namespace hpcReact

src/constitutive/unitTests/CMakeLists.txt

@@ -0,0 +1,22 @@
+# Specify list of tests
+set( testSourceFiles
+     testBdot.cpp
+     testIonicStrength.cpp )
+
+
+set( dependencyList hpcReact gtest )
+
+if( ENABLE_CUDA )
+    list( APPEND dependencyList cuda )
+endif()
+
+# Add gtest C++ based tests
+foreach(test ${testSourceFiles})
+    get_filename_component( test_name ${test} NAME_WE )
+    blt_add_executable( NAME ${test_name}
+                        SOURCES ${test}
+                        OUTPUT_DIR ${TEST_OUTPUT_DIRECTORY}
+                        DEPENDS_ON ${dependencyList} )
+    blt_add_test( NAME ${test_name}
+                   COMMAND ${test_name} )
+endforeach()

src/constitutive/unitTests/testBdot.cpp

@@ -0,0 +1,49 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: (BSD-3-Clause)
+ *
+ * Copyright (c) 2025- Lawrence Livermore National Security LLC
+ * All rights reserved
+ *
+ * See top level LICENSE files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+
+#include "constitutive/activity/Bdot.hpp"
+#include "reactions/reactionsSystems/Parameters.hpp"
+#include "constitutive/ionicStrength/SpeciatedIonicStrength.hpp"
+#include "common/pmpl.hpp"
+
+#include <gtest/gtest.h>
+
+using namespace hpcReact;
+
+
+
+
+constexpr SpeciatedIonicStrength<double, int>::Params<3> testParams
+{
+  // Species charge
+  { 1.0, -1.0, 2.0 }
+};
+
+TEST( testBdot, testIonicStrength )
+{
+  double speciesConcentration[ testParams.numSpecies() ] = { 0.1, 0.2, 0.3 };
+
+  double I = SpeciatedIonicStrength< double, int >::calculate( testParams,
+                                                               speciesConcentration );
+
+  double expectedI = 0.5 * ( 0.1 * 1.0 * 1.0 + 0.2 * (-1.0) * (-1.0) + 0.3 * 2.0 * 2.0 );
+  EXPECT_DOUBLE_EQ( I, expectedI );
+  
+}
+
+
+int main( int argc, char * * argv )
+{
+  ::testing::InitGoogleTest( &argc, argv );
+  int const result = RUN_ALL_TESTS();
+  return result;
+}

src/constitutive/unitTests/testIonicStrength.cpp

@@ -0,0 +1,48 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: (BSD-3-Clause)
+ *
+ * Copyright (c) 2025- Lawrence Livermore National Security LLC
+ * All rights reserved
+ *
+ * See top level LICENSE files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+
+#include "constitutive/ionicStrength/SpeciatedIonicStrength.hpp"
+#include "common/CArrayWrapper.hpp"
+#include "common/pmpl.hpp"
+
+#include <gtest/gtest.h>
+
+using namespace hpcReact;
+
+
+using IonicStrength = SpeciatedIonicStrength<double, int, 3>;
+
+constexpr IonicStrength::Params testParams
+{
+  // Species charge
+  { 1.0, -1.0, 2.0 }
+};
+
+TEST( testBdot, testIonicStrength )
+{
+  double speciesConcentration[ testParams.numSpecies() ] = { 0.1, 0.2, 0.3 };
+
+  double I = IonicStrength::calculate( testParams,
+                                       speciesConcentration );
+
+  double expectedI = 0.5 * ( 0.1 * 1.0 * 1.0 + 0.2 * (-1.0) * (-1.0) + 0.3 * 2.0 * 2.0 );
+  EXPECT_DOUBLE_EQ( I, expectedI );
+  
+}
+
+
+int main( int argc, char * * argv )
+{
+  ::testing::InitGoogleTest( &argc, argv );
+  int const result = RUN_ALL_TESTS();
+  return result;
+}

src/reactions/reactionsSystems/Parameters.hpp

@@ -116,6 +116,10 @@ struct KineticReactionsParameters
 };
 
 
+
+
+
+
 template< typename REAL_TYPE,
           typename INT_TYPE,
           typename INDEX_TYPE,
@@ -136,12 +140,18 @@ struct MixedReactionsParameters
                                       CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantForward,
                                       CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantReverse,
                                       CArrayWrapper< IntType, NUM_REACTIONS > mobileSecondarySpeciesFlag,
+                                      CArrayWrapper< RealType, NUM_SPECIES > const & speciesCharge,
+                                      CArrayWrapper< RealType, NUM_SPECIES > const & ionSizeParameter,
+                                      CArrayWrapper< RealType, NUM_SPECIES > const & bdotParameter,
                                       IntType const reactionRatesUpdateOption = 1 ):
     m_stoichiometricMatrix( stoichiometricMatrix ),
     m_equilibriumConstant( equilibriumConstant ),
     m_rateConstantForward( rateConstantForward ),
     m_rateConstantReverse( rateConstantReverse ),
     m_mobileSecondarySpeciesFlag( mobileSecondarySpeciesFlag ),
+    m_speciesCharge( speciesCharge ),
+    m_ionSizeParameter( ionSizeParameter ),
+    m_bdotParameter( bdotParameter ),
     m_reactionRatesUpdateOption( reactionRatesUpdateOption )
   {}
 
@@ -250,6 +260,10 @@ struct MixedReactionsParameters
   CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse;
   CArrayWrapper< IntType, NUM_REACTIONS > m_mobileSecondarySpeciesFlag;
 
+  CArrayWrapper< RealType, NUM_SPECIES > m_speciesCharge;
+  CArrayWrapper< RealType, NUM_SPECIES > m_ionSizeParameter;
+  CArrayWrapper< RealType, NUM_SPECIES > m_bdotParameter;
+
   IntType m_reactionRatesUpdateOption; // 0: forward and reverse rate. 1: quotient form.
 };