[Temp] V2 Review DAE 3 Rosenbrock integration#927
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 #927 +/- ##
==========================================
+ Coverage 94.63% 94.66% +0.03%
==========================================
Files 68 68
Lines 3185 3487 +302
==========================================
+ Hits 3014 3301 +287
- Misses 171 186 +15 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
There was a problem hiding this comment.
Pull request overview
This PR introduces initial plumbing for algebraic constraints (DAE-style) into the CPU solver stack, including a new Constraint/ConstraintSet API, solver-builder wiring, and an integration test that exercises the new builder API.
Changes:
- Add
Constraint(variant wrapper),ConstraintSet, and integrate constraint Jacobian/forcing contributions into Rosenbrock solve setup. - Extend
State/StateParameterswith constraint-related metadata and build Jacobians sized for species + constraint rows. - Add a new integration test target (
equilibrium) and update existing unit tests to use the new constraint container types.
Reviewed changes
Copilot reviewed 20 out of 21 changed files in this pull request and generated 8 comments.
Show a summary per file
| File | Description |
|---|---|
include/micm/constraint/constraint.hpp |
Introduces Constraint as a std::variant wrapper (currently for EquilibriumConstraint). |
include/micm/constraint/constraint_set.hpp |
Stores constraints by value; computes residuals/Jacobian contributions for constraint rows. |
include/micm/constraint/equilibrium_constraint.hpp |
Refactors equilibrium constraint into a standalone concrete type (no longer inherits Constraint). |
include/micm/solver/solver_builder.hpp |
Adds SetConstraints(...) API and stores constraints in the builder. |
include/micm/solver/solver_builder.inl |
Builds ConstraintSet, merges constraint Jacobian sparsity, and passes constraints into solver policies. |
include/micm/solver/state.hpp |
Adds number_of_constraints_, constraint_size_, and upper_left_identity_diagonal_ fields. |
include/micm/solver/state.inl |
Builds Jacobian sized state_size_ + constraint_size_; initializes identity-diagonal metadata. |
include/micm/solver/rosenbrock.hpp |
Threads ConstraintSet through Rosenbrock solver constructors. |
include/micm/solver/rosenbrock.inl |
Adds constraint forcing/Jacobian terms and modifies diagonal regularization behavior for DAE support. |
include/micm/solver/rosenbrock_temporary_variables.hpp |
Allocates Rosenbrock temporaries sized to species + constraint rows. |
include/micm/solver/backward_euler.hpp |
Adds ConstraintSet to constructor signature for API compatibility (not yet used). |
include/micm/solver/solver.hpp |
Removes GetNumberOfSpecies() / GetNumberOfReactions() accessors. |
test/integration/test_equilibrium.cpp |
New integration tests verifying ConstraintSet API and builder acceptance of constraints. |
test/integration/CMakeLists.txt |
Registers the new equilibrium integration test. |
test/unit/constraint/test_constraint_set.cpp |
Updates tests to use std::vector<Constraint> and std::unordered_map. |
test/unit/process/test_process_configuration.cpp |
Comment typo adjustment (still contains one remaining typo). |
test/unit/solver/test_rosenbrock.cpp |
Minor formatting fix. |
test/unit/solver/test_*_policy*.hpp |
Adds <iomanip> and minor formatting fixes. |
Comments suppressed due to low confidence (1)
include/micm/solver/solver.hpp:105
GetNumberOfSpecies()/GetNumberOfReactions()were removed from this public header. If this is part of the public API, this is a breaking change; consider keeping them (possibly delegating tostate_parameters_) or providing a replacement accessor to avoid downstream breakage.
StatePolicy GetState(const std::size_t number_of_grid_cells = 1) const
{
auto state = std::move(StatePolicy(state_parameters_, number_of_grid_cells));
if constexpr (std::is_convertible_v<typename SolverPolicy::ParametersType, RosenbrockSolverParameters>)
{
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| if constexpr (LuDecompositionInPlaceConcept<LuDecompositionPolicy, SparseMatrixPolicy>) | ||
| { | ||
| jacobian_ = | ||
| BuildJacobian<SparseMatrixPolicy>(parameters.nonzero_jacobian_elements_, number_of_grid_cells, state_size_, true); | ||
| BuildJacobian<SparseMatrixPolicy>(parameters.nonzero_jacobian_elements_, number_of_grid_cells, state_size_ + constraint_size_, true); | ||
| auto lu = LuDecompositionPolicy::template GetLUMatrix<SparseMatrixPolicy>(jacobian_, 0, false); | ||
| jacobian_ = std::move(lu); | ||
| } | ||
| else | ||
| { | ||
| jacobian_ = | ||
| BuildJacobian<SparseMatrixPolicy>(parameters.nonzero_jacobian_elements_, number_of_grid_cells, state_size_, false); | ||
| BuildJacobian<SparseMatrixPolicy>(parameters.nonzero_jacobian_elements_, number_of_grid_cells, state_size_ + constraint_size_, false); | ||
| auto lu = LuDecompositionPolicy::template GetLUMatrices<SparseMatrixPolicy, LMatrixPolicy, UMatrixPolicy>( |
There was a problem hiding this comment.
jacobian_ is now built with size state_size_ + constraint_size_, but variables_ is still allocated from parameters.variable_names_.size() (species only). Several solver paths assume these dimensions match (e.g., linear solves over the Jacobian’s dimension), so enabling constraints risks out-of-bounds access/termination unless the state/temporary matrices are made dimension-consistent or solver code is updated to explicitly operate on the extended size.
There was a problem hiding this comment.
@copilot open a new pull request to apply changes based on this feedback
|
@boulderdaze I've opened a new pull request, #928, to work on those changes. Once the pull request is ready, I'll request review from you. |
#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>
|
@boulderdaze I've opened a new pull request, #929, to work on those changes. Once the pull request is ready, I'll request review from you. |
…CAR/micm into v2-review-dae-3-rosenbrock-integration
This draft provides the full diffs between DAE 3 and main. Will be closed without review process