add initial constitutive relations api

src/CMakeLists.txt

@@ -2,6 +2,8 @@
 set( hpcReact_headers
      common/macros.hpp
      common/CArrayWrapper.hpp
+     constitutive/activity/activity.hpp
+     constitutive/activity/Bdot.hpp
      reactions/exampleSystems/BulkGeneric.hpp
      reactions/geochemistry/Carbonate.hpp
      reactions/geochemistry/Forge.hpp
@@ -67,10 +69,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 )
 
 if( NOT is_submodule )

src/common/constants.hpp

@@ -20,5 +20,10 @@ constexpr double R = 8.31446261815324; // J/(mol K)
 constexpr double F = 96485.3321233100184; // C/mol
 constexpr double NA = 6.02214076e23; // 1/mol
 
+constexpr double metersPerAngstrom = 1.0e-10; // m/Angstrom
+
+constexpr double ln10 = 2.302585092994046e+00;
+constexpr double invln10 = 4.342944819032518e-01;
+
 } // namespace constants
 } // namespace hpcReact

src/common/nonlinearSolvers.hpp

@@ -14,6 +14,7 @@
 #include "macros.hpp"
 #include "DirectSystemSolve.hpp"
 #include <math.h>
+#include <stdio.h>
 
 namespace hpcReact
 {
@@ -170,11 +171,17 @@ bool newtonRaphson( REAL_TYPE (& x)[N],
 
     double const norm = internal::norm< N >( residual );
 
-    printf( "--Iter %d: Residual norm = %.12e\n", iter, norm );
+    if( do_print )
+    {
+      printf( "--Iter %d: Residual norm = %.12e\n", iter, norm ); // LCOV_EXCL_LINE
+    }
 
     if( norm < tol )
     {
-      printf( "--Converged.\n" );
+      if( do_print )
+      {
+        printf( "--Converged.\n" ); // LCOV_EXCL_LINE
+      }
       isConverged = true;
       break;
     }
@@ -190,7 +197,7 @@ bool newtonRaphson( REAL_TYPE (& x)[N],
 
   }
 
-  if( !isConverged )
+  if( !isConverged && do_print )
   {
     printf( "--Newton solver error: Max iterations reached without convergence.\n" ); // LCOV_EXCL_LINE
   }

src/constitutive/activity/Bdot.hpp

@@ -0,0 +1,113 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 "DebyeHuckel.hpp"
+#include "common/CArrayWrapper.hpp"
+#include "common/constants.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
+  {
+    /// Ion size parameter in ANGSTROM (as tabulated by phreeqc.dat).
+    CArrayWrapper< RealType, IONIC_STRENGTH_TYPE::Params::numSpecies() > m_ionSizeParameter;
+
+    /// B-dot parameter in kg/mol, so that b*I is dimensionless.
+    CArrayWrapper< RealType, IONIC_STRENGTH_TYPE::Params::numSpecies() > m_bdotParameter;
+  };
+
+
+
+  /**
+   * @brief Compute ln(gamma) for every species, and its derivatives wrt linear concentration.
+   * @param params activity model parameters
+   * @param speciesConcentrations linear concentrations c_i
+   * @param logActivityCoefficients [out] ln(gamma_i)
+   * @param dLogActivityCoefficients_dConcentrations [out] d ln(gamma_i) / d c_j
+   *
+   * The caller composes the activity as a = c * gamma. Returning gamma rather than the activity
+   * keeps gamma available to callers that need to invert it (e.g. converting a secondary species'
+   * activity back to a concentration for the mole balance).
+   */
+  template< typename ARRAY_1D_TO_CONST,
+            typename ARRAY_1D,
+            typename ARRAY_2D >
+  static inline HPCREACT_HOST_DEVICE
+  void
+  calculateLogActivityCoefficients( Params const & params,
+                                    ARRAY_1D_TO_CONST const & speciesConcentrations,
+                                    ARRAY_1D & logActivityCoefficients,
+                                    ARRAY_2D & dLogActivityCoefficients_dConcentrations )
+  {
+
+    RealType dIonicStrength_dConcentration[ Params::numSpecies() ];
+    RealType const ionicStrength = IONIC_STRENGTH_TYPE::calculate( params,
+                                                                   speciesConcentrations,
+                                                                   dIonicStrength_dConcentration );
+    RealType const sqrtI = sqrt( ionicStrength );
+    RealType const rho_w = 997.0479; // kg/m3
+    RealType const eps_r = 78.54; // dimensionless
+    RealType const T_K = 298.15;
+    RealType const A_gamma = DebyeHuckel< RealType >::A_gamma( T_K, rho_w, eps_r );
+    // A_gamma is returned in its natural-log form, while the log10_gamma equation below is
+    // evaluated in log10. Convert it to the log10 scale.
+    RealType const A_gamma_log10 = A_gamma * constants::invln10;
+
+    // B_gamma*sqrt(I) is an inverse Debye length in 1/m, while m_ionSizeParameter is specified
+    // in Angstrom in the parameter files (e.g. Carbonate.hpp). Scale B_gamma so that the
+    // B*a*sqrt(I) group is dimensionless.
+    RealType const B_gamma = DebyeHuckel< RealType >::B_gamma( T_K, rho_w, eps_r ) * constants::metersPerAngstrom;
+    auto const & speciesCharge = params.m_speciesCharge;
+    auto const & a = params.m_ionSizeParameter;
+    auto const & b = params.m_bdotParameter;
+
+    const IndexType numSpecies = params.numSpecies();
+    for( IndexType i=0; i<numSpecies; ++i )
+    {
+      RealType dlog10_gamma_dI;
+      RealType const DebyeHuckel_term = DebyeHuckel< RealType >::log10_gamma( sqrtI,
+                                                                              speciesCharge[i],
+                                                                              a[i],
+                                                                              A_gamma_log10,
+                                                                              B_gamma,
+                                                                              dlog10_gamma_dI );
+      logActivityCoefficients[i] = ( DebyeHuckel_term + b[i] * ionicStrength ) * constants::ln10;
+
+      // d ln(gamma_i)/dc_j = ln(10) * dlog10(gamma_i)/dI * dI/dc_j.
+      // dlog10_gamma_dI is singular at I = 0, where the ionic strength term is dropped.
+      RealType const dLogGamma_dIonicStrength =
+        ionicStrength > 0.0 ?
+        constants::ln10 * ( dlog10_gamma_dI + b[i] ) :
+        0.0;
+      for( IndexType j=0; j<numSpecies; ++j )
+      {
+        dLogActivityCoefficients_dConcentrations[i][j] = dLogGamma_dIonicStrength * dIonicStrength_dConcentration[j];
+      }
+    }
+  }
+
+};
+
+
+} // namespace hpcReact

