Add algebraic constraint support to CUDA ProcessSet kernels (#944)

* Add constraint base classes for DAE support

Introduces the constraint class hierarchy for algebraic constraints
in DAE (Differential-Algebraic Equation) systems:

- Constraint: Abstract base class defining the interface
- EquilibriumConstraint: Implements equilibrium relationships (K_eq = products/reactants)
- ConstraintSet: Manager for collections of constraints with Jacobian support

Includes comprehensive unit tests for all constraint functionality.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Add state and builder infrastructure for DAE constraints

Prepares the solver infrastructure to support algebraic constraints:

State changes:
- Add number_of_constraints_ to StateParameters
- Add constraint_size_ and upper_left_identity_diagonal_ to State
- Size jacobian for species + constraints

SolverBuilder changes:
- Add SetConstraintCount() and SetConstraintNames() methods
- Update Build() to include constraint variables in state

This is backward compatible - when constraint_count_=0 (default),
behavior is unchanged from before.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Integrate DAE constraints into Rosenbrock solver

Completes the DAE constraint system by wiring constraints into
the Rosenbrock solver loop:

Solver changes:
- Add ConstraintSet member to RosenbrockSolver
- Call AddForcingTerms() and SubtractJacobianTerms() in solve loop
- Add SetConstraints() to SolverBuilder with full Build() integration

Key implementation details:
- AlphaMinusJacobian adds alpha to ALL diagonals (ODE and algebraic)
- This treats algebraic constraints as stiff ODEs (ε*z' = g with ε=hγ)
- Constraint Jacobian elements merged with ODE elements at build time

Integration tests:
- test_equilibrium.cpp demonstrates working DAE equilibrium solving

Documentation:
- ARCHITECTURE.md: Implementation details and sign conventions
- TESTS.md: Test case specifications
- CLAUDE.md: Project context and build instructions
- Custom Claude Code skills for testing and debugging

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Add documentation for constraint classes and fix param comments

- Add educational walkthrough of Constraint, EquilibriumConstraint, and
  ConstraintSet classes (constraint_classes_walkthrough.md)
- Add line-by-line explanation of reactant loop indexing logic
  (reactant_loop_explanation.md)
- Add comparison of ProcessSet vs ConstraintSet Jacobian computation
  (forcing_jacobian_parallel.md)
- Fix incorrect @param comments in Constraint base class
- Remove unused <map> include from constraint.hpp

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Standardize mathematical notation in docs and comments

Use consistent notation across documentation and code:
- y for state variables
- F(y) for forcing (uppercase)
- G(y) for constraints (uppercase)

Add note on DAE notation conventions explaining the y/z distinction
from DAE literature and why MICM uses a unified state vector.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Rename and expand indexing patterns documentation

Rename reactant_loop_explanation.md to indexing_patterns.md to better
reflect its broader scope.

Add detailed comparison of 2-level indirection (ProcessSet) vs 3-level
indirection (Constraint classes), explaining:
- How ProcessSet pre-computes global indices at construction
- How Constraint receives indices as a parameter for decoupling
- The ConstraintSet bridge that connects them
- Performance implications and potential optimization paths
- Design philosophy trade-offs (OO vs data-oriented)

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Add architecture overview and test coverage to walkthrough

Architecture Overview section:
- Data flow diagram showing how constraints integrate with solver
- Component responsibilities for Constraint, EquilibriumConstraint,
  ConstraintSet, and Rosenbrock Solver
- Parallel structure comparison of ProcessSet and ConstraintSet

Test Coverage section:
- Summary table of EquilibriumConstraint tests
- Summary table of ConstraintSet tests
- Key scenarios covered by each test suite
- Notes on what's not yet tested (future work)
- Instructions for running constraint tests

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Add flat array memory layout documentation to walkthrough

New section "Flat Array Memory Layout and Strides" explains:
- How MICM uses flat arrays with stride-based access
- Dense matrix access pattern (rows=cells, cols=species)
- Sparse Jacobian access via AsVector() and FlatBlockSize()
- Pre-computed jacobian_flat_ids_ for efficient sparse access
- Note that vectorized variants exist in ProcessSet but not yet
  in ConstraintSet (planned for future optimization)

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Rename forcing_jacobian_parallel.md for clarity

Rename to constraintset_mirrors_processset.md to better describe
the document's purpose: showing how ConstraintSet mirrors ProcessSet's
design for Jacobian computation.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Address PR review feedback: move docs and use Yield type

- Move documentation files from root to docs/design/
- Remove branch name reference from walkthrough intro
- Replace std::pair<std::string, double> with Yield type for
  reactants_ and products_ in EquilibriumConstraint
- Fix dead code in Jacobian calculation (removed unused assignment)
- Update test files to use Yield(Species(...), coeff) syntax

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Use micm-specific error types for constraint validation

- Add MICM_ERROR_CATEGORY_CONSTRAINT and error codes to error.hpp
- Create constraint_error.hpp with MicmConstraintErrc enum and error category
- Replace std::invalid_argument with std::system_error in EquilibriumConstraint
- Replace std::runtime_error with std::system_error in ConstraintSet
- Update tests to expect std::system_error

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Optimize Constraint interface to avoid per-cell allocations

Address PR feedback from K20shores and Copilot regarding performance:
- K20shores (line 199): "make Residual operate on the dense matrix directly"
- K20shores (line 247): "agreed" with Copilot on SubtractJacobianTerms

Refactor Residual() and Jacobian() to use pointers instead of vectors,
eliminating temporary allocations in the inner grid cell loop:

- Change Constraint interface from vector-based to pointer-based
- Add dependency_offset_ and jacobian_flat_offset_ to ConstraintInfo
  for O(1) access to pre-computed index arrays
- Add max_dependencies_ for reusable Jacobian buffer allocation
- AddForcingTerms: access row data directly via &cell_state[0]
- SubtractJacobianTerms: use single reusable jac_buffer across all
  constraints and grid cells

Performance impact:
- Eliminates per-cell vector copy of concentrations (was N species)
- Eliminates per-constraint vector copy of indices
- Eliminates per-constraint Jacobian vector allocation

Update documentation to reflect the new interface.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Add constructor validation for EquilibriumConstraint

Address Copilot review suggestions:
- Validate reactants is not empty (EmptyReactants error)
- Validate products is not empty (EmptyProducts error)
- Validate all stoichiometric coefficients are positive (InvalidStoichiometry error)

Adds corresponding error codes and test cases.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Cleanup.

* Added detailed Codex reviews and recommendations in TODO.

* Added Codex 5.2 at Extra High reasoning depth review pass.

* merge main

* fix the tests

* constraint check is done out of loop

* remove unused function

* Review dae 1 constraint classes (#914)

remove md files, and update the test names

* fix the typo

* add a missing header

* remove unnecessary setter for constraints

* change the type of unique ptr to vector

* remove unncessary docs

* add deep copy for constraint

* temp

* replace unique ptr to variant

* fix the error

* fix test

* code clean up

* revert the previous change for get solver

* make builder copyable

* fix test

* add comments on test

* revert back

* state size

* code cleanup

* address copilot review

* Fix missing constraint Jacobian terms on Rosenbrock rejected-step path (#928)

* Initial plan

* Add constraint Jacobian terms on rejected-step path

Co-authored-by: boulderdaze <55209567+boulderdaze@users.noreply.github.com>

* Fix comment: should say 'Subtract' not 'Add' for constraint Jacobian

Co-authored-by: boulderdaze <55209567+boulderdaze@users.noreply.github.com>

---------

Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com>
Co-authored-by: boulderdaze <55209567+boulderdaze@users.noreply.github.com>

* address copilot review

* nvhpc error

* add params for constructor for cuda

* fix a bug

* Implement full mass-matrix DAE row replacement

* Update DAE review and move design docs

* Fix review issues and add DAE integration tests

- Guard std::pow() against negative concentrations in EquilibriumConstraint
- Selective Max(0.0) clamping: only ODE variables, not algebraic
- Add self-diagonal guarantee for constraint rows in sparse Jacobian
- Remove unused constructor params from Rosenbrock and BackwardEuler
- Fix Dervied->Derived typo, stray semicolon, stale comments
- Remove std::cout from tests and unused <iostream> include
- Add 7 new integration tests: DAE parameter variants, coupled constraints,
  conservation law, stiff coupling, non-unit stoichiometry, selective clamping

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* Fix State copy slicing, Jacobian duplicate deps, and unused-species check

- Add virtual Clone() to TemporaryVariables base class; override in
  RosenbrockTemporaryVariables and BackwardEulerTemporaryVariables.
  State copy ctor/assignment now use Clone() to preserve derived type.
- Change ConstraintSet::SubtractJacobianTerms replace mode from = to -=
  so duplicate dependency columns accumulate correctly.
- Include constraint species in UnusedSpeciesCheck() so species used
  only by constraints are not falsely flagged as unused.
- Move updated REVIEW_CLAUDE.md to docs/design/ with fixed/open status.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* Add multi-grid-cell DAE integration test

3 grid cells with different initial conditions exercise the per-cell
iteration in AddForcingTerms, SubtractJacobianTerms, AlphaMinusJacobian,
and NormalizedError. Verifies per-cell constraint satisfaction, mass
conservation, non-negativity, and independent evolution.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* Fix Jacobian concentration guards and Solve overload clamping

- Apply std::max(0.0, ...) consistently in EquilibriumConstraint::Jacobian
  for aggregate product terms and special-case fallbacks, matching the
  Residual policy to prevent NaN with negative transient concentrations.
- Extract PostSolveClamp() helper in Solver so both Solve overloads
  apply identical post-solve clamping (ODE-only for DAE systems).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* Remove design docs from tracking and add to .gitignore

Design review and architecture docs are working documents not
intended for the main repository.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* Resolve all LOW review items: return type, dead code, sparsity filter, regression tests

- Fix GetAlgebraicSpecies() to return const std::string& instead of copy
- Remove dead append-rows code paths from ConstraintSet, rosenbrock temps, state.inl
- Simplify ConstraintSet constructor to 2 params (always replace mode)
- Filter kinetic sparsity from algebraic rows in solver builder
- Add regression tests: State copy+solve, overlapping species Jacobian, constraint-only unused species
- Update all unit tests to use replace mode exclusively

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

* Fix DAE algebraic-row Jacobian flat-id pruning regression

* Add DAE reorder-state and algebraic-reactant regression tests

* Track architecture doc and untrack internal review notes

* Untrack architecture doc from PR scope

* Fix the test failure for DAE constraint enforcement PR (#933)

* merge main

* fixes missing inline and template args

* fix the cuda constructor signature

* fixes the syntax error

* Review and rework DAE Constraint enforcement PR (#937)

* review and rework

* add missing file

* add unit test

* address the copilot review

* add support for cuda algebraic constrait

* replace bool to uint8_t to supports .data()

* add a missing header

* update set algebraic functions for cuda

* address code review

* add setting algebraic ids tests in process set

---------

Co-authored-by: David Fillmore <fillmore@ucar.edu>
Co-authored-by: Claude Opus 4.5 <noreply@anthropic.com>
Co-authored-by: Copilot <198982749+Copilot@users.noreply.github.com>

Author: 4 people

include/micm/cuda/process/cuda_process_set.cuh

@@ -21,6 +21,11 @@ namespace micm
     /// to calculated Jacobian terms to the device
     void CopyJacobianParams(ProcessSetParam& hoststruct, ProcessSetParam& devstruct);
 
+    /// @brief Copies algebraic variable data from host to device memory
+    /// @param hoststruct Host-side struct containing pointers to algebraic variable data
+    /// @param devstruct Device-side struct with pre-allocated memory for algebraic variables
+    void CopyAlgebraicVariableParams(ProcessSetParam& hoststruct, ProcessSetParam& devstruct);
+
     /// This is the host function that will call the CUDA kernel
     ///   to calculate the forcing terms
     void AddForcingTermsKernelDriver(

include/micm/cuda/process/cuda_process_set.hpp

@@ -62,6 +62,13 @@ namespace micm
     /// @param matrix
     void SetJacobianFlatIds(const SparseMatrixPolicy& matrix);
 
+    /// @brief Marks species rows that should be treated as algebraic (constraints replace ODE rows).
+    ///        Updates algebraic variable IDs after `ProcessSetParam` construction.
+    ///        If algebraic variable IDs are not set post-construction, then this function may not be
+    ///        necessary.
+    /// @param variable_ids Set of variable ids whose forcing/Jacobian rows should not receive kinetic contributions
+    void SetAlgebraicVariableIds(const std::set<std::size_t>& variable_ids);
+
     void AddForcingTerms(const auto& state, const DenseMatrixPolicy& state_variables, DenseMatrixPolicy& forcing)
     const requires(VectorizableDense<DenseMatrixPolicy>);
 
@@ -91,12 +98,14 @@ namespace micm
     hoststruct.number_of_products_ = this->number_of_products_.data();
     hoststruct.product_ids_ = this->product_ids_.data();
     hoststruct.yields_ = this->yields_.data();
+    hoststruct.is_algebraic_variable_ = this->is_algebraic_variable_.data();
 
     hoststruct.number_of_reactants_size_ = this->number_of_reactants_.size();
     hoststruct.reactant_ids_size_ = this->reactant_ids_.size();
     hoststruct.number_of_products_size_ = this->number_of_products_.size();
     hoststruct.product_ids_size_ = this->product_ids_.size();
     hoststruct.yields_size_ = this->yields_.size();
+    hoststruct.algebraic_variable_size_ = this->is_algebraic_variable_.size();
 
     // Copy the data from host struct to device struct
     this->devstruct_ = micm::cuda::CopyConstData(hoststruct);
@@ -128,12 +137,14 @@ namespace micm
     hoststruct.number_of_products_ = this->number_of_products_.data();
     hoststruct.product_ids_ = this->product_ids_.data();
     hoststruct.yields_ = this->yields_.data();
+    hoststruct.is_algebraic_variable_ = this->is_algebraic_variable_.data();
 
     hoststruct.number_of_reactants_size_ = this->number_of_reactants_.size();
     hoststruct.reactant_ids_size_ = this->reactant_ids_.size();
     hoststruct.number_of_products_size_ = this->number_of_products_.size();
     hoststruct.product_ids_size_ = this->product_ids_.size();
     hoststruct.yields_size_ = this->yields_.size();
+    hoststruct.algebraic_variable_size_ = this->is_algebraic_variable_.size();
 
     // Copy the data from host struct to device struct
     this->devstruct_ = micm::cuda::CopyConstData(hoststruct);
@@ -172,6 +183,22 @@ namespace micm
     micm::cuda::CopyJacobianParams(hoststruct, this->devstruct_);
   }
 
+  template<typename DenseMatrixPolicy, typename SparseMatrixPolicy>
+    requires(CudaMatrix<DenseMatrixPolicy> && CudaMatrix<SparseMatrixPolicy>)
+  inline void CudaProcessSet<DenseMatrixPolicy, SparseMatrixPolicy>::SetAlgebraicVariableIds(const std::set<std::size_t>& variable_ids)
+  {
+    // Update the host-side is_algebraic_variable_ array
+    ProcessSet<DenseMatrixPolicy, SparseMatrixPolicy>::SetAlgebraicVariableIds(variable_ids);
+
+    // Then update the device memory
+    ProcessSetParam hoststruct;
+    hoststruct.is_algebraic_variable_ = this->is_algebraic_variable_.data();
+    hoststruct.algebraic_variable_size_ = this->is_algebraic_variable_.size();
+  
+    // Copy the data from host struct to device struct
+    micm::cuda::CopyAlgebraicVariableParams(hoststruct, this->devstruct_);
+  }
+
   template<typename DenseMatrixPolicy, typename SparseMatrixPolicy>
     requires(CudaMatrix<DenseMatrixPolicy> && CudaMatrix<SparseMatrixPolicy>)
   inline void CudaProcessSet<DenseMatrixPolicy, SparseMatrixPolicy>::AddForcingTerms(

include/micm/cuda/util/cuda_param.hpp

@@ -3,6 +3,7 @@
 #pragma once
 
 #include <cstddef>
+#include <cstdint>
 #include <utility>
 
 // To make the NormalizedError function works properly on GPU,
@@ -33,6 +34,7 @@ struct ProcessSetParam
   std::size_t* jacobian_product_ids_ = nullptr;
   double* jacobian_yields_ = nullptr;
   std::size_t* jacobian_flat_ids_ = nullptr;
+  uint8_t* is_algebraic_variable_ = nullptr;
   std::size_t number_of_reactants_size_;
   std::size_t reactant_ids_size_;
   std::size_t number_of_products_size_;
@@ -43,6 +45,7 @@ struct ProcessSetParam
   std::size_t jacobian_product_ids_size_;
   std::size_t jacobian_yields_size_;
   std::size_t jacobian_flat_ids_size_;
+  std::size_t algebraic_variable_size_;
 };
 
 /// This struct holds the (1) pointer to, and (2) size of

include/micm/process/process_set.hpp

@@ -11,6 +11,7 @@
 
 #include <algorithm>
 #include <cassert>
+#include <cstdint>
 #include <unordered_map>
 #include <vector>
 
@@ -43,7 +44,7 @@ namespace micm
     std::vector<std::size_t> jacobian_product_ids_;
     std::vector<double> jacobian_yields_;
     std::vector<std::size_t> jacobian_flat_ids_;
-    std::vector<bool> is_algebraic_variable_;
+    std::vector<uint8_t> is_algebraic_variable_;  // uint8_t instead of bool for CUDA compatibility
     std::unordered_map<std::string, std::size_t> variable_map_;
     std::vector<ExternalModelProcessSet<DenseMatrixPolicy, SparseMatrixPolicy>> external_models_;
     std::vector<std::function<void(

src/process/process_set.cu

@@ -23,6 +23,7 @@ namespace micm
       const std::size_t* const __restrict__ d_number_of_products = devstruct.number_of_products_;
       const std::size_t* __restrict__ d_product_ids = devstruct.product_ids_;
       const double* __restrict__ d_yields = devstruct.yields_;
+      const uint8_t* __restrict__ d_is_algebraic_variable = devstruct.is_algebraic_variable_;
       const std::size_t number_of_grid_cells = rate_constants_param.number_of_grid_cells_;
       const double* __restrict__ d_rate_constants = rate_constants_param.d_data_;
       const double* __restrict__ d_state_variables = state_variables_param.d_data_;
@@ -52,12 +53,20 @@ namespace micm
           }
           for (int i_react = 0; i_react < number_of_reactants; ++i_react)
           {
-            d_forcing[(d_reactant_ids[i_react] * cuda_matrix_vector_length) + local_tid] -= rate;
+            const std::size_t row_id = d_reactant_ids[i_react];
+            if (!d_is_algebraic_variable[row_id])
+            {
+              d_forcing[(row_id * cuda_matrix_vector_length) + local_tid] -= rate;
+            }
           }
           const std::size_t number_of_products = d_number_of_products[i_rxn];
           for (int i_prod = 0; i_prod < number_of_products; ++i_prod)
           {
-            d_forcing[d_product_ids[i_prod] * cuda_matrix_vector_length + local_tid] += d_yields[i_prod] * rate;
+            const std::size_t row_id = d_product_ids[i_prod];
+            if (!d_is_algebraic_variable[row_id])
+            {
+              d_forcing[row_id * cuda_matrix_vector_length + local_tid] += d_yields[i_prod] * rate;
+            }
           }
           d_reactant_ids += d_number_of_reactants[i_rxn];
           d_product_ids += d_number_of_products[i_rxn];
@@ -79,8 +88,10 @@ namespace micm
       /// Local device variables
       const ProcessInfoParam* const __restrict__ d_jacobian_process_info = devstruct.jacobian_process_info_;
       const std::size_t* __restrict__ d_reactant_ids = devstruct.jacobian_reactant_ids_;
+      const std::size_t* __restrict__ d_product_ids = devstruct.jacobian_product_ids_;
       const double* __restrict__ d_yields = devstruct.jacobian_yields_;
       const std::size_t* __restrict__ d_jacobian_flat_ids = devstruct.jacobian_flat_ids_;
+      const uint8_t* __restrict__ d_is_algebraic_variable = devstruct.is_algebraic_variable_;
       const std::size_t number_of_grid_cells = rate_constants_param.number_of_grid_cells_;
       const std::size_t number_of_process_infos = devstruct.jacobian_process_info_size_;
       const double* __restrict__ d_rate_constants = rate_constants_param.d_data_;
@@ -109,18 +120,32 @@ namespace micm
           {
             d_rate_d_ind *= d_state_variables[(d_reactant_ids[i_react] * cuda_matrix_vector_length) + local_tid];
           }
-          for (int i_dep = 0; i_dep < process_info.number_of_dependent_reactants_ + 1; ++i_dep)
+          for (int i_dep = 0; i_dep < process_info.number_of_dependent_reactants_; ++i_dep)
           {
-            d_jacobian[*d_jacobian_flat_ids + local_tid] += d_rate_d_ind;
+            const std::size_t row_id = d_reactant_ids[i_dep];
+            if (!d_is_algebraic_variable[row_id])
+            {
+              d_jacobian[*d_jacobian_flat_ids + local_tid] += d_rate_d_ind;
+            }
             ++d_jacobian_flat_ids;
           }
+          if (!d_is_algebraic_variable[process_info.independent_id_])
+          {
+            d_jacobian[*d_jacobian_flat_ids + local_tid] += d_rate_d_ind;
+          }
+          ++d_jacobian_flat_ids;
           for (int i_dep = 0; i_dep < process_info.number_of_products_; ++i_dep)
           {
-            std::size_t jacobian_idx = *d_jacobian_flat_ids + local_tid;
-            d_jacobian[jacobian_idx] -= d_yields[i_dep] * d_rate_d_ind;
+            const std::size_t row_id = d_product_ids[i_dep];
+            if (!d_is_algebraic_variable[row_id])
+            {
+              std::size_t jacobian_idx = *d_jacobian_flat_ids + local_tid;
+              d_jacobian[jacobian_idx] -= d_yields[i_dep] * d_rate_d_ind;
+            }
             ++d_jacobian_flat_ids;
           }
           d_reactant_ids += process_info.number_of_dependent_reactants_;
+          d_product_ids += process_info.number_of_products_;
           d_yields += process_info.number_of_products_;
         }  // end of loop over reactions in a grid cell
       }  // end of checking a CUDA thread id
@@ -137,6 +162,7 @@ namespace micm
       std::size_t number_of_products_bytes = sizeof(std::size_t) * hoststruct.number_of_products_size_;
       std::size_t product_ids_bytes = sizeof(std::size_t) * hoststruct.product_ids_size_;
       std::size_t yields_bytes = sizeof(double) * hoststruct.yields_size_;
+      std::size_t algebraic_variable_bytes = sizeof(uint8_t) * hoststruct.algebraic_variable_size_;
 
       /// Create a struct whose members contain the addresses in the device memory.
       ProcessSetParam devstruct;
@@ -151,6 +177,7 @@ namespace micm
           cudaMallocAsync(&(devstruct.number_of_products_), number_of_products_bytes, cuda_stream_id), "cudaMalloc");
       CHECK_CUDA_ERROR(cudaMallocAsync(&(devstruct.product_ids_), product_ids_bytes, cuda_stream_id), "cudaMalloc");
       CHECK_CUDA_ERROR(cudaMallocAsync(&(devstruct.yields_), yields_bytes, cuda_stream_id), "cudaMalloc");
+      CHECK_CUDA_ERROR(cudaMallocAsync(&(devstruct.is_algebraic_variable_), algebraic_variable_bytes, cuda_stream_id), "cudaMalloc");
 
       /// Copy the data from host to device
       CHECK_CUDA_ERROR(
@@ -180,12 +207,16 @@ namespace micm
       CHECK_CUDA_ERROR(
           cudaMemcpyAsync(devstruct.yields_, hoststruct.yields_, yields_bytes, cudaMemcpyHostToDevice, cuda_stream_id),
           "cudaMemcpy");
+      CHECK_CUDA_ERROR(
+          cudaMemcpyAsync(devstruct.is_algebraic_variable_, hoststruct.is_algebraic_variable_, algebraic_variable_bytes, cudaMemcpyHostToDevice, cuda_stream_id),
+          "cudaMemcpy");
 
       devstruct.number_of_reactants_size_ = hoststruct.number_of_reactants_size_;
       devstruct.reactant_ids_size_ = hoststruct.reactant_ids_size_;
       devstruct.number_of_products_size_ = hoststruct.number_of_products_size_;
       devstruct.product_ids_size_ = hoststruct.product_ids_size_;
       devstruct.yields_size_ = hoststruct.yields_size_;
+      devstruct.algebraic_variable_size_ = hoststruct.algebraic_variable_size_;
 
       return devstruct;
     }
@@ -263,6 +294,24 @@ namespace micm
       devstruct.jacobian_flat_ids_size_ = hoststruct.jacobian_flat_ids_size_;
     }
 
+    void CopyAlgebraicVariableParams(ProcessSetParam& hoststruct, ProcessSetParam& devstruct)
+    {
+      std::size_t algebraic_variable_bytes = sizeof(uint8_t) * hoststruct.algebraic_variable_size_;
+
+      auto cuda_stream_id = micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0);
+
+      CHECK_CUDA_ERROR(
+          cudaMemcpyAsync(
+              devstruct.is_algebraic_variable_,
+              hoststruct.is_algebraic_variable_,
+              algebraic_variable_bytes,
+              cudaMemcpyHostToDevice,
+              cuda_stream_id),
+          "cudaMemcpy");
+
+      devstruct.algebraic_variable_size_ = hoststruct.algebraic_variable_size_;
+    }
+
     /// This is the function that will delete the constant data
     ///   members of class "CudaProcessSet" on the device
     void FreeConstData(ProcessSetParam& devstruct)
@@ -279,6 +328,8 @@ namespace micm
         CHECK_CUDA_ERROR(cudaFreeAsync(devstruct.product_ids_, cuda_stream_id), "cudaFree");
       if (devstruct.yields_ != nullptr)
         CHECK_CUDA_ERROR(cudaFreeAsync(devstruct.yields_, cuda_stream_id), "cudaFree");
+      if (devstruct.is_algebraic_variable_ != nullptr)
+        CHECK_CUDA_ERROR(cudaFreeAsync(devstruct.is_algebraic_variable_, cuda_stream_id), "cudaFree");
       if (devstruct.jacobian_process_info_ != nullptr)
         CHECK_CUDA_ERROR(cudaFreeAsync(devstruct.jacobian_process_info_, cuda_stream_id), "cudaFree");
       if (devstruct.jacobian_reactant_ids_ != nullptr)

test/unit/cuda/process/test_cuda_process_set.cpp

@@ -79,3 +79,12 @@ TEST(CudaProcessSet, CudaMatrix)
   testProcessSet<Group739CudaDenseMatrix, Group739CudaSparseMatrix, micm::CudaProcessSet<Group739CudaDenseMatrix, Group739CudaSparseMatrix>>();
   testProcessSet<Group1130CudaDenseMatrix, Group1130CudaSparseMatrix, micm::CudaProcessSet<Group1130CudaDenseMatrix, Group1130CudaSparseMatrix>>();
 }
