Merge branch 'master' into Clean-up-code_59

SimVascular · Nov 7, 2023 · 7182a1f · 7182a1f
2 parents 2087819 + a75882c
commit 7182a1f
Show file tree

Hide file tree

Showing 50 changed files with 476 additions and 252 deletions.
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
@@ -12,7 +12,7 @@ jobs:
         run: mkdir docs/build
       - name: Build doxygen documentation
         continue-on-error: true
-        uses: mattnotmitt/doxygen-action@v1.9.5
+        uses: mattnotmitt/doxygen-action@edge
         with:
           working-directory: '.'
           doxyfile-path: 'docs/Doxyfile'

diff --git a/.gitignore b/.gitignore
@@ -2,6 +2,9 @@
 __pycache__
 *.pyc
 
+# simulation output
+*.csv
+
 # macOS file system
 .DS_Store
 
@@ -25,4 +28,4 @@ Release*/
 Debug*/
 externals/
 *.so
-build
+build
diff --git a/docs/Doxyfile b/docs/Doxyfile
@@ -189,7 +189,7 @@ ENUM_VALUES_PER_LINE   = 1
 TREEVIEW_WIDTH         = 250
 EXT_LINKS_IN_WINDOW    = NO
 FORMULA_FONTSIZE       = 10
-USE_MATHJAX            = NO
+USE_MATHJAX            = YES
 MATHJAX_FORMAT         = HTML-CSS
 MATHJAX_RELPATH        = http://cdn.mathjax.org/mathjax/latest
 MATHJAX_EXTENSIONS     = 

diff --git a/docs/pages/developer_guide.md b/docs/pages/developer_guide.md
@@ -75,33 +75,53 @@ inspiration, you can look at the existing source files and how they use
 documentation. In the following you can find a short recap of the most important
 commands:
 
-* **Latex equations**: For inline equations use `\f$a+b=c\f$` and for block equations use:
+### Latex equations
+For inline equations use `\f$a+b=c\f$` and for block equations use:
 ```
 \f[
 a+b=c
 \f]
 ```