src/constitutive/activity/DebyeHuckel.hpp

@@ -0,0 +1,153 @@
+#pragma once
+
+#include "common/constants.hpp"
+#include "common/macros.hpp"
+
+#include <cmath>
+
+/**
+ * @file DebyeHuckel.hpp
+ * @brief Debye–Hückel A^γ and B parameters for aqueous electrolytes.
+ *
+ * This header provides helper functions to compute the Debye–Hückel
+ * parameters A^γ and B in their "native" (natural-log) form for
+ * molal (mol/kg) activity-coefficient models.
+ *
+ * The functions are expressed in terms of fundamental physical constants
+ * and water properties (density and relative permittivity). They can be
+ * used directly in Debye–Hückel or extended Debye–Hückel/B-dot models.
+ */
+
+template< typename REAL_TYPE >
+class DebyeHuckel
+{
+public:
+  using RealType = REAL_TYPE;
+
+  /// π (pi).
+  static constexpr RealType pi        = 3.141592653589793e+00;
+
+  /// Vacuum permittivity ε₀ [F/m].
+  static constexpr RealType e0        = 8.854187812800001e-12;
+
+  /// Elementary charge e [C].
+  static constexpr RealType eChg      = 1.602176634000000e-19;
+
+  /// Boltzmann constant k_B [J/K].
+  static constexpr RealType kB        = 1.380649000000000e-23;
+
+  /// Avogadro constant N_A [1/mol].
+  static constexpr RealType NA        = 6.022140760000000e+23;
+
+
+  // -------------------------------------------------------------
+  // Debye–Hückel A^γ (natural log, molal scale)
+  // -------------------------------------------------------------
+
+  /**
+   * @brief Debye–Hückel A^γ parameter in natural-log form.
+   *
+   * Computes the coefficient A^γ(T,ρ,ε_r) used in the Debye–Hückel
+   * expression for the natural logarithm of the activity coefficient:
+   *
+   * \f[
+   *  \ln \gamma_i =
+   *    - A^\gamma_{\ln}(T,P) \, z_i^2
+   *      \frac{\sqrt{I}}{1 + B(T,P)\, a_i \sqrt{I}}
+   * \f]
+   *
+   * where:
+   *   - \f$ I \f$ is ionic strength in mol/kg (molal),
+   *   - \f$ z_i \f$ is the ionic charge,
+   *   - \f$ a_i \f$ is the ion-size parameter (length),
+   *   - A^γ is independent of the log base (this function is for ln).
+   *
+   * The implementation follows the "native" Debye–Hückel form,
+   * using fundamental physical constants without any 1/ln(10) factors.
+   *
+   * @param T_K          Temperature in kelvin [K].
+   * @param rho_w  Density of water in g/L (≈ kg/m³ numerically).
+   * @param eps_r        Relative permittivity (dielectric constant) of water.
+   * @return A^γ in units consistent with molal ionic strength, for use
+   *         in ln(γ) expressions.
+   */
+  static inline HPCREACT_HOST_DEVICE
+  RealType A_gamma( RealType const T_K,
+                    RealType const rho_w,
+                    RealType const eps_r )
+  {
+    RealType const num = ::pow( eChg, 3.0 ) * ::sqrt( 2.0 * pi * NA * rho_w );
+    RealType const den = ::pow( 4.0 * pi * e0 * eps_r * kB * T_K, 1.5 );
+    return num / den;
+  }
+
+
+  // -------------------------------------------------------------
+  // Debye–Hückel B (natural log, molal scale)
+  // -------------------------------------------------------------
+
+  /**
+   * @brief Debye–Hückel B parameter in natural-log form.
+   *
+   * Computes the Debye–Hückel length-scale parameter B(T,ρ,ε_r) used
+   * in the extended Debye–Hückel law:
+   *
+   * \f[
+   *  \ln \gamma_i =
+   *    - A^\gamma_{\ln}(T,P) \, z_i^2
+   *      \frac{\sqrt{I}}{1 + B(T,P)\, a_i \sqrt{I}} \; ,
+   * \f]
+   *
+   * where:
+   *   - \f$ I \f$ is ionic strength in mol/kg,
+   *   - \f$ a_i \f$ is an ion-size parameter (length).
+   *
+   * The combination \f$ B a_i \sqrt{I} \f$ is dimensionless; the
+   * absolute units of B therefore depend on the length units chosen
+   * for \f$ a_i \f$.
+   *
+   * @param T_K          Temperature in kelvin [K].
+   * @param rho_w  Density of water in g/L (≈ kg/m³ numerically).
+   * @param eps_r        Relative permittivity (dielectric constant) of water.
+   * @return B parameter for use in ln(γ) expressions.
+   */
+  static inline HPCREACT_HOST_DEVICE
+  RealType B_gamma( RealType const T_K,
+                    RealType const rho_w,
+                    RealType const eps_r )
+  {
+    RealType const num = 2.0 * NA * rho_w * eChg * eChg;
+    RealType const den = e0 * eps_r * kB * T_K;
+    return ::sqrt( num / den );
+  }
+
+
+  static inline HPCREACT_HOST_DEVICE
+  RealType log10_gamma( RealType const sqrtI,
+                        RealType const zi,
+                        RealType const ai,
+                        RealType const T_K,
+                        RealType & dlog10_gamma_dI )
+  {
+    RealType const A       = A_gamma( T_K );
+    RealType const B       = B_gamma( T_K );
+    RealType const denom   = 1 + B * ai * sqrtI;
+    dlog10_gamma_dI = -0.5 * A * zi * zi / ( sqrtI * denom * denom );
+    return -A * zi * zi * sqrtI / denom;
+  }
+
+
+  static inline HPCREACT_HOST_DEVICE
+  RealType log10_gamma( RealType const sqrtI,
+                        RealType const zi,
+                        RealType const ai,
+                        RealType const A,
+                        RealType const B,
+                        RealType & dlog10_gamma_dI )
+  {
+    RealType const denom   = 1 + B * ai * sqrtI;
+    dlog10_gamma_dI = -0.5 * A * zi * zi / ( sqrtI * denom * denom );
+    return -A * zi * zi * sqrtI / denom;
+  }
+
+};

