Add algebraic constraint support to CUDA ProcessSet kernels#944
Add algebraic constraint support to CUDA ProcessSet kernels#944boulderdaze merged 76 commits intomainfrom
Conversation
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>
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>
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 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>
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 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>
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>
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 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>
- 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>
- 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>
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>
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>
remove md files, and update the test names
…nbrock-integration
* merge main * fixes missing inline and template args * fix the cuda constructor signature
* review and rework * add missing file * add unit test
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #944 +/- ##
=======================================
Coverage 94.61% 94.61%
=======================================
Files 69 69
Lines 3973 3973
=======================================
Hits 3759 3759
Misses 214 214 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
There was a problem hiding this comment.
Pull request overview
Adds end-to-end algebraic constraint (DAE) support with a “replace-species-row” mass-matrix formulation, including honoring algebraic row masks in CUDA ProcessSet kernels (closing the gap described in #938) and expanding test coverage for constraints.
Changes:
- Add constraint plumbing through
SolverBuilder→State(mass-matrix diagonal / algebraic IDs) → solvers (Rosenbrock forcing/Jacobian + normalization + clamping). - Update CPU/CUDA
ProcessSetto skip kinetic forcing/Jacobian writes for algebraic rows (including CUDA device-side mask propagation). - Add/extend unit, integration, and regression tests for equilibrium + linear constraints, and for DAE-specific behaviors (normalization, clamping, state copy).
Reviewed changes
Copilot reviewed 36 out of 36 changed files in this pull request and generated 2 comments.
Show a summary per file
| File | Description |
|---|---|
src/process/process_set.cu |
CUDA kernels now honor algebraic-row masks for forcing/Jacobian updates; device memory now carries is_algebraic_variable_. |
include/micm/process/process_set.hpp |
CPU ProcessSet gains SetAlgebraicVariableIds and skips algebraic rows; Jacobian flat-id generation supports skipped rows. |
include/micm/cuda/util/cuda_param.hpp |
Adds uint8_t* is_algebraic_variable_ to CUDA params for device-side masking. |
include/micm/cuda/process/cuda_process_set.hpp |
Propagates algebraic-row mask to device and overrides SetAlgebraicVariableIds. |
include/micm/cuda/solver/cuda_rosenbrock.hpp |
Updates constructor signature to accept ConstraintSet. |
include/micm/constraint/constraint.hpp |
Reworks constraints into a std::variant wrapper supporting equilibrium + linear constraints. |
include/micm/constraint/equilibrium_constraint.hpp |
Refactors equilibrium constraint to concrete type with dependency bookkeeping + algebraic target selection. |
include/micm/constraint/linear_constraint.hpp |
Introduces LinearConstraint for conservation/linear relations. |
include/micm/constraint/constraint_set.hpp |
Implements “replace species rows” constraint set (residuals + Jacobian contributions) + algebraic variable ID tracking. |
include/micm/solver/solver_builder.hpp |
Adds SetConstraints API and stores constraints on the builder. |
include/micm/solver/solver_builder.inl |
Builds ConstraintSet, marks algebraic rows, prunes kinetic sparsity for those rows, and merges constraint Jacobian sparsity. |
include/micm/solver/state.hpp |
Adds constraint_size_ and upper_left_identity_diagonal_ (mass-matrix diagonal) to state. |
include/micm/solver/state.inl |
Initializes mass-matrix diagonal from parameters (or defaults to ODE-only). |
include/micm/solver/temporary_variables.hpp |
Makes TemporaryVariables polymorphically clonable via pure-virtual Clone(). |
include/micm/solver/rosenbrock_temporary_variables.hpp |
Implements Clone() for Rosenbrock temp vars. |
include/micm/solver/backward_euler_temporary_variables.hpp |
Implements Clone() for Backward Euler temp vars. |
include/micm/solver/rosenbrock.hpp |
Rosenbrock solver now stores a ConstraintSet and accepts it via constructor. |
include/micm/solver/backward_euler.hpp |
Backward Euler solver stores ConstraintSet for API compatibility (not used in algorithm). |
include/micm/solver/rosenbrock.inl |
Integrates constraints into forcing/Jacobian, mass-matrix scaling in linear stage combinations, and error normalization ignoring algebraic vars. |
include/micm/solver/solver.hpp |
Post-solve clamping now avoids clamping algebraic variables; supports constraint-aware state handling. |
test/unit/solver/test_rosenbrock.cpp |
Adds unit tests ensuring NormalizedError ignores algebraic variables. |
test/unit/constraint/test_linear.cpp |
New unit tests for LinearConstraint residual/Jacobian/algebraic target behavior. |
test/unit/constraint/test_constraint_set.cpp |
Updates tests for “replace-row” constraint formulation; adds vectorized sparse ordering coverage. |
test/unit/constraint/CMakeLists.txt |
Adds linear unit test target. |
test/integration/test_equilibrium.cpp |
Adds broad integration coverage for equilibrium constraints across solver configs (DAE solving, reorder state, stiff coupling, etc.). |
test/integration/test_linear_constraint.cpp |
Adds integration test for a linear conservation constraint in a DAE solve. |
test/integration/CMakeLists.txt |
Adds new integration test targets for equilibrium + linear constraints. |
test/regression/RosenbrockChapman/regression_test_solve.cpp |
Adjusts Chapman regression to use six-stage DAE Rosenbrock parameters. |
test/unit/process/test_process_configuration.cpp |
Fixes comment spelling (“stoichiometric coefficients”). |
test/unit/solver/test_lu_decomposition_policy.hpp |
Adds <iomanip> include / minor formatting fix. |
test/unit/solver/test_lu_decomposition_in_place_policy.hpp |
Adds <iomanip> include / minor formatting fix. |
test/unit/solver/test_linear_solver_policy.hpp |
Adds <iomanip> include. |
test/unit/solver/test_linear_solver_in_place_policy.hpp |
Adds <iomanip> include. |
test/integration/analytical_policy.hpp |
Removes a stray standalone ;. |
docs/design/FULL_REVIEW.md |
Deleted design review artifact. |
docs/design/CODEX_TODO.md |
Deleted design TODO artifact. |
Comments suppressed due to low confidence (1)
include/micm/constraint/equilibrium_constraint.hpp:12
EquilibriumConstraintusesstd::maxbut this header doesn’t include<algorithm>. This will fail to compile in translation units that don’t already include<algorithm>before this header. Add#include <algorithm>here (preferred) or avoidstd::maxin the header.
#include <micm/constraint/constraint_error.hpp>
#include <micm/system/stoich_species.hpp>
#include <cmath>
#include <cstddef>
#include <string>
#include <system_error>
#include <vector>
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 5 out of 5 changed files in this pull request and generated 4 comments.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
You can also share your feedback on Copilot code review. Take the survey.
This PR
process_info.number_of_dependent_reactants_ + 1, combining both dependent reactants and the independent variable in one pass. This change helps improve clarity and readibility.