add initial constitutive relations api

[rrsettgast]

The primary goal of this PR is to enable a constitutive relation interface in HPCReact for computation of activity coefficients, ionic strength, and water activity.

Specific changes include:

• Adds a constitutive activity layer — a new src/constitutive/ module with a common calculateActivities API (activity/activity.hpp) and four models (Identity, DebyeHuckel, Drummond, Bdot). Drummond is for certain neutral species with a different activity coefficient calculation, e.g., CO2.

• Threads ACTIVITY_MODEL as a template parameter through the reaction stack — EquilibriumReactions, KineticReactions, and MixedEquilibriumKineticReactions now take activity parameters alongside reaction parameters, and every example system and unit test is templated on the model.

• Adds water as an explicit participant — activity models return ln(a_w) with derivatives, and the parameter struct carries a waterStoichiometry array so reactions with water (e.g. OH⁻ + H⁺ = H₂O, CO₂ + H₂O = H⁺ + HCO₃⁻) accounts for the water activity term.

• Refactors the secondary species calculation as a self-consistent inner Newton solve — calculateLogSecondarySpeciesConcentration now iterates on log C_sec until concentrations, ionic strength, activity coefficients and a_w are mutually consistent at the given primary concentrations, and d log C_sec/d log C_prim is recovered analytically from the converged constraint.

• Accounts for the activity model in both mass action and kinetic rate evaluation — MassActions converts secondary activities back to concentrations for the mole balance of given aggregate/total concentration, while kinetic rates are evaluated on activities in both the elementary and affinity forms.

• Refactors the reactionRatesUpdateOption integer with a ReactionRateLawOption enum (Elementary vs Affinity) and renames the rate-evaluation entry point to textbook-consistent naming.

• Adds an ideal-solution presolve option (HPCREACT_IDEAL_PRESOLVE) that solves the Identity speciation first (activity coefficient and water activity assumed to be unity) and uses it as the Newton initial guess for advanced models (e.g., BDot), which is what makes the B-dot carbonate cases converge reliably.

• Adds EQ3/6 verification — the geochemistry carbonate system is populated with data0.com.V8.R6 logK values and phreeqc ion-size/b-dot parameters, with EQ3NR/EQ6 input-output decks checked in as reference, new tests for equilibrium, kinetic, and mixed carbonate cases, standalone testBdot / testIonicStrength unit tests, and Python scripts for the reference values of the simple example system.

Note: This PR needs #21 and GEOS-DEV/GEOS#3996 as prerequisites as the activity model uses molality concentration by default as the input, while current GEOS develop considers molarity concentration.

[azibitsker]

The sum of charges needs to be 0. Suggesting to change 1.0 to -1.0

[frankfeifan]

Bdot tests are expected to fail now. Need to update the expected solution.

[frankfeifan]

I've heavily modified and refactored this file to account for the activity model during mass action solve. The root cause is secondary concentrations and activity coefficients are coupled (i.e., ionic strength depends on the secondary concentrations). So we should solve them together, and therefore, concentrations, activity coefficients and activities are mutually consistent.

The specific changes include:

• Refactored calculateLogSecondarySpeciesConcentration and calculateLogSecondarySpeciesConcentrationWrtLogC: an inner Newton solve on log(C_sec) returning concentrations, activity coefficients, activities, a convergence flag, and the exact d log(C_sec)/d log(C_prim).

• Renamed the old fixed-activity-coefficient versions to ...NoActivityUpdate and ...WrtLogCNoActivityUpdate.

• Removed the two calculateAggregatePrimaryConcentrationsWrtLogC overloads and the calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC, which reconstructed secondary concentrations internally under gamma == 1. The two remaining aggregate functions are a pure mole balance, T_i = C_prim,i + sum_j nu_ji C_sec,j, taking the solved secondary arrays as input.

[frankfeifan]

New tests added:

Self-consistent solve, checking mass action holds and a = C * gamma for all species:

• test_calculateLogSecondarySpeciesConcentration_allEquilibrium_identity
• test_calculateLogSecondarySpeciesConcentration_allEquilibrium_bdot
• test_calculateLogSecondarySpeciesConcentration_mixedSystem_identity
• test_calculateLogSecondarySpeciesConcentration_mixedSystem_bdot
• test_calculateLogSecondarySpeciesConcentration_identityActivityModel — reproduces the reference
  concentrations and confirms unit activity coefficients.

Derivative d log(C_sec)/d log(C_prim) against a central difference of the solve itself:

• test_calculateLogSecondarySpeciesConcentrationWrtLogC_allEquilibrium_identity
• test_calculateLogSecondarySpeciesConcentrationWrtLogC_allEquilibrium_bdot
• test_calculateLogSecondarySpeciesConcentrationWrtLogC_mixedSystem_bdot

[codecov]

Codecov Report

❌ Patch coverage is 98.32402% with 12 lines in your changes missing coverage. Please review.
✅ Project coverage is 99.08%. Comparing base (7423e3e) to head (bd34d2f).

Files with missing lines	Patch %	Lines
...stry/unitTests/testGeochemicalKineticReactions.cpp	89.18%	4 Missing ⚠️
src/reactions/massActions/MassActions.hpp	96.87%	4 Missing ⚠️
...eactions/massActions/unitTests/testMassActions.cpp	97.82%	3 Missing ⚠️
...actions/reactionsSystems/KineticReactions_impl.hpp	98.70%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##              main      #17      +/-   ##
===========================================
- Coverage   100.00%   99.08%   -0.92%     
===========================================
  Files           25       34       +9     
  Lines          830     1309     +479     
===========================================
+ Hits           830     1297     +467     
- Misses           0       12      +12     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.