SimVascular / svZeroDSolver

A C++ lumped-parameter solver for blood flow and pressure in hemodynamic networks
Other
6 stars 18 forks source link

Clean up residual and linearization #53

Closed mrp089 closed 10 months ago

mrp089 commented 11 months ago

While writing up the Methods section of the upcoming 0D-3D-UQ paper with @richterjakob, @menon-karthik, and others, I realized that our definition of the residual and its linearization is a bit more complicated than it needs to be. With $\mathbf\alpha$ being the 0D element parameters, the residual is

$$\textbf{r}\left(\dot{\textbf{y}}, \textbf{y}, \mathbf\alpha\right) = \textbf{E}\left(\textbf{y}, \mathbf\alpha \right) \cdot \dot{\textbf{y}} + \textbf{F}\left( \textbf{y}, \mathbf\alpha \right) \cdot \textbf{y} + \textbf{c}\left(\textbf{y}, \mathbf\alpha \right),$$

and the linearization

$$\textbf{K}\left(\dot{\textbf{y}}, \textbf{y}, \mathbf\alpha \right) = \underset{\text{Term 1}}{\underbrace{\frac{\partial \textbf{E}\left(\textbf{y}, \mathbf\alpha \right)}{\partial \textbf{y}}\cdot\dot{\textbf{y}}}} + \frac{\alpha_{m}}{\Delta t\alpha_{f}\gamma}\textbf{E}\left(\textbf{y}, \mathbf\alpha \right) + \underset{\text{Term 2}}{\underbrace{\frac{\partial \textbf{F}\left(\textbf{y}, \mathbf\alpha \right)}{\partial \textbf{y}}\cdot\textbf{y}}} + \textbf{F}\left(\textbf{y}, \mathbf\alpha \right) + \underset{\text{Term 3}}{\underbrace{\frac{\partial \textbf{c}\left(\textbf{y}, \mathbf\alpha \right)}{\partial \textbf{y}}}}.$$

However, this adds confusion on where to put nonlinear terms ($\textbf{E}$, $\textbf{F}$, or $\textbf{c}$) and also adds three partial derivatives. Therefore, I'm proposing the following residual

$$\textbf{r}\left(\dot{\textbf{y}}, \textbf{y}, \mathbf\alpha\right) = \textbf{E}\left(\mathbf\alpha \right) \cdot \dot{\textbf{y}} + \textbf{F}\left( \mathbf\alpha \right) \cdot \textbf{y} + \textbf{c}\left(\textbf{y}, \mathbf\alpha \right).$$

Here, $\textbf{E}$ and $\textbf{F}$ are constant, and all nonlinear terms are in $\textbf{c}$. The linearzation is also simplified as

$$\mathbf K (\mathbf y, \dot{\mathbf y}, \boldsymbol \alpha) = \frac{\alpha_m}{\alpha_f\gamma\Delta t} \left( \mathbf E + \frac{\partial\mathbf c}{\partial\dot{\mathbf y}} \right) + \mathbf F + \frac{\partial\mathbf c}{\partial\mathbf y}.$$

This requires only the two partial derivatives of the nonlinear term $\textbf{c}$.

I think, in the end, the assembled matrices should look the same in the solver (but this will serve as another sanity check). This will only affect blocks with nonlinear elements, which I think are BloodVessel, BloodVesselJunction, and ClosedLoopHeartPulmonary. I'll make the changes to those blocks and SparseSystem, and update the Doxygen documentation. I'll also move any (non-redundant) content from the ROM Simulation Guide to Doxygen, remove the documentation in simvascular.github.io, and link to Doxygen.