-* **Citation**: If you want to cite a piece literature in your documentation, add
-    a respective BibTeX citation to `docs/cpp/references.bib` and use `\cite name_of_citation` to
-    cite the document.
-* **Drawing circuits**: As the elements of the svZeroDSolver are often represented
-    in the form of electrical circuits, we make use of [CircuiTikZ](https://ctan.org/pkg/circuitikz?lang=en)
-    to draw circuits in the documentation. The start a CircuitTikZ drawing use
-    the following command:
+
+### Citations
+If you want to cite a piece literature in your documentation, add
+a respective BibTeX citation to `docs/cpp/references.bib` and use `\cite name_of_citation` to
+cite the document.
+
+
+### Drawing circuits
+As the elements of the svZeroDSolver are often represented
+in the form of electrical circuits, we use [CircuiTikZ](https://ctan.org/pkg/circuitikz?lang=en)
+to draw circuits in the documentation (see blocks in Block for examples). 
+To start a CircuitTikZ drawing use the following command:
 ```
 \f[
 \begin{circuitikz}
 ...
 \end{circuitikz}
 \f]
 ```
+We currently use MathJax, which only supports [a couple of LaTeX packages](https://docs.mathjax.org/en/latest/input/tex/extensions/index.html). In our `Doxyfile`, this is set as `USE_MATHJAX = YES`. The equations look nicer than without MathJax. Unfortunately, CircuiTikZ is currently not supported by MathJax. Thus, we store all current schematics already compiled in `docs/png`.
+
+If you are adding new schematics, you can follow these steps:
+1. Set `USE_MATHJAX = NO` in `docs/Doxyfile`
+2. Locally build Doxygen (see next section)
+3. Copy `png` output from `docs/build/latex` to `docs/png`
+4. Include image with `\image html FILENAME_dark.png` in `.h`/`.cpp` Doxygen
+5. Set `USE_MATHJAX = YES` in `docs/Doxyfile`
+6. Locally build Doxygen again
+7. Check if LaTeX and schematic look nice
+8. Commit your changes and new `png`
+
 
-The documentation is automatically build in the GitHub CI/CD and published
+### Build
+The documentation is automatically built in the GitHub CI/CD and published
 on GitHub pages. If you want to build the documentation locally, you can use:
 
 ```
 doxygen docs/Doxyfile
 ```
+You can then view the documentation locally in your browser by opening `docs/build/html/index.html`.
 
 If you do not have Doxygen install you can do that with `brew install doxygen`
 on macOS or with `sudo apt-get install doxygen` on Linux.

diff --git a/docs/pages/main.md b/docs/pages/main.md
@@ -34,12 +34,19 @@ reflect the hemodynamics and physiology of different cardiovascular
 anatomies. These 0D models are governed by differential algebraic equations
 (DAEs).
 
-For more background information on 0D models, have a look at SimVascular's
-[ROM Simulation Guide](http://simvascular.github.io/docsROMSimulation.html).
+For more background information on 0D models, have a look at the detailed documentation:
+- System of equations: SparseSystem
+- Time integration: Integrator
+- Overview of available 0D elements (blocks): Block
 
-* <a href="https://github.com/StanfordCBCL/svZeroDPlus">Source
-repository</a>
-* <a href="https://simvascular.github.io">About SimVascular</a>
+You can find more details about governing equations in individual blocks, for example:
+- BloodVessel
+- BloodVesselJunction
+- WindkesselBC
+
+For implementation details, have a look at the [source code](https://github.com/StanfordCBCL/svZeroDPlus).
+
+[About SimVascular](https://simvascular.github.io)
 
 # Installation
 

diff --git a/docs/png/blood_vessel.png b/docs/png/blood_vessel.png
diff --git a/docs/png/blood_vessel_dark.png b/docs/png/blood_vessel_dark.png
diff --git a/docs/png/blood_vessel_junction.png b/docs/png/blood_vessel_junction.png
diff --git a/docs/png/blood_vessel_junction_dark.png b/docs/png/blood_vessel_junction_dark.png
diff --git a/docs/png/blood_vessel_junction_individual.png b/docs/png/blood_vessel_junction_individual.png
diff --git a/docs/png/blood_vessel_junction_individual_dark.png b/docs/png/blood_vessel_junction_individual_dark.png
diff --git a/docs/png/closed_loop_coronary_b_c.png b/docs/png/closed_loop_coronary_b_c.png
diff --git a/docs/png/closed_loop_coronary_b_c_dark.png b/docs/png/closed_loop_coronary_b_c_dark.png
diff --git a/docs/png/closed_loop_r_c_r_b_c.png b/docs/png/closed_loop_r_c_r_b_c.png
diff --git a/docs/png/closed_loop_r_c_r_b_c_dark.png b/docs/png/closed_loop_r_c_r_b_c_dark.png
diff --git a/docs/png/flow_reference_b_c.png b/docs/png/flow_reference_b_c.png
diff --git a/docs/png/flow_reference_b_c_dark.png b/docs/png/flow_reference_b_c_dark.png
diff --git a/docs/png/junction.png b/docs/png/junction.png
diff --git a/docs/png/junction_dark.png b/docs/png/junction_dark.png
diff --git a/docs/png/open_loop_coronary_b_c.png b/docs/png/open_loop_coronary_b_c.png
diff --git a/docs/png/open_loop_coronary_b_c_dark.png b/docs/png/open_loop_coronary_b_c_dark.png
diff --git a/docs/png/pressure_reference_b_c.png b/docs/png/pressure_reference_b_c.png
diff --git a/docs/png/pressure_reference_b_c_dark.png b/docs/png/pressure_reference_b_c_dark.png
diff --git a/docs/png/resistance_b_c.png b/docs/png/resistance_b_c.png
diff --git a/docs/png/resistance_b_c_dark.png b/docs/png/resistance_b_c_dark.png
diff --git a/docs/png/resistive_junction.png b/docs/png/resistive_junction.png
diff --git a/docs/png/resistive_junction_dark.png b/docs/png/resistive_junction_dark.png
diff --git a/docs/png/windkessel_b_c.png b/docs/png/windkessel_b_c.png
diff --git a/docs/png/windkessel_b_c_dark.png b/docs/png/windkessel_b_c_dark.png
diff --git a/docs/references.bib b/docs/references.bib
@@ -6,11 +6,11 @@ @article{JANSEN2000305
 pages = {305-319},
 year = {2000},
 issn = {0045-7825},
-doi = {https://doi.org/10.1016/S0045-7825(00)00203-6},
-url = {https://www.sciencedirect.com/science/article/pii/S0045782500002036},
+url = {https://doi.org/10.1016/S0045-7825(00)00203-6},
 author = {Kenneth E. Jansen and Christian H. Whiting and Gregory M. Hulbert},
 abstract = {A generalized-α method is developed and analyzed for linear, first-order systems. The method is then extended to the filtered Navier–Stokes equations within the context of a stabilized finite element method. The formulation is studied through the application to laminar flow past a circular cylinder and turbulent flow past a long, transverse groove. The method is formulated to obtain a second-order accurate family of time integrators whose high frequency amplification factor is the sole free parameter. Such an approach allows the replication of midpoint rule (zero damping), Gear's method (maximal damping), or anything in between.}
 }
+
 @article{kim_coronary,
 author = {Kim, H and Vignon-Clementel, Irene and Coogan, Jessica and Figueroa, Carlos and Jansen, Kenneth and Taylor, CA},
 year = {2010},
@@ -19,5 +19,44 @@ @article{kim_coronary
 title = {Patient-Specific Modeling of Blood Flow and Pressure in Human Coronary Arteries},
 volume = {38},
 journal = {Annals of biomedical engineering},
-doi = {10.1007/s10439-010-0083-6}
+doi = {10.1007/s10439-010-0083-6},
+url = {https://doi.org/10.1007/s10439-010-0083-6}
+}
+
+@Article{pfaller22,
+  author           = {Martin R. Pfaller and Jonathan Pham and Aekaansh Verma and Luca Pegolotti and Nathan M. Wilson and David W. Parker and Weiguang Yang and Alison L. Marsden},
+  journal          = {International Journal for Numerical Methods in Biomedical Engineering},
+  title            = {Automated generation of {0D} and {1D} reduced-order models of patient-specific blood flow},
+  year             = {2022},
+  month            = {aug},
+  number           = {10},
+  volume           = {38},
+  doi              = {10.1002/cnm.3639},
+  url              = {https://doi.org/10.1002/cnm.3639},
+  publisher        = {Wiley},
+}
+
+@Article{pfaller21,
+  author           = {Martin R. Pfaller and Jonathan Pham and Nathan M. Wilson and David W. Parker and Alison L. Marsden},
+  journal          = {Annals of Biomedical Engineering},
+  title            = {On the Periodicity of Cardiovascular Fluid Dynamics Simulations},
+  year             = {2021},
+  month            = {jun},
+  doi              = {10.1007/s10439-021-02796-x},
+  url              = {https://doi.org/10.1007/s10439-021-02796-x},
+  publisher        = {Springer Science and Business Media {LLC}},
+}
+
+@Article{vignon04,
+  author    = {Irene E. Vignon and Charles A. Taylor},
+  journal   = {Wave Motion},
+  title     = {Outflow boundary conditions for one-dimensional finite element modeling of blood flow and pressure waves in arteries},
+  year      = {2004},
+  month     = {apr},
+  number    = {4},
+  pages     = {361--374},
+  volume    = {39},
+  doi       = {10.1016/j.wavemoti.2003.12.009},
+  url       = {https://doi.org/10.1016/j.wavemoti.2003.12.009},
+  publisher = {Elsevier {BV}},
 }
diff --git a/src/algebra/Integrator.cpp b/src/algebra/Integrator.cpp
@@ -92,7 +92,13 @@ State Integrator::step(const State& old_state, double time) {
   // Update time-dependent element contributions in system
   model->update_time(system, new_time);
 
+  // Count total number of step calls
+  n_iter++;
+
   for (size_t i = 0; i < max_iter; i++) {
+    // Count total number of nonlinear iterations
+    n_nonlin_iter++;
+
     // Update solution-dependent element contribitions
     model->update_solution(system, y_af, ydot_am);
 
@@ -126,3 +132,7 @@ State Integrator::step(const State& old_state, double time) {
 
   return new_state;
 }
+
+double Integrator::avg_nonlin_iter() {
+  return (double)n_nonlin_iter / (double)n_iter;
+}
diff --git a/src/algebra/Integrator.h b/src/algebra/Integrator.h
@@ -43,22 +43,66 @@
  * @brief Generalized-alpha integrator
  *
  * This class handles the time integration scheme for solving 0D blood
- * flow system.
+ * flow system using the generalized-\f$\alpha\f$ method \cite JANSEN2000305.
  *
- * Flow rate, pressure, and other hemodynamic quantities in 0D models of
- * vascular anatomies are governed by a system of nonlinear
- * differential-algebraic equations (DAEs):
+ * We are interested in solving the DAE system defined in SparseSystem for the
+ * solutions, \f$\mathbf{y}_{n+1}\f$ and \f$\dot{\mathbf{y}}_{n+1}\f$, at the
+ * next time, \f$t_{n+1}\f$, using the known solutions, \f$\mathbf{y}_{n}\f$
+ * and \f$\dot{\mathbf{y}}_{n}\f$, at the current time, \f$t_{n}\f$. Note that
+ * \f$t_{n+1} = t_{n} + \Delta t\f$, where \f$\Delta t\f$ is the time step
+ * size.
  *
+ * Using the generalized-\f$\alpha\f$ method, we launch a predictor step
+ * and a series of multi-corrector steps to solve for \f$\mathbf{y}_{n+1}\f$
+ * and \f$\dot{\mathbf{y}}_{n+1}\f$. Similar to other predictor-corrector
+ * schemes, we evaluate the solutions at intermediate times between \f$t_{n}\f$
+ * and \f$t_{n + 1}\f$. However, in the generalized-\f$\alpha\f$ method, we
+ * evaluate \f$\mathbf{y}\f$ and \f$\dot{\mathbf{y}}\f$ at different
+ * intermediate times. Specifically, we evaluate \f$\mathbf{y}\f$ at
+ * \f$t_{n+\alpha_{f}}\f$ and \f$\dot{\mathbf{y}}\f$ at
+ * \f$t_{n+\alpha_{m}}\f$, where \f$t_{n+\alpha_{f}} = t_{n} +
+ * \alpha_{f}\Delta t\f$ and \f$t_{n+\alpha_{m}} = t_{n} + \alpha_{m}\Delta
+ * t\f$. Here, \f$\alpha_{m}\f$ and \f$\alpha_{f}\f$ are the
+ * generalized-\f$\alpha\f$ parameters, where \f$\alpha_{m} = \frac{3 - \rho}{2
+ * + 2\rho}\f$ and \f$\alpha_{f} = \frac{1}{1 + \rho}\f$. In the 0D solver, we
+ * set the spectral radius, \f$\rho\f$, to be \f$0.1\f$. For each time step, the
+ * procedure works as follows.
+ *
+ * 1. \f$\textbf{Predictor step}\f$: First, we make an initial guess for
+ * \f$\mathbf{y}_{n+1}\f$ and \f$\dot{\mathbf{y}}_{n+1}\f$,
  * \f[
- * \mathbf{E}(\mathbf{y}, t) \cdot \dot{\mathbf{y}}+\mathbf{F}(\mathbf{y}, t)
- * \cdot \mathbf{y}+\mathbf{c}(\mathbf{y}, t)=\mathbf{0} \f]
+ * \mathbf{y}_{n+1} = \mathbf{y}_{n},\\
+ * \dot{\mathbf{y}}_{n+1} = \frac{\gamma - 1}{\gamma}\dot{\mathbf{y}}_{n},
+ * \f]
+ * where \f$\gamma = 0.5 + \alpha_{m} - \alpha_{f}\f$.
+ *
+ * 2. \f$\textbf{Initiator step}\f$: Then, we initialize the values of
+ * \f$\dot{\mathbf{y}}_{n+\alpha_{m}}\f$ and
+ * \f$\mathbf{y}_{n+\alpha_{f}}\f$,
+ * \f[\dot{\mathbf{y}}_{n+\alpha_{m}}^{k=0} = \dot{\mathbf{y}}_{n} +
+ * \alpha_{m}\left(\dot{\mathbf{y}}_{n+1} - \dot{\mathbf{y}}_{n}\right),\\
+ * \mathbf{y}_{n+\alpha_{f}}^{k=0} = \mathbf{y}_{n} +
+ * \alpha_{f}\left(\mathbf{y}_{n+1} - \mathbf{y}_{n}\right).\f]
+ *
+ * 3. \f$\textbf{Multi-corrector step}\f$: Then, for \f$k \in \left[0, N_{int}
+ * - 1\right]\f$, we iteratively update our guess of
+ * \f$\dot{\mathbf{y}}_{n+\alpha_{m}}^{k}\f$ and
+ * \f$\mathbf{y}_{n+\alpha_{f}}^{k}\f$. We desire the residual,
+ * \f$\textbf{r}\left(\dot{\mathbf{y}}_{n+\alpha_{m}}^{k + 1},
+ * \mathbf{y}_{n+\alpha_{f}}^{k + 1}, t_{n+\alpha_{f}}\right)\f$, to be
+ * \f$\textbf{0}\f$. We solve this system using Newton's method. For details,see
+ * SparseSystem.
  *
- * Here, \f$y\f$ is the vector of solution quantities and \f$\dot{y}\f$ is the
- * time derivative of \f$y\f$. \f$N\f$ is the total number of equations and the
- * total number of global unknowns. The DAE system is solved implicitly using
- * the generalized-\f$\alpha\f$ method \cite JANSEN2000305.
+ * 4. \f$\textbf{Update step}\f$: Finally, we update \f$\mathbf{y}_{n+1}\f$ and
+ * \f$\dot{\mathbf{y}}_{n+1}\f$ using our final value of
+ * \f$\dot{\mathbf{y}}_{n+\alpha_{m}}\f$ and
+ * \f$\mathbf{y}_{n+\alpha_{f}}\f$. \f[
+ * \mathbf{y}_{n+1} = \mathbf{y}_{n} +
+ * \frac{\mathbf{y}_{n+\alpha_{f}}^{N_{int}} -
+ * \mathbf{y}_{n}}{\alpha_{f}},\\ \dot{\mathbf{y}}_{n+1} =
+ * \dot{\mathbf{y}}_{n} + \frac{\dot{\mathbf{y}}_{n+\alpha_{m}}^{N_{int}} -
+ * \dot{\mathbf{y}}_{n}}{\alpha_{m}} \f]
  *
- * `SparseSystem`)
  */
 class Integrator {
  private:
@@ -81,6 +125,7 @@ class Integrator {
   SparseSystem system;
   Model* model{nullptr};
 
+
  public:
   /**
    * @brief Construct a new Integrator object
@@ -91,7 +136,7 @@ class Integrator {
    * @param atol Absolut tolerance for non-linear iteration termination
    * @param max_iter Maximum number of non-linear iterations
    */
-  Integrator(Model* model, double time_step_size, double rho, double atol,
+  Integrator(Model *model, double time_step_size, double rho, double atol,
              int max_iter);
 
   /**
@@ -127,7 +172,14 @@ class Integrator {
    * @param time Current time
    * @return New state
    */
+
   State step(const State& state, double time);
+
+  /**
+   * @brief Get average number of nonlinear iterations in all step calls
+   *
+   */
+  double avg_nonlin_iter();
 };
 
 #endif  // SVZERODSOLVER_ALGEBRA_INTEGRATOR_HPP_
diff --git a/src/algebra/SparseSystem.cpp b/src/algebra/SparseSystem.cpp
@@ -37,7 +37,8 @@ SparseSystem::SparseSystem() {}
 SparseSystem::SparseSystem(int n) {
   F = Eigen::SparseMatrix<double>(n, n);
   E = Eigen::SparseMatrix<double>(n, n);
-  D = Eigen::SparseMatrix<double>(n, n);
+  dC_dy = Eigen::SparseMatrix<double>(n, n);
+  dC_dydot = Eigen::SparseMatrix<double>(n, n);
   C = Eigen::Matrix<double, Eigen::Dynamic, 1>::Zero(n);
 
   jacobian = Eigen::SparseMatrix<double>(n, n);
@@ -58,6 +59,9 @@ void SparseSystem::reserve(Model *model) {
   F.reserve(num_triplets.F);
   E.reserve(num_triplets.E);
   D.reserve(num_triplets.D);
+  dC_dy.reserve(num_triplets.D);
+  dC_dydot.reserve(num_triplets.D);
+
   model->update_constant(*this);
   model->update_time(*this, 0.0);
 
@@ -71,8 +75,10 @@ void SparseSystem::reserve(Model *model) {
 
   F.makeCompressed();
   E.makeCompressed();
-  D.makeCompressed();
+  dC_dy.makeCompressed();
+  dC_dydot.makeCompressed();
   jacobian.reserve(num_triplets.F + num_triplets.E);  // Just an estimate
+
   update_jacobian(1.0);  // Update it once to have sparsity pattern
   jacobian.makeCompressed();
   solver->analyzePattern(jacobian);  // Let solver analyze pattern
@@ -89,7 +95,7 @@ void SparseSystem::update_residual(
 
 void SparseSystem::update_jacobian(double e_coeff) {
   jacobian.setZero();
-  jacobian += F + D + E * e_coeff;
+  jacobian += F + dC_dy + (E + dC_dydot) * e_coeff;
 }
 
 void SparseSystem::solve() {