+
+TEST(CudaProcessSetAlgebraicVariables, CudaMatrix)
+{
+  testAlgebraicMasking<Group1CudaDenseMatrix, Group1CudaSparseMatrix, micm::CudaProcessSet<Group1CudaDenseMatrix, Group1CudaSparseMatrix>>();
+  testAlgebraicMasking<Group3CudaDenseMatrix, Group3CudaSparseMatrix, micm::CudaProcessSet<Group3CudaDenseMatrix, Group3CudaSparseMatrix>>();
+  testAlgebraicMasking<Group32CudaDenseMatrix, Group32CudaSparseMatrix, micm::CudaProcessSet<Group32CudaDenseMatrix, Group32CudaSparseMatrix>>();
+  testAlgebraicMasking<Group43CudaDenseMatrix, Group43CudaSparseMatrix, micm::CudaProcessSet<Group43CudaDenseMatrix, Group43CudaSparseMatrix>>();
+  testAlgebraicMasking<Group512CudaDenseMatrix, Group512CudaSparseMatrix, micm::CudaProcessSet<Group512CudaDenseMatrix, Group512CudaSparseMatrix>>();
+}

test/unit/process/test_process_set.cpp

@@ -41,4 +41,12 @@ TEST(RandomProcessSet, Matrix)
   testRandomSystem<micm::Matrix<double>, SparseMatrixTest, micm::ProcessSet<micm::Matrix<double>, SparseMatrixTest>>(200, 50, 40);
   testRandomSystem<micm::Matrix<double>, SparseMatrixTest, micm::ProcessSet<micm::Matrix<double>, SparseMatrixTest>>(300, 30, 20);
   testRandomSystem<micm::Matrix<double>, SparseMatrixTest, micm::ProcessSet<micm::Matrix<double>, SparseMatrixTest>>(400, 100, 80);
+}
+
+TEST(ProcessSetAlgebraicVariables, CudaMatrix)
+{
+  testAlgebraicMasking<Group1VectorMatrix, Group1SparseVectorMatrix, micm::ProcessSet<Group1VectorMatrix, Group1SparseVectorMatrix>>();
+  testAlgebraicMasking<Group2VectorMatrix, Group2SparseVectorMatrix, micm::ProcessSet<Group2VectorMatrix, Group2SparseVectorMatrix>>();
+  testAlgebraicMasking<Group3VectorMatrix, Group3SparseVectorMatrix, micm::ProcessSet<Group3VectorMatrix, Group3SparseVectorMatrix>>();
+  testAlgebraicMasking<Group4VectorMatrix, Group4SparseVectorMatrix, micm::ProcessSet<Group4VectorMatrix, Group4SparseVectorMatrix>>();
 }