src/constitutive/activity/Identity.hpp

@@ -0,0 +1,60 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 Identity
+{
+public:
+  using RealType = REAL_TYPE;
+  using IndexType = INDEX_TYPE;
+
+  struct Params : public IONIC_STRENGTH_TYPE::Params
+  {};
+
+  /**
+   * @brief Ideal solution: gamma = 1 for every species, so ln(gamma) = 0 and all derivatives vanish.
+   */
+  template< typename ARRAY_1D_TO_CONST,
+            typename ARRAY_1D,
+            typename ARRAY_2D,
+            typename PARAMS >
+  static inline HPCREACT_HOST_DEVICE
+  void
+  calculateLogActivityCoefficients( PARAMS const &,
+                                    ARRAY_1D_TO_CONST const & speciesConcentrations,
+                                    ARRAY_1D & logActivityCoefficients,
+                                    ARRAY_2D & dLogActivityCoefficients_dConcentrations )
+  {
+    HPCREACT_UNUSED_VAR( speciesConcentrations );
+
+    constexpr IndexType numSpecies = PARAMS::numSpecies();
+    for( IndexType i=0; i<numSpecies; ++i )
+    {
+      logActivityCoefficients[i] = 0.0;
+      for( IndexType j=0; j<numSpecies; ++j )
+      {
+        dLogActivityCoefficients_dConcentrations[i][j] = 0.0;
+      }
+    }
+  }
+
+};
+
+
+} // namespace hpcReact