Implement DAE solver

[K20shores]

DAE solving is needed for cloud chemistry

Acceptance criteria

• Can create a 4stageDAE parameter solver when constraints are supplied
• Can create a 6stageDAE parameter solver when constraints are supplied
• A unit test with an exemplary example that can be analytically calculated exists
• An integration test for a carbonate system exists
• GPU and CPU options are available
• Configuring a system that needs constraints without DAE parameters throws an error

Ideas

• If needed, add a new temporary variable for matrix inversions is included to calculate the constraints
  • Place them at the top of rosenbrock
• Naming convention, current system size of $N$ concentrations
  • $f(Y)$, forcing function
  • $Y_N$, our jacobian is $N \times N$
  • $M$ constrained parameters (things like $\mathrm{H^+}$ concentrations)
  • Also $M$ constraints
  • This means we solve a constraint function, $G(Y_N, Z_M) = 0$
    • equation 4.4 in Harier et.al needs to compute
      • $\frac{\partial f}{\partial Y}$ ($N \times N$), $\frac{\partial f}{\partial Z}$ ($N \times M$)
      • $\frac{\partial G}{\partial Y}$ ($N \times N$), $\frac{\partial G}{\partial Z}$ ($N \times M$)
      • Then form a block matrix with (skipping other paramters like $\gamma$ and $I$
        • upper left $I - \gamma h \frac{\partial f}{\partial Y}$
        • upper right $- \gamma h \frac{\partial f}{\partial Z}$
        • lower left $- \gamma h\frac{\partial G}{\partial Y}$
        • lower right $- \gamma h\frac{\partial G}{\partial Z}$
      • This new matrix must be inverted with LU decomposition
• Update the state to appropriately size the jacobian_ matrix when constraints are applied
  • The jacobian is set here. Figure out how to make this be the correct size
• Create a new process set, similar to the Cuda process set,
  • Override the vectorized and non vectorized SubtractJacobianTerms to calculate the constraints
• Configure solver builder to use the new constrained process set
• In the rosenbrock implementation update AlphaMinusJacobian
  • Right now, it seems to compute $I-h\gamma$
    • Add logic to limit this calculation only to the upper left elements
    • Add logic to compute the upper right, lower left, and lower right parts of the jacobian when we have a constrained solver
• Write tests
  • One analytical test
  • One integration test
    • carbonate example, precompute the expected answer and get as close as we can