From 4a7a3101b61ec264ce2364f6a1836e5b0f338b43 Mon Sep 17 00:00:00 2001 From: Alex Zibitsker Date: Tue, 25 Nov 2025 15:37:23 -0800 Subject: [PATCH 1/4] feat: updated the forge experiment chemical system using the data prepared by Chaoyi --- src/reactions/geochemistry/Forge.hpp | 143 ++++++++++++++++----------- 1 file changed, 84 insertions(+), 59 deletions(-) diff --git a/src/reactions/geochemistry/Forge.hpp b/src/reactions/geochemistry/Forge.hpp index 276577c..5fe288b 100644 --- a/src/reactions/geochemistry/Forge.hpp +++ b/src/reactions/geochemistry/Forge.hpp @@ -22,57 +22,67 @@ namespace geochemistry namespace forge { - // Stoichiometric matrix [13 reactions × 23 species] - // Columns 0–12: secondary species (must be -1 on diagonal) - // Columns 13–22: primary species -constexpr CArrayWrapper< double, 19, 26 > soichMatrix = -{// CaCO₃ CaHCO₃⁺ CaSO₄ CaCl⁺ CaCl₂ MgHCO₃⁺ MgCO₃ MgCl⁺ CO₂(aq) HSO₄⁻ KHSO₄ HSiO₃⁻ NaHSilO₃ NaCl KCl KSO₄⁻ | H⁺ Ca²⁺ Mg²⁺ Na⁺ K⁺ Al³⁺ HCO₃⁻ SO₄²⁻ Cl⁻ SiO₂(aq) - { -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 1, 0, 0, 0 }, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻ - { 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0 }, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻ - { 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0 }, // CaSO₄ ⇌ Ca²⁺ + SO₄²⁻ - { 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0 }, // CaCl⁺ ⇌ Ca²⁺ + Cl⁻ - { 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0 }, // CaCl₂ ⇌ Ca²⁺ + 2Cl⁻ - { 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0 }, // MgHCO₃⁺ ⇌ Mg²⁺ + HCO₃⁻ - { 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 1, 0, 0, 0 }, // MgCO₃(aq) + H⁺⇌ Mg²⁺ + HCO₃⁻ - { 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0 }, // MgCl⁺ ⇌ Mg²⁺ + Cl⁻ - { 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // CO₂(aq) + H₂O ⇌ H⁺ + HCO₃⁻ - { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // HSO₄⁻ ⇌ H⁺ + SO₄²⁻ - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0 }, // KHSO₄ ⇌ H⁺ + K⁺ + SO₄²⁻ - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // HSiO₃⁻ ⇌ H⁺ + SiO₂(aq) - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1 }, // NaHSilO₃ ⇌ H⁺ + Na⁺ + SiO₂(aq) - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0 }, // NaCl ⇌ Na⁺ + Cl⁻ - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0 }, // KCl ⇌ K⁺ + Cl⁻ - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0 }, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻ - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1, 1, 0, 0, 0, 2, 0, 0, 0 }, // Dolomite: CaMg(CO₃)₂(s) + 2H⁺ ⇌ Ca²⁺ + Mg²⁺ + 2HCO₃⁻ - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, 0, 1, 1, 0, 0, 0, 3 }, // Microcline: KAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + K⁺ + 3SiO₂(aq) - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, 1, 0, 1, 0, 1, 0, 3 } // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq) + // Stoichiometric matrix [24 reactions × 29 species] + // Columns 1–19: secondary species (must be -1 on diagonal) + // Columns 20–29: primary species +constexpr CArrayWrapper< double, 24, 29 > stoichMatrix = +{// CaCO₃ CaHCO₃⁺ CaSO₄ CaCl⁺ CaCl₂ MgHCO₃⁺ MgCO₃ MgCl⁺ CO₂(aq) HSO₄⁻ KHSO₄ HSiO₃⁻ NaHSilO₃ NaCl KCl KSO₄⁻ AlOH²⁺ Al(OH)₂ OH⁻ | H⁺ Ca²⁺ Mg²⁺ Na⁺ K⁺ Al³⁺ HCO₃⁻ SO₄²⁻ Cl⁻ SiO₂(aq) + { -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 1, 0, 0, 0 }, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻ + { 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0 }, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻ + { 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0 }, // CaSO₄(aq) ⇌ Ca²⁺ + SO₄²⁻ + { 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0 }, // CaCl⁺ ⇌ Ca²⁺ + Cl⁻ + { 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0 }, // CaCl₂(aq) ⇌ Ca²⁺ + 2Cl⁻ + { 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0 }, // MgHCO₃⁺ ⇌ Mg²⁺ + HCO₃⁻ + { 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 1, 0, 0, 0 }, // MgCO₃(aq) + H⁺⇌ Mg²⁺ + HCO₃⁻ + { 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0 }, // MgCl⁺ ⇌ Mg²⁺ + Cl⁻ + { 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // CO₂(aq) + H₂O ⇌ H⁺ + HCO₃⁻ + { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // HSO₄⁻ ⇌ H⁺ + SO₄²⁻ + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0 }, // KHSO₄(aq) ⇌ H⁺ + K⁺ + SO₄²⁻ + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // HSiO₃⁻ ⇌ H⁺ + SiO₂(aq) + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1 }, // NaHSiO₃(aq) ⇌ H⁺ + Na⁺ + SiO₂(aq) + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0 }, // NaCl(aq) ⇌ Na⁺ + Cl⁻ + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0 }, // KCl(aq) ⇌ K⁺ + Cl⁻ + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0 }, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻ + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // AlOH²⁺ + H⁺ ⇌ Al³⁺ + H₂O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -2, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // Al(OH)₂⁺ + 2H⁺ ⇌ Al³⁺ + 2H₂O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // OH⁻ + H⁺ ⇌ H₂O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // SiO2(s) ⇌ SiO2(aq) + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -10, 0, 3, 0, 1, 1, 0, 0, 0, 3 }, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -8, 1, 0, 0, 0, 2, 0, 0, 0, 2 }, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, 0, 1, 1, 0, 0, 0, 3 }, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0, 0, 0, 0, 1, 0, 0, 0, 2 } // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O }; // Must convert these. They should not be the log. -constexpr CArrayWrapper< double, 19 > equilibriumConstants = +constexpr CArrayWrapper< double, 24 > equilibriumConstants = { - 5.9636, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻ - -1.4181, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻ - -2.5111, // CaSO₄ ⇌ Ca²⁺ + SO₄²⁻ - 0.3811, // CaCl⁺ ⇌ Ca²⁺ + Cl⁻ - 0.3811, // CaCl₂ ⇌ Ca²⁺ + 2Cl⁻ (approximate, same source) - -1.4355, // MgHCO₃⁺ ⇌ Mg²⁺ + HCO₃⁻ - -6.5632, // MgCO₃(aq) + H⁺ ⇌ Mg²⁺ + HCO₃⁻ - -0.1820, // MgCl⁺ ⇌ Mg²⁺ + Cl⁻ - -6.3882, // CO₂(aq) + H₂O ⇌ H⁺ + HCO₃⁻ - -3.0020, // HSO₄⁻ ⇌ H⁺ + SO₄²⁻ - -2.2935, // KHSO₄ ⇌ H⁺ + K⁺ + SO₄²⁻ - -9.0844, // HSiO₃⁻ ⇌ H⁺ + SiO₂(aq) - -7.8291, // NaHSilO₃ ⇌ H⁺ + Na⁺ + SiO₂(aq) - 0.4730, // NaCl ⇌ Na⁺ + Cl⁻ - 0.9240, // KCl ⇌ K⁺ + Cl⁻ - -1.1946, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻ - 0.0944, // Dolomite: CaMg(CO₃)₂(s) + 2H⁺ ⇌ Ca²⁺ + Mg²⁺ + 2HCO₃⁻ - -1.8683, // Microcline: KAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + K⁺ + 3SiO₂(aq) - 0.2236 // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq) + 104351.8133, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻ + 3.98107E-03, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻ + 3.6915E-04, // CaSO₄ ⇌ Ca²⁺ + SO₄²⁻ + 0.2930, // CaCl⁺ ⇌ Ca²⁺ + Cl⁻ + 0.1228, // CaCl₂ ⇌ Ca²⁺ + 2Cl⁻ (approximate, same source) + 0.0038, // MgHCO₃⁺ ⇌ Mg²⁺ + HCO₃⁻ + 7.4559E+05, // MgCO₃(aq) + H⁺ ⇌ Mg²⁺ + HCO₃⁻ + 0.0726, // MgCl⁺ ⇌ Mg²⁺ + Cl⁻ + 6.3548E-08, // CO₂(aq) + H₂O ⇌ H⁺ + HCO₃⁻ + 3.4017E-05, // HSO₄⁻ ⇌ H⁺ + SO₄²⁻ + 3.7068E-05, // KHSO₄ ⇌ H⁺ + K⁺ + SO₄²⁻ + 6.9088E+08, // HSiO₃⁻ ⇌ H⁺ + SiO₂(aq) + 4.5520E+07, // NaHSiO₃ ⇌ H⁺ + Na⁺ + SiO₂(aq) + 0.8067, // NaCl ⇌ Na⁺ + Cl⁻ + 1.6398, // KCl ⇌ K⁺ + Cl⁻ + 0.0121, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻ + 14.6994, // AlOH²⁺ + H⁺ ⇌ Al³⁺ + H₂O + 3.8530E+03, // Al(OH)₂⁺ + 2H⁺ ⇌ Al³⁺ + 2H₂O + 1.9213E+11, // OH⁻ + H⁺ ⇌ H₂O + 0.0036, // SiO2(s) ⇌ SiO2(aq) + 6.6130E+15, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O + 1.3047E+05, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O + 1.7669E-04, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O + 3.1463E-05 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O }; -constexpr CArrayWrapper< double, 19 > fwRateConstant = +constexpr CArrayWrapper< double, 24 > fwRateConstant = { 0.0, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻ 0.0, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻ @@ -90,13 +100,18 @@ constexpr CArrayWrapper< double, 19 > fwRateConstant = 0.0, // NaCl ⇌ Na⁺ + Cl⁻ 0.0, // KCl ⇌ K⁺ + Cl⁻ 0.0, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻ - 1.0, // Dolomite: CaMg(CO₃)₂(s) + 2H⁺ ⇌ Ca²⁺ + Mg²⁺ + 2HCO₃⁻ - 1.0, // Microcline: KAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + K⁺ + 3SiO₂(aq) - 1.0 // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq) + 0.0, // AlOH²⁺ + H⁺ ⇌ Al³⁺ + H₂O + 0.0, // Al(OH)₂⁺ + 2H⁺ ⇌ Al³⁺ + 2H₂O + 0.0, // OH⁻ + H⁺ ⇌ H₂O + 2.7043E-08, // SiO2(s) ⇌ SiO2(aq) + 3.0145E-11, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O + 1.0801E-08, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O + 1.1283e-10, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O + 1.8137e-12 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O }; -constexpr CArrayWrapper< double, 19 > reverseRateConstant = +constexpr CArrayWrapper< double, 24 > reverseRateConstant = { 0.0, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻ 0.0, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻ @@ -114,12 +129,17 @@ constexpr CArrayWrapper< double, 19 > reverseRateConstant = 0.0, // NaCl ⇌ Na⁺ + Cl⁻ 0.0, // KCl ⇌ K⁺ + Cl⁻ 0.0, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻ - 1.0, // Dolomite: CaMg(CO₃)₂(s) + 2H⁺ ⇌ Ca²⁺ + Mg²⁺ + 2HCO₃⁻ - 1.0, // Microcline: KAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + K⁺ + 3SiO₂(aq) - 1.0 // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq) + 0.0, // AlOH²⁺ + H⁺ ⇌ Al³⁺ + H₂O + 0.0, // Al(OH)₂⁺ + 2H⁺ ⇌ Al³⁺ + 2H₂O + 0.0, // OH⁻ + H⁺ ⇌ H₂O + 7.5119E-06, // SiO2(s) ⇌ SiO2(aq) + 4.5584E-27, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O + 8.2785E-14, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O + 6.3858E-07, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O + 5.7645E-08 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O }; -constexpr CArrayWrapper< int, 19 > mobileSpeciesFlag = +constexpr CArrayWrapper< int, 24 > mobileSpeciesFlag = { 1, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻ 1, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻ @@ -137,19 +157,24 @@ constexpr CArrayWrapper< int, 19 > mobileSpeciesFlag = 1, // NaCl ⇌ Na⁺ + Cl⁻ 1, // KCl ⇌ K⁺ + Cl⁻ 1, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻ - 1, // Dolomite: CaMg(CO₃)₂(s) + 2H⁺ ⇌ Ca²⁺ + Mg²⁺ + 2HCO₃⁻ - 1, // Microcline: KAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + K⁺ + 3SiO₂(aq) - 1 // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq) + 1, // AlOH²⁺ + H⁺ ⇌ Al³⁺ + H₂O + 1, // Al(OH)₂⁺ + 2H⁺ ⇌ Al³⁺ + 2H₂O + 1, // OH⁻ + H⁺ ⇌ H₂O + 1, // SiO2(s) ⇌ SiO2(aq) + 1, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O + 1, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O + 1, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O + 1 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O }; } -using forgeSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 26, 19, 16 >; +using forgeSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 29, 24, 19 >; -constexpr forgeSystemType forgeSystem( forge::soichMatrix, forge::equilibriumConstants, forge::fwRateConstant, forge::reverseRateConstant, forge::mobileSpeciesFlag ); +constexpr forgeSystemType forgeSystem( forge::stoichMatrix, forge::equilibriumConstants, forge::fwRateConstant, forge::reverseRateConstant, forge::mobileSpeciesFlag ); // *****UNCRUSTIFY-ON****** } -} +} \ No newline at end of file From 3aac5c670eb582ee4230966652c3a8e4019b099f Mon Sep 17 00:00:00 2001 From: Alex Zibitsker Date: Wed, 3 Dec 2025 11:02:57 -0800 Subject: [PATCH 2/4] Added a unit test for the FORGE reaction system, currently just to test the convergence and not validate against a known solution --- .../geochemistry/GeochemicalSystems.hpp | 4 +- .../geochemistry/unitTests/CMakeLists.txt | 1 + .../unitTests/testForgeReactions.cpp | 73 +++++++++++++++++++ .../mixedReactionsTestUtilities.hpp | 2 +- 4 files changed, 77 insertions(+), 3 deletions(-) create mode 100644 src/reactions/geochemistry/unitTests/testForgeReactions.cpp diff --git a/src/reactions/geochemistry/GeochemicalSystems.hpp b/src/reactions/geochemistry/GeochemicalSystems.hpp index 00c0731..c9ae53f 100644 --- a/src/reactions/geochemistry/GeochemicalSystems.hpp +++ b/src/reactions/geochemistry/GeochemicalSystems.hpp @@ -22,9 +22,9 @@ namespace hpcReact namespace geochemistry { -using systemTypes = std::variant< ultramaficSystemType, - carbonateSystemType, +using systemTypes = std::variant< carbonateSystemType, carbonateSystemAllEquilibriumType, + ultramaficSystemType, forgeSystemType >; } // namespace geochemistry diff --git a/src/reactions/geochemistry/unitTests/CMakeLists.txt b/src/reactions/geochemistry/unitTests/CMakeLists.txt index 69e5522..b0e1ccc 100644 --- a/src/reactions/geochemistry/unitTests/CMakeLists.txt +++ b/src/reactions/geochemistry/unitTests/CMakeLists.txt @@ -3,6 +3,7 @@ set( testSourceFiles testGeochemicalEquilibriumReactions.cpp testGeochemicalKineticReactions.cpp testGeochemicalMixedReactions.cpp + testForgeReactions.cpp ) set( dependencyList hpcReact gtest ) diff --git a/src/reactions/geochemistry/unitTests/testForgeReactions.cpp b/src/reactions/geochemistry/unitTests/testForgeReactions.cpp new file mode 100644 index 0000000..bda6435 --- /dev/null +++ b/src/reactions/geochemistry/unitTests/testForgeReactions.cpp @@ -0,0 +1,73 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp" +#include "../GeochemicalSystems.hpp" + + +using namespace hpcReact; +using namespace hpcReact::unitTest_utilities; + + +TEST( testForgeReactions, testTimeStep_forgeSystem ) +{ + using namespace hpcReact::geochemistry; + + static constexpr int numPrimarySpecies = forgeSystemType::numPrimarySpecies(); + + double const surfaceArea[forgeSystemType::numKineticReactions()] = + { + 1.0, 1.0, 1.0, 1.0, 1.0 // SiO2(s), KAlMg3Si3O10(OH)2(s), CaAl2(SiO4)2(s), KAlSi3O8(s), Al2Si2O5(OH)4(s) + }; + + double const initialAggregateSpeciesConcentration[numPrimarySpecies] = + { + 5e-08, // H+ + 8.3991e-0, // Ca+2 + 0.00014146, // Mg+2 + 0.0034061, // Na+ + 0.00020949, // K+ + 2.3494e-08, // Al+++ + 0.0016047, // HCO3- + 0.00038069, // SO4-2 + 0.0010018, // Cl- + 4.4047e-06 //SiO2(aq) + }; + + double const expectedSpeciesConcentrations[numPrimarySpecies] = + { + 5e-08, // H+ + 8.3991e-0, // Ca+2 + 0.00014146, // Mg+2 + 0.0034061, // Na+ + 0.00020949, // K+ + 2.3494e-08, // Al+++ + 0.0016047, // HCO3- + 0.00038069, // SO4-2 + 0.0010018, // Cl- + 4.4047e-06 //SiO2(aq) + }; + + timeStepTest< double, true >( forgeSystem, + 1.e-3, + 10, + initialAggregateSpeciesConcentration, + surfaceArea, + expectedSpeciesConcentrations ); + +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + int const result = RUN_ALL_TESTS(); + return result; +} diff --git a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp index 7426c03..2227e53 100644 --- a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp @@ -130,7 +130,7 @@ void timeStepTest( PARAMS_DATA const & params, } }; - nonlinearSolvers::newtonRaphson< numPrimarySpecies >( logPrimarySpeciesConcentration, computeResidualAndJacobian ); + nonlinearSolvers::newtonRaphson< numPrimarySpecies >( logPrimarySpeciesConcentration, computeResidualAndJacobian, 25 ); time += dt; } From 181359fb2a64144fedc3d6ff77937827ea4f401c Mon Sep 17 00:00:00 2001 From: Alex Zibitsker Date: Fri, 12 Dec 2025 08:32:08 -0800 Subject: [PATCH 3/4] Updated reaction constans (EQ and kb) to use molarity-based units until activity model is implemented --- src/reactions/geochemistry/Forge.hpp | 48 +++++++++---------- .../unitTests/testForgeReactions.cpp | 42 ++++++++-------- .../mixedReactionsTestUtilities.hpp | 2 +- 3 files changed, 46 insertions(+), 46 deletions(-) diff --git a/src/reactions/geochemistry/Forge.hpp b/src/reactions/geochemistry/Forge.hpp index 5fe288b..95bee9a 100644 --- a/src/reactions/geochemistry/Forge.hpp +++ b/src/reactions/geochemistry/Forge.hpp @@ -57,24 +57,24 @@ constexpr CArrayWrapper< double, 24, 29 > stoichMatrix = constexpr CArrayWrapper< double, 24 > equilibriumConstants = { 104351.8133, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻ - 3.98107E-03, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻ - 3.6915E-04, // CaSO₄ ⇌ Ca²⁺ + SO₄²⁻ - 0.2930, // CaCl⁺ ⇌ Ca²⁺ + Cl⁻ - 0.1228, // CaCl₂ ⇌ Ca²⁺ + 2Cl⁻ (approximate, same source) - 0.0038, // MgHCO₃⁺ ⇌ Mg²⁺ + HCO₃⁻ + 3.9739, //3.98107E-03, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻ + 0.3685, //3.6915E-04, // CaSO₄ ⇌ Ca²⁺ + SO₄²⁻ + 292.4726, //0.2930 // CaCl⁺ ⇌ Ca²⁺ + Cl⁻ + 1.2236e+05, //0.1228 // CaCl₂ ⇌ Ca²⁺ + 2Cl⁻ (approximate, same source) + 3.7932, //0.0038 // MgHCO₃⁺ ⇌ Mg²⁺ + HCO₃⁻ 7.4559E+05, // MgCO₃(aq) + H⁺ ⇌ Mg²⁺ + HCO₃⁻ - 0.0726, // MgCl⁺ ⇌ Mg²⁺ + Cl⁻ - 6.3548E-08, // CO₂(aq) + H₂O ⇌ H⁺ + HCO₃⁻ - 3.4017E-05, // HSO₄⁻ ⇌ H⁺ + SO₄²⁻ - 3.7068E-05, // KHSO₄ ⇌ H⁺ + K⁺ + SO₄²⁻ - 6.9088E+08, // HSiO₃⁻ ⇌ H⁺ + SiO₂(aq) - 4.5520E+07, // NaHSiO₃ ⇌ H⁺ + Na⁺ + SiO₂(aq) - 0.8067, // NaCl ⇌ Na⁺ + Cl⁻ - 1.6398, // KCl ⇌ K⁺ + Cl⁻ - 0.0121, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻ - 14.6994, // AlOH²⁺ + H⁺ ⇌ Al³⁺ + H₂O - 3.8530E+03, // Al(OH)₂⁺ + 2H⁺ ⇌ Al³⁺ + 2H₂O - 1.9213E+11, // OH⁻ + H⁺ ⇌ H₂O + 72.4693, //0.0726 // MgCl⁺ ⇌ Mg²⁺ + Cl⁻ + 6.3434e-05, //6.3548E-08 // CO₂(aq) + H₂O ⇌ H⁺ + HCO₃⁻ + 0.0340, //3.4017E-05, // HSO₄⁻ ⇌ H⁺ + SO₄²⁻ + 36.9347, //3.7068E-05, // KHSO₄ ⇌ H⁺ + K⁺ + SO₄²⁻ + 6.8964e+11, // 6.9088E+08, // HSiO₃⁻ ⇌ H⁺ + SiO₂(aq) + 4.5356e+13, //4.5520E+07, // NaHSiO₃ ⇌ H⁺ + Na⁺ + SiO₂(aq) + 805.2479, //0.8067, // NaCl ⇌ Na⁺ + Cl⁻ + 1.6368e+03, //1.6398, // KCl ⇌ K⁺ + Cl⁻ + 12.0782, //0.0121, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻ + 0.0147, //14.6994, // AlOH²⁺ + H⁺ ⇌ Al³⁺ + H₂O + 0.0039, //3.8530E+03, // Al(OH)₂⁺ + 2H⁺ ⇌ Al³⁺ + 2H₂O + 1.9282e+05, //1.9213E+11, // OH⁻ + H⁺ ⇌ H₂O 0.0036, // SiO2(s) ⇌ SiO2(aq) 6.6130E+15, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O 1.3047E+05, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O @@ -106,8 +106,8 @@ constexpr CArrayWrapper< double, 24 > fwRateConstant = 2.7043E-08, // SiO2(s) ⇌ SiO2(aq) 3.0145E-11, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O 1.0801E-08, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O - 1.1283e-10, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O - 1.8137e-12 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O + 1.1283e-10, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O + 1.8137e-12 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O }; @@ -132,11 +132,11 @@ constexpr CArrayWrapper< double, 24 > reverseRateConstant = 0.0, // AlOH²⁺ + H⁺ ⇌ Al³⁺ + H₂O 0.0, // Al(OH)₂⁺ + 2H⁺ ⇌ Al³⁺ + 2H₂O 0.0, // OH⁻ + H⁺ ⇌ H₂O - 7.5119E-06, // SiO2(s) ⇌ SiO2(aq) - 4.5584E-27, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O - 8.2785E-14, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O - 6.3858E-07, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O - 5.7645E-08 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O + 7.4240e-09, //7.5119E-06, // SiO2(s) ⇌ SiO2(aq) + 4.5420e-21, //4.5584E-27, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O + 8.2341e-05, //8.2785E-14, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O + 6.3975e-10, //6.3858E-07, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O + 57.3348 //5.7645E-08 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O }; constexpr CArrayWrapper< int, 24 > mobileSpeciesFlag = diff --git a/src/reactions/geochemistry/unitTests/testForgeReactions.cpp b/src/reactions/geochemistry/unitTests/testForgeReactions.cpp index bda6435..4fda1c4 100644 --- a/src/reactions/geochemistry/unitTests/testForgeReactions.cpp +++ b/src/reactions/geochemistry/unitTests/testForgeReactions.cpp @@ -30,34 +30,34 @@ TEST( testForgeReactions, testTimeStep_forgeSystem ) double const initialAggregateSpeciesConcentration[numPrimarySpecies] = { - 5e-08, // H+ - 8.3991e-0, // Ca+2 - 0.00014146, // Mg+2 - 0.0034061, // Na+ - 0.00020949, // K+ - 2.3494e-08, // Al+++ - 0.0016047, // HCO3- - 0.00038069, // SO4-2 - 0.0010018, // Cl- - 4.4047e-06 //SiO2(aq) + 1.9964e-04, // H+ + 0.1050, // Ca+2 + 0.1209, // Mg+2 + 3.4000, // Na+ + 0.1250, // K+ + 5.4367e-04, // Al+++ + 2.60, // HCO3- + 0.380, // SO4-2 + 1.0, // Cl- + 0.0103 //SiO2(aq) }; double const expectedSpeciesConcentrations[numPrimarySpecies] = { - 5e-08, // H+ - 8.3991e-0, // Ca+2 - 0.00014146, // Mg+2 - 0.0034061, // Na+ - 0.00020949, // K+ - 2.3494e-08, // Al+++ - 0.0016047, // HCO3- - 0.00038069, // SO4-2 - 0.0010018, // Cl- - 4.4047e-06 //SiO2(aq) + 1.9964e-04, // H+ + 0.1050, // Ca+2 + 0.1209, // Mg+2 + 3.4000, // Na+ + 0.1250, // K+ + 5.4367e-04, // Al+++ + 2.60, // HCO3- + 0.380, // SO4-2 + 1.0, // Cl- + 0.0103 //SiO2(aq) }; timeStepTest< double, true >( forgeSystem, - 1.e-3, + 1.e-5, 10, initialAggregateSpeciesConcentration, surfaceArea, diff --git a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp index 2227e53..96c9af5 100644 --- a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp @@ -130,7 +130,7 @@ void timeStepTest( PARAMS_DATA const & params, } }; - nonlinearSolvers::newtonRaphson< numPrimarySpecies >( logPrimarySpeciesConcentration, computeResidualAndJacobian, 25 ); + nonlinearSolvers::newtonRaphson< numPrimarySpecies >( logPrimarySpeciesConcentration, computeResidualAndJacobian, 50 ); time += dt; } From 0672f1c5c3dfd4e6de70d2b1ce366d46a2e847cf Mon Sep 17 00:00:00 2001 From: Alex Zibitsker Date: Fri, 12 Dec 2025 14:51:17 -0800 Subject: [PATCH 4/4] Fixed bug in stoicheometry coeff in the last kinetic reaction and updated the rates --- src/reactions/geochemistry/Forge.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/reactions/geochemistry/Forge.hpp b/src/reactions/geochemistry/Forge.hpp index 95bee9a..24a8a6e 100644 --- a/src/reactions/geochemistry/Forge.hpp +++ b/src/reactions/geochemistry/Forge.hpp @@ -50,7 +50,7 @@ constexpr CArrayWrapper< double, 24, 29 > stoichMatrix = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -10, 0, 3, 0, 1, 1, 0, 0, 0, 3 }, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -8, 1, 0, 0, 0, 2, 0, 0, 0, 2 }, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, 0, 1, 1, 0, 0, 0, 3 }, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0, 0, 0, 0, 1, 0, 0, 0, 2 } // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0, 0, 0, 0, 2, 0, 0, 0, 2 } // Al2Si2O5(OH)4(s) + 6H+ ⇌ 2Al3+ + 2 SiO2 + 5H2O }; // Must convert these. They should not be the log. @@ -79,7 +79,7 @@ constexpr CArrayWrapper< double, 24 > equilibriumConstants = 6.6130E+15, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O 1.3047E+05, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O 1.7669E-04, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O - 3.1463E-05 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O + 3.1463E-05 // Al2Si2O5(OH)4(s) + 6H+ ⇌ 2Al3+ + 2 SiO2 + 5H2O }; constexpr CArrayWrapper< double, 24 > fwRateConstant = @@ -107,7 +107,7 @@ constexpr CArrayWrapper< double, 24 > fwRateConstant = 3.0145E-11, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O 1.0801E-08, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O 1.1283e-10, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O - 1.8137e-12 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O + 1.8137e-12 // Al2Si2O5(OH)4(s) + 6H+ ⇌ 2Al3+ + 2 SiO2 + 5H2O }; @@ -136,7 +136,7 @@ constexpr CArrayWrapper< double, 24 > reverseRateConstant = 4.5420e-21, //4.5584E-27, // KAlMg3Si3O10(OH)2(s) + 10H+ ⇌ Al3+ + K+ + 3Mg2+ + 3SiO2(aq) + 6H2O 8.2341e-05, //8.2785E-14, // CaAl2(SiO4)2(s) + 8H+ ⇌ Ca2+ + 2 Al3+ + 2 SiO2(aq) + 4H2O 6.3975e-10, //6.3858E-07, // KAlSi3O8(s) + 4H+ ⇌ Al3+ + K+ + 3 SiO2 + 2H2O - 57.3348 //5.7645E-08 // Al2Si2O5(OH)4(s) + 6H+ ⇌ Al3+ + 2 SiO2 + 5H2O + 0.0574 //5.7645E-08 // Al2Si2O5(OH)4(s) + 6H+ ⇌ 2Al3+ + 2 SiO2 + 5H2O }; constexpr CArrayWrapper< int, 24 > mobileSpeciesFlag =