-
Notifications
You must be signed in to change notification settings - Fork 7
Add constraint base classes for DAE support #900
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
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>
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #900 +/- ##
==========================================
- Coverage 94.66% 94.36% -0.31%
==========================================
Files 66 70 +4
Lines 3245 3407 +162
==========================================
+ Hits 3072 3215 +143
- Misses 173 192 +19 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
- 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>
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Pull request overview
This pull request introduces the foundational constraint class hierarchy for DAE (Differential-Algebraic Equation) support in MICM. It is the first of three PRs that will enable the Rosenbrock solver to handle algebraic constraints alongside differential equations for atmospheric chemistry modeling.
Changes:
- Introduces abstract
Constraintbase class defining the interface for algebraic constraints withResidual()andJacobian()methods - Implements
EquilibriumConstraintfor chemical equilibrium relationships (K_eq = [products]/[reactants]) - Adds
ConstraintSetto manage collections of constraints with Jacobian support, mirroring the existingProcessSetpattern - Includes comprehensive unit tests (16 tests) covering various constraint configurations and edge cases
- Provides extensive documentation explaining the design, architecture, and integration patterns
Reviewed changes
Copilot reviewed 11 out of 11 changed files in this pull request and generated 12 comments.
Show a summary per file
| File | Description |
|---|---|
include/micm/constraint/constraint.hpp |
Abstract base class defining the constraint interface with virtual methods for residual and Jacobian computation |
include/micm/constraint/equilibrium_constraint.hpp |
Concrete implementation of equilibrium constraints with support for arbitrary stoichiometry and special handling for zero concentrations |
include/micm/constraint/constraint_set.hpp |
Collection manager that maps species names to indices, pre-computes sparse Jacobian structure, and provides batch operations over grid cells |
include/micm/Constraint.hpp |
Umbrella header providing convenient access to all constraint-related classes |
test/unit/constraint/test_constraint.cpp |
Unit tests for EquilibriumConstraint covering equilibrium conditions, Jacobian calculations, and error cases |
test/unit/constraint/test_constraint_set.cpp |
Comprehensive tests for ConstraintSet including multi-constraint systems, sparse Jacobian operations, and coupled constraints |
test/unit/constraint/CMakeLists.txt |
Build configuration for constraint unit tests |
test/unit/CMakeLists.txt |
Integration of constraint test subdirectory into main test suite |
constraint_classes_walkthrough.md |
Detailed educational documentation with architecture diagrams, API walkthrough, and design patterns |
constraintset_mirrors_processset.md |
Side-by-side comparison showing how ConstraintSet follows ProcessSet patterns for solver integration |
indexing_patterns.md |
Technical deep-dive explaining the 3-level indirection pattern used for species concentration access |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| // Convert cell state to vector for constraint evaluation | ||
| std::vector<double> concentrations(state_variables.NumColumns()); | ||
| for (std::size_t j = 0; j < state_variables.NumColumns(); ++j) | ||
| { | ||
| concentrations[j] = cell_state[j]; | ||
| } |
Copilot
AI
Jan 20, 2026
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The method allocates a temporary std::vector<double> concentrations for each grid cell in the inner loop. This could impact performance for large numbers of grid cells. Consider reusing a pre-allocated buffer or passing a view directly if the DenseMatrixPolicy supports it. The documentation mentions this mirrors ProcessSet, so this pattern may be intentional, but it's worth noting for potential future optimization.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree with this. If possible, we should make Residual operate on the dense matrix directly and elide the copy here
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fixed in commit 2079709 - refactored to pointer-based interface. The Constraint interface now takes raw pointers (concentrations, indices, jacobian buffer), and ConstraintSet reuses a single jac_buffer allocated once per SubtractJacobianTerms call. No per-cell allocations.
| EquilibriumConstraint( | ||
| const std::string& name, | ||
| const std::vector<std::pair<std::string, double>>& reactants, | ||
| const std::vector<std::pair<std::string, double>>& products, | ||
| double equilibrium_constant) | ||
| : Constraint(name), | ||
| reactants_(reactants), | ||
| products_(products), | ||
| equilibrium_constant_(equilibrium_constant) | ||
| { | ||
| if (equilibrium_constant_ <= 0) | ||
| { | ||
| throw std::invalid_argument("Equilibrium constant must be positive"); | ||
| } | ||
|
|
||
| // Build species dependencies list (reactants first, then products) | ||
| std::size_t idx = 0; | ||
| for (const auto& r : reactants_) | ||
| { | ||
| species_dependencies_.push_back(r.first); | ||
| reactant_dependency_indices_.push_back(idx++); | ||
| } | ||
| for (const auto& p : products_) | ||
| { | ||
| species_dependencies_.push_back(p.first); | ||
| product_dependency_indices_.push_back(idx++); | ||
| } | ||
| } |
Copilot
AI
Jan 20, 2026
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The constructor validates that equilibrium_constant is positive but does not validate that reactants and products are non-empty. An equilibrium constraint with no reactants or no products is physically meaningless and would lead to incorrect calculations (e.g., residual would always be -prod([products]) if reactants is empty). Consider adding validation to throw an exception if either vector is empty.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fixed in commit 070a19d - constructor now validates that reactants and products are non-empty, throwing EmptyReactants or EmptyProducts error.
| EquilibriumConstraint( | ||
| const std::string& name, | ||
| const std::vector<std::pair<std::string, double>>& reactants, | ||
| const std::vector<std::pair<std::string, double>>& products, | ||
| double equilibrium_constant) | ||
| : Constraint(name), | ||
| reactants_(reactants), | ||
| products_(products), | ||
| equilibrium_constant_(equilibrium_constant) | ||
| { | ||
| if (equilibrium_constant_ <= 0) | ||
| { | ||
| throw std::invalid_argument("Equilibrium constant must be positive"); | ||
| } | ||
|
|
||
| // Build species dependencies list (reactants first, then products) | ||
| std::size_t idx = 0; | ||
| for (const auto& r : reactants_) | ||
| { | ||
| species_dependencies_.push_back(r.first); | ||
| reactant_dependency_indices_.push_back(idx++); | ||
| } | ||
| for (const auto& p : products_) | ||
| { | ||
| species_dependencies_.push_back(p.first); | ||
| product_dependency_indices_.push_back(idx++); | ||
| } | ||
| } |
Copilot
AI
Jan 20, 2026
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The constructor does not validate that stoichiometric coefficients are positive. Negative or zero stoichiometry values would lead to incorrect equilibrium calculations. Consider adding validation to ensure all stoichiometric coefficients in both reactants and products are positive (> 0).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fixed in commit 070a19d - constructor now validates that all stoichiometric coefficients are positive, throwing InvalidStoichiometry error.
| // Build species dependencies list (reactants first, then products) | ||
| std::size_t idx = 0; | ||
| for (const auto& r : reactants_) | ||
| { | ||
| species_dependencies_.push_back(r.first); | ||
| reactant_dependency_indices_.push_back(idx++); | ||
| } | ||
| for (const auto& p : products_) | ||
| { | ||
| species_dependencies_.push_back(p.first); | ||
| product_dependency_indices_.push_back(idx++); | ||
| } | ||
| } |
Copilot
AI
Jan 20, 2026
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The constructor does not check for duplicate species within reactants, within products, or across both lists. For example, specifying the same species twice in the reactants list (e.g., [{"A", 1.0}, {"A", 2.0}]) would create duplicate entries in species_dependencies_ and lead to incorrect Jacobian calculations. Consider adding validation to detect and either reject duplicates or merge them by summing stoichiometric coefficients.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As K20shores noted, duplicate species are valid (e.g., 2X represented as [X, X]). No change needed.
K20shores
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree with copilot that we should reduce the number of copies if possible
| // Convert cell state to vector for constraint evaluation | ||
| std::vector<double> concentrations(state_variables.NumColumns()); | ||
| for (std::size_t j = 0; j < state_variables.NumColumns(); ++j) | ||
| { | ||
| concentrations[j] = cell_state[j]; | ||
| } |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree with this. If possible, we should make Residual operate on the dense matrix directly and elide the copy here
| // Convert cell state to vector for constraint evaluation | ||
| std::vector<double> concentrations(state_variables.NumColumns()); | ||
| for (std::size_t j = 0; j < state_variables.NumColumns(); ++j) | ||
| { | ||
| concentrations[j] = cell_state[j]; | ||
| } | ||
|
|
||
| auto dep_id = dependency_ids_.begin(); | ||
| auto flat_id = jacobian_flat_ids_.begin(); | ||
|
|
||
| for (const auto& info : constraint_info_) | ||
| { | ||
| // Build indices vector for this constraint | ||
| std::vector<std::size_t> indices(dep_id, dep_id + info.number_of_dependencies_); | ||
|
|
||
| // Compute constraint Jacobian | ||
| std::vector<double> jac = constraints_[info.constraint_index_]->Jacobian(concentrations, indices); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
agreed
| for (const auto& r : reactants_) | ||
| { | ||
| species_dependencies_.push_back(r.first); | ||
| reactant_dependency_indices_.push_back(idx++); | ||
| } | ||
| for (const auto& p : products_) | ||
| { | ||
| species_dependencies_.push_back(p.first); | ||
| product_dependency_indices_.push_back(idx++); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Our lists of products and reactants can have repeated species, as in 2X is represented as [X, X]. I don't think this correctly calculate the stoichiometric coefficient
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh I see, the input requires that to be true. I guess this isn't really an issue
- 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>
|
Addressed review feedback in commits aa71b5d and 89cffe3: @boulderdaze comments:
@copilot comment:
New files:
All 15 constraint tests pass. |
|
@davidfillmore I've opened a new pull request, #901, to work on those changes. Once the pull request is ready, I'll request review from you. |
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>
|
Addressed performance feedback in commit 2079709: @K20shores + @copilot comments (lines 199, 247):
Changes to Constraint interface: // Old (vector-based, caused allocations)
virtual double Residual(const std::vector<double>& concentrations,
const std::vector<std::size_t>& indices) const = 0;
// New (pointer-based, zero allocations in inner loop)
virtual double Residual(const double* concentrations,
const std::size_t* indices) const = 0;Performance impact:
All 54 tests pass. Documentation updated to reflect new interface. |
|
@davidfillmore I've opened a new pull request, #902, to work on those changes. Once the pull request is ready, I'll request review from you. |
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>
boulderdaze
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you for addressing the comments. I have a couple of additional questions, while reviewing the updated PR.
| /// @brief Name of the constraint (for identification/debugging) | ||
| std::string name_; | ||
|
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| /// @brief Name of the constraint (for identification/debugging) | |
| std::string name_; | |
| /// @brief Name of the constraint | |
| std::string name_; | |
| // SPDX-License-Identifier: Apache-2.0 | ||
|
|
||
| #include <micm/constraint/constraint.hpp> | ||
| #include <micm/constraint/equilibrium_constraint.hpp> |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should this test file be named test_equilibrium_constraint.cpp?
constraint.hpp contains the base abstract class.
| // Test at equilibrium: [A] = 0.001, [B] = 0.001, [AB] = 0.001 | ||
| // K_eq * [A] * [B] = 1000 * 0.001 * 0.001 = 0.001 = [AB] | ||
| // Residual should be 0 | ||
| std::vector<double> concentrations = { 0.001, 0.001, 0.001 }; | ||
| std::vector<std::size_t> indices = { 0, 1, 2 }; | ||
|
|
||
| double residual = constraint.Residual(concentrations.data(), indices.data()); | ||
| EXPECT_NEAR(residual, 0.0, 1e-10); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The Residual() function takes pointers to concentration data and indices. The concentration data is stored in the state object. I’m unclear how the specific set of the concentrations used for this constraint is created in the integration test.
Would it make sense for the function to instead take a const reference to the state object, along with indices derived from that state?
| /// @param concentrations Pointer to species concentrations (row of state matrix) | ||
| /// @param indices Pointer to indices mapping species_dependencies_ to concentrations | ||
| /// @return Residual value | ||
| double Residual(const double* concentrations, const std::size_t* indices) const override |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have a question regarding the types of the function parameters. Please refer to the comments in the test files for this class.
The same question applies to Jacobian(...) function.
Summary
Introduces the constraint class hierarchy for algebraic constraints in DAE (Differential-Algebraic Equation) systems:
Constraint: Abstract base class defining the interface for algebraic constraintsEquilibriumConstraint: Implements equilibrium relationships (K_eq = [products]/[reactants])ConstraintSet: Manager for collections of constraints with Jacobian supportThis is PR 1 of 3 in a series that adds DAE constraint solving to MICM.
Documentation for Reviewers
We've prepared detailed documentation to assist with review. Please start here:
constraint_classes_walkthrough.mdconstraintset_mirrors_processset.mdProcessSetandConstraintSetJacobian computation, showing how constraints mirror the existing reaction patternindexing_patterns.mdKey Sections in the Walkthrough
PR Series Overview
dae-1-constraint-classesdae-2-state-infrastructuredae-3-rosenbrock-integrationEach PR builds on the previous and can be reviewed independently. All branches have passing tests.
What's in This PR
New Files
Documentation Files
Changes to Existing Files
test/unit/CMakeLists.txt: Addadd_subdirectory(constraint)Review Focus
Constraintbase classEquilibriumConstraintformulation:G(y) = K_eq * [reactants] - [products] = 0ConstraintSetmanagement of variable indices and Jacobian elementsTest Plan
🤖 Generated with Claude Code