Skip to content

Dae 3 rosenbrock integration#909

Closed
dwfncar wants to merge 22 commits intomainfrom
dae-3-rosenbrock-integration
Closed

Dae 3 rosenbrock integration#909
dwfncar wants to merge 22 commits intomainfrom
dae-3-rosenbrock-integration

Conversation

@dwfncar
Copy link
Collaborator

@dwfncar dwfncar commented Jan 29, 2026

No description provided.

dwfncar and others added 14 commits January 17, 2026 20:18
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>
@codecov-commenter
Copy link

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.79%. Comparing base (4654f26) to head (1f7ba7d).

Additional details and impacted files
@@                      Coverage Diff                       @@
##           dae-2-state-infrastructure     #909      +/-   ##
==============================================================
+ Coverage                       94.20%   94.79%   +0.59%     
==============================================================
  Files                              70       69       -1     
  Lines                            3432     3440       +8     
==============================================================
+ Hits                             3233     3261      +28     
+ Misses                            199      179      -20     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@boulderdaze boulderdaze changed the base branch from dae-2-state-infrastructure to main January 30, 2026 02:28
@boulderdaze boulderdaze changed the base branch from main to main-before-dev1-constraint February 5, 2026 03:37
@boulderdaze boulderdaze changed the base branch from main-before-dev1-constraint to main February 5, 2026 03:48
@dwfncar
Copy link
Collaborator Author

dwfncar commented Feb 11, 2026

Closing, this PR is superceded by #931

1 similar comment
@dwfncar
Copy link
Collaborator Author

dwfncar commented Feb 11, 2026

Closing, this PR is superceded by #931

@dwfncar dwfncar closed this Feb 11, 2026
@boulderdaze boulderdaze deleted the dae-3-rosenbrock-integration branch March 19, 2026 20:58
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants