Skip to content

Correct external process algebraic rows jacobian - #1095

Merged
K20shores merged 4 commits into
mainfrom
develop-1094-external-process-algebraic-rows
Oct 8, 2026
Merged

K20shores merged 4 commits into
mainfrom
develop-1094-external-process-algebraic-rows

Conversation

@K20shores

@K20shores K20shores commented Sep 30, 2026 •

Copy link
Copy Markdown
Collaborator

Closes #1094

Before this PR, external model processes added their jacobian terms (SubtractJacobianTerms) to the rows of algebraic variables, and the constraints (SubtractConstraintJacobian) then added their terms on top. However, the row of an algebraic variable in the jacobian should hold only the constraint's Jacobian terms. After this PR, micm still calls SubtractJacobianTerms for all of the rows that a process affects, but it removes the entries in algebraic rows before the constraints add their terms.

In math:

Split the state into two sets of indices, $D$ is the set of differential variables, and $A$ is the set of algebraic variables. Below, $M$ is the mass matrix (defined in solver builder in micm), $y$ is the state variables (concentrations, number concentration, radius, charge, etc), $y'$ is the time derivative of each variable. Without constraints $M=I$, $I$ being the identity matrix, and then $y'$ is just the forcing.

$$F(y)=My', \quad M = \mathrm{diag}(m), \quad m_i = 1 \;\; \forall\, i \in D, \quad m_i = 0 \;\; \forall\, i \in A,$$

We can describe each process that can affect the jacobian or forcing function as one of three things

  • $r(y)$: the built-in reaction rates
  • $p(y)$: external model process rates
  • $g(y)$: the constraint residuals. There is one residual $\forall i \in A$

The correct right-hand sides of the forcing uses the rate only in differential rows and the constraints only in algebraic rows

$$\begin{aligned} F_i(y) &= r_i(y) + p_i(y) \quad &\forall i\in D \\\ F_i(y) &= g_i(y) \quad &\forall i\in A \end{aligned}$$

and then the correct jacobian similarly uses only the constraints in algebraic rows

$$\begin{aligned} J_i(y) &= \partial r_i / \partial y + \partial p_i / \partial y \quad &\forall i\in D \\\ J_i(y) &= \partial g_i / \partial y \quad &\forall i\in A \end{aligned}$$

Before this PR, we were incorrectly computing the jacobian on algebraic rows as

$$J_i(y) = \partial g_i / \partial y + \partial p_i / \partial y \quad \forall i\in A$$

After this PR, we have the corrected jacobian

$$\begin{aligned} J_i(y) &= \partial r_i / \partial y + \partial p_i / \partial y \quad &\forall i\in D \\\ J_i(y) &= \partial g_i / \partial y \quad &\forall i\in A \end{aligned}$$

@K20shores K20shores changed the title Develop 1094 external process algebraic rows Correct external process algebraic rows jacobian Sep 30, 2026
@K20shores
K20shores force-pushed the develop-1094-external-process-algebraic-rows branch 2 times, most recently from d666720 to 04e64da Compare September 30, 2026 20:00
@codecov-commenter

codecov-commenter commented Oct 1, 2026 •

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 97.05882% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 96.14%. Comparing base (f2b7b0a) to head (adcaaaa).

Files with missing lines Patch % Lines
include/micm/solver/external_model_dispatcher.hpp 96.15% 1 Missing ⚠️
Additional details and impacted files
@@           Coverage Diff           @@
##             main    #1095   +/-   ##
=======================================
  Coverage   96.13%   96.14%           
=======================================
  Files          58       58           
  Lines        4791     4825   +34     
=======================================
+ Hits         4606     4639   +33     
- Misses        185      186    +1     

☔ View full report in Codecov by Harness.
📢 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.

K20shores added a commit to NCAR/musica that referenced this pull request Oct 2, 2026
Move the micm pin to the NCAR/micm#1095 commit that reads the
normalized error into a Real before std::pow. With C++23, nvc++ found
std::pow(ScalarView<Real>, Real) ambiguous.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@K20shores
K20shores added this pull request to stack #1099 October 5, 2026 16:08
K20shores and others added 4 commits October 7, 2026 14:06
Make the Clear* methods public, because nvcc does not allow extended
lambdas in private member functions. Make the cloud stub conserve mass.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Remove the cloud chemistry sweep and the fast-reaction test to keep the
PR small.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
@K20shores
K20shores force-pushed the develop-1094-external-process-algebraic-rows branch from 8f7f04a to adcaaaa Compare October 7, 2026 19:14

@dwfncar dwfncar left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All OK

@K20shores
K20shores merged commit 28d84c0 into main Oct 8, 2026
32 of 34 checks passed
@K20shores
K20shores deleted the develop-1094-external-process-algebraic-rows branch October 8, 2026 14:49
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.

External process models write Jacobian terms into algebraic rows

4 participants