address copilot review

Author: boulderdaze

include/micm/constraint/constraint_error.hpp

@@ -14,6 +14,7 @@ enum class MicmConstraintErrc
   EmptyReactants = MICM_CONSTRAINT_ERROR_CODE_EMPTY_REACTANTS,
   EmptyProducts = MICM_CONSTRAINT_ERROR_CODE_EMPTY_PRODUCTS,
   InvalidStoichiometry = MICM_CONSTRAINT_ERROR_CODE_INVALID_STOICHIOMETRY,
+  DuplicateAlgebraicSpecies = MICM_CONSTRAINT_ERROR_CODE_DUPLICATE_ALGEBRAIC_SPECIES,
 };
 
 namespace std
@@ -41,6 +42,7 @@ class MicmConstraintErrorCategory : public std::error_category
       case MicmConstraintErrc::EmptyReactants: return "Equilibrium constraint requires at least one reactant";
       case MicmConstraintErrc::EmptyProducts: return "Equilibrium constraint requires at least one product";
       case MicmConstraintErrc::InvalidStoichiometry: return "Stoichiometric coefficients must be positive";
+      case MicmConstraintErrc::DuplicateAlgebraicSpecies: return "Multiple constraints map to the same algebraic species";
       default: return "Unknown error";
     }
   }

include/micm/constraint/constraint_set.hpp

@@ -85,7 +85,8 @@ namespace micm
 
         if (!algebraic_variable_ids_.insert(info.row_index_).second)
         {
-          throw std::runtime_error(
+          throw std::system_error(
+              make_error_code(MicmConstraintErrc::DuplicateAlgebraicSpecies),
               "Multiple constraints map to the same algebraic species row '" + algebraic_species + "'");
         }
 

include/micm/constraint/types/equilibrium_constraint.hpp

@@ -183,7 +183,7 @@ namespace micm
 
             state.ForEachRow([stoich](const double& conc, double& product)
             {
-              product *= std::pow(conc, stoich);
+              product *= std::pow(std::max(0.0, conc), stoich);
             }, state.GetConstColumnView(species_idx), product_product);
           }
 

include/micm/constraint/types/linear_constraint.hpp

@@ -3,6 +3,7 @@
 #pragma once
 
 #include <micm/constraint/constraint_error.hpp>
+#include <micm/constraint/constraint_info.hpp>
 #include <micm/system/stoich_species.hpp>
 
 #include <cstddef>
@@ -99,6 +100,12 @@ namespace micm
           // Create a variable for accumulating the linear sum
           auto linear_sum = force.GetRowVariable();
 
+          // Initialize linear_sum to 0.0 before accumulation
+          state.ForEachRow([](double& sum)
+          {
+            sum = 0.0;
+          }, linear_sum);
+
           for (std::size_t i = 0; i < coeffs.size(); ++i)
           {
             const double coeff = coeffs[i];

include/micm/util/error.hpp

@@ -71,3 +71,4 @@
 #define MICM_CONSTRAINT_ERROR_CODE_EMPTY_REACTANTS               3
 #define MICM_CONSTRAINT_ERROR_CODE_EMPTY_PRODUCTS                4
 #define MICM_CONSTRAINT_ERROR_CODE_INVALID_STOICHIOMETRY         5
+#define MICM_CONSTRAINT_ERROR_CODE_DUPLICATE_ALGEBRAIC_SPECIES   6

test/unit/constraint/test_constraint_set_policy.hpp

@@ -174,7 +174,7 @@ void testAddForcingTerms()
   // State with 2 grid cells
   DenseMatrixPolicy state(2, num_species);
   state[0] = { 0.2, 0.4, 0.6 };  // Away from equilibrium
-  state[1] = { 0.1, 0.3, 0.7 };  // At equilibrium
+  state[1] = { 0.1, 0.3, 0.7 };  // Away from equilibrium
 
   // Forcing vector (same size as state; constraint replaces AB row at index 2)
   DenseMatrixPolicy forcing(2, num_species, 0.0);