DAE constraint restructure to full mass-matrix row replacement#931
DAE constraint restructure to full mass-matrix row replacement#931boulderdaze merged 70 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
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #931 +/- ##
==========================================
+ Coverage 94.04% 94.61% +0.56%
==========================================
Files 68 69 +1
Lines 3797 3973 +176
==========================================
+ Hits 3571 3759 +188
+ Misses 226 214 -12 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
* merge main * fixes missing inline and template args * fix the cuda constructor signature
There was a problem hiding this comment.
Pull request overview
This PR restructures the DAE constraint implementation to use full species-row replacement (mass-matrix diagonal with algebraic rows) and expands unit/regression/integration coverage for constraint-related solver edge cases.
Changes:
- Integrates
ConstraintSetinto solver build + Rosenbrock solve loop (constraint residuals/Jacobian terms; algebraic-row mass-matrix handling; clamp behavior adjusted for DAEs). - Updates
ConstraintSetsemantics to map constraints onto existing species rows (instead of appending extra constraint rows), and updates Jacobian sparsity/flat-id handling accordingly. - Adds new unit/regression/integration tests covering constraint/DAE builds, error normalization behavior, reordering, and solve-time constraint enforcement.
Reviewed changes
Copilot reviewed 29 out of 30 changed files in this pull request and generated 4 comments.
Show a summary per file
| File | Description |
|---|---|
| test/unit/solver/test_solver_builder.cpp | Adds regression tests ensuring DAE builds don’t throw when constrained species appears in kinetics (reactant/product). |
| test/unit/solver/test_rosenbrock.cpp | Adds normalized-error tests ensuring constraint-only columns don’t affect normalization. |
| test/unit/solver/test_lu_decomposition_policy.hpp | Adds <iomanip> include (test infra tweak). |
| test/unit/solver/test_lu_decomposition_in_place_policy.hpp | Adds <iomanip> include (test infra tweak). |
| test/unit/solver/test_linear_solver_policy.hpp | Adds <iomanip> include (test infra tweak). |
| test/unit/solver/test_linear_solver_in_place_policy.hpp | Adds <iomanip> include (test infra tweak). |
| test/unit/process/test_process_configuration.cpp | Fixes a spelling issue in a test comment. |
| test/unit/constraint/test_constraint_set.cpp | Updates tests for row-replacement constraint semantics; adds vectorized-matrix indexing test. |
| test/regression/RosenbrockChapman/regression_test_solve.cpp | Switches a regression to SixStageDifferentialAlgebraicRosenbrockParameters(). |
| test/integration/test_equilibrium.cpp | New integration suite exercising constraint APIs and DAE solve behavior (incl. reorder + multi-cell cases). |
| test/integration/CMakeLists.txt | Registers new equilibrium integration test. |
| include/micm/solver/temporary_variables.hpp | Makes TemporaryVariables polymorphically clonable (Clone()), enabling correct State copy behavior. |
| include/micm/solver/state.inl | Adds mass-matrix diagonal initialization + constraint metadata into State construction. |
| include/micm/solver/state.hpp | Extends StateParameters/State with constraint metadata + mass-matrix diagonal; uses TemporaryVariables::Clone() for deep copy. |
| include/micm/solver/solver_builder.inl | Adds SetConstraints, incorporates constraints into sparsity/Jacobian IDs, sets algebraic-row masks and mass-matrix diagonal. |
| include/micm/solver/solver_builder.hpp | Adds constraints support to builder API and stores constraints in the builder. |
| include/micm/solver/solver.hpp | Adds post-solve clamp logic that avoids clamping algebraic variables; removes species/reaction count accessors. |
| include/micm/solver/rosenbrock_temporary_variables.hpp | Implements Clone() for Rosenbrock temporary variables. |
| include/micm/solver/rosenbrock.inl | Adds constraint residual/Jacobian integration + mass-matrix scaling; adjusts error norm implementation. |
| include/micm/solver/rosenbrock.hpp | Plumbs ConstraintSet through solver constructors; minor template naming fix. |
| include/micm/solver/backward_euler_temporary_variables.hpp | Implements Clone() for Backward Euler temporary variables. |
| include/micm/solver/backward_euler.hpp | Plumbs ConstraintSet through Backward Euler solver constructor for API compatibility. |
| include/micm/process/process_set.hpp | Adds algebraic-row masking: skips kinetic forcing/Jacobian updates for algebraic rows and aligns flat-id vectors. |
| include/micm/cuda/solver/cuda_rosenbrock.hpp | Updates CUDA Rosenbrock constructor signature to accept ConstraintSet. |
| include/micm/constraint/equilibrium_constraint.hpp | Refactors equilibrium constraint to a concrete type used via Constraint variant; adds negative-concentration guarding. |
| include/micm/constraint/constraint_set.hpp | Reworks constraint bookkeeping for row replacement (algebraic row mapping, Jacobian nonzeros/flat-ids, vectorized handling). |
| include/micm/constraint/constraint.hpp | Refactors constraints to a std::variant wrapper (currently EquilibriumConstraint). |
| docs/design/FULL_REVIEW.md | Removes internal design/review artifact doc. |
| docs/design/CODEX_TODO.md | Removes internal TODO artifact doc. |
| .gitignore | Ignores docs/design/ and CLAUDE.md. |
Comments suppressed due to low confidence (1)
include/micm/constraint/equilibrium_constraint.hpp:12
EquilibriumConstraintusesstd::maxbut this header doesn't include<algorithm>. Since this is a public header (pulled in viaconstraint.hpp/solver_builder.hpp), it should be self-contained; add the missing include to avoid compile failures depending on include order.
#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.
* review and rework * add missing file * add unit test
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 33 out of 33 changed files in this pull request and generated 4 comments.
Comments suppressed due to low confidence (1)
include/micm/constraint/equilibrium_constraint.hpp:12
EquilibriumConstraintusesstd::max(...)but this header doesn't include<algorithm>. Relying on transitive includes is fragile and can break compilation on some standard library implementations; please include<algorithm>here explicitly.
#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.
| for (const auto& info : constraint_info_) | ||
| { | ||
| const std::size_t* indices = dependency_ids_.data() + info.dependency_offset_; | ||
| double residual = constraints_[info.constraint_index_].Residual(concentrations.data(), indices); |
There was a problem hiding this comment.
Should we consider moving the indexing logic into Residual? I don't think we'd need the temporary concentrations vector above. That would reduce our memory footprint a little
| { | ||
| // Get pointer to indices for this constraint | ||
| const std::size_t* indices = dependency_ids_.data() + info.dependency_offset_; | ||
| constraints_[info.constraint_index_].Jacobian(concentration_ptr, indices, jac_buffer.data()); |
| /// @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 |
There was a problem hiding this comment.
I get that turning concentrations into one of our matrix types makes this implementation more complicated, but especially in large models when we might be solving a million grid cells in one state, the copy of concentrations could get large
mattldawson
left a comment
There was a problem hiding this comment.
In general, I like it, particularly that it sets up a structure similar to what exists for processes.
A couple general suggestions:
- maybe move the specific constraints (linear, equilibrium) into a sub-folder, similar to how rate constants are in a sub-folder of
process/. - I agree with Kyle about avoiding allocating memory in the functions called by the linear solver. The functions that calculate forcing and Jacobian terms need to be as performant as possible. So, specifically:
- allocate any working memory needed in the Forcing or Jacobian functions as part of the
State(trying to minimize what's needed as much as possible) - avoid using the square-bracket syntax for matrix objects (this is not performant enough to use in the linear solver workflow). You could specialize functions for the various matrix types (as done in
ProcessSet), but this might be a good opportunity to use the newDenseMatrixPolicy::Function()andSparseMatrixPolicy::Function()that were created for the external aerosol/cloud model API. These should be as performant as the specialized linear solver functions inProcessSetbut don't require specialization for the different matrix orderings.
- allocate any working memory needed in the Forcing or Jacobian functions as part of the
I'm happy to do another review after you have a chance to work through the existing comments.
|
Kyle suggested merging this PR and addressing the comments in a separate PR. I'm happy to do that since resolving the review comments requires changes in multiple places. Let me know what you think! @mattldawson I learned from Kyle that @mattldawson suggested merging this PR, so I'll go ahead and merge it. |
Constraintscan either be solver variables (i.e. species), or parameterized species that aren't solver variables.This implementation supports the first case where the constraints are solver variables.
RosenbrocksolverEquilibriumConstraintfor chemical equilibrium relationshipsLinearConstraintfor mass conservation.A + B + C = constant, and Terminator reaction without any constraint (test/integration/test_linear_constraint.cpp)