Spherical ventricle block (#176)

Co-authored-by: Emma Van Epps <[email protected]>
Co-authored-by: ncdorn <[email protected]>

Author: 3 people

.github/workflows/codechecks.yml

@@ -3,13 +3,12 @@ name: Codechecks
 on: [push, pull_request]
 jobs:
   clang-format:
-    runs-on: ubuntu-latest
+    runs-on: macos-latest
     steps:
       - uses: actions/checkout@v3
       - name: Install dependencies
         run: |
-          wget -O - https://apt.llvm.org/llvm-snapshot.gpg.key | sudo apt-key add -
-          sudo apt-get install clang-format-19
+          brew install clang-format
       - name: Run clang-format
         run: |
           mkdir Release

docs/pages/add_block.md

@@ -2,7 +2,7 @@
 
 [TOC]
 
-Below are details on the steps required to implement a new block in svZeroDSolver.
+Below are details on the steps required to implement a new block in svZeroDSolver. To help with defining the block matrices, you can use the [Jacobian Generator](@ref jacobian).
 
 *Note: The best way to implement a new block is to look at examples of existing block classes. See the `ValveTanh` class for an example.*
 
@@ -68,6 +68,8 @@ Below are details on the steps required to implement a new block in svZeroDSolve
 
 ## 3. Set up the governing equations for the block.
 
+Use the [Jacobian Generator](@ref jacobian) to calculate matrix contributions symbolically.
+
 ### State vector
 
 * The local state vector for each block is always arranged as `y = [P_in, Q_in, P_out, Q_out, InternalVariable_1, ..., InternalVariable_N]`.   
@@ -111,7 +113,7 @@ Below are details on the steps required to implement a new block in svZeroDSolve
 
 ## 4. Implement the matrix equations for the block.
 
-* Implement the `update_constant`, `update_time` and `update_solution` functions.
+* Implement the `update_constant`, `update_time` and `update_solution` functions (which can be output by the [Jacobian Generator](@ref jacobian)).
   * All matrix elements that are constant must be specified in `update_constant`.
   * Matrix elements that depend only on time (not the state variables) must be specified in `update_time`.
   * Matrix elements that change with the solution (i.e. depend on the state variables themselves) must be specified in `update_solution`. 

docs/pages/jacobian.md

@@ -0,0 +1,84 @@
+@page jacobian Jacobian Generator for svZeroDSolver
+
+This tool generates C++ code for block implementations in the svZeroDSolver framework using symbolic mathematics.
+
+## Overview
+
+The script `script/jacobian.py` reads block definitions from YAML files and generates C++ code for implementing the mathematical models in the solver. It uses symbolic differentiation (via SymPy) to automatically derive the necessary Jacobian matrices and system contributions.
+
+## Usage
+
+```bash
+python jacobian.py <yaml_file>
+```
+
+## YAML File Format
+
+The YAML files define the mathematical model of a block with the following structure:
+
+### Required Sections
+
+- `variables`: List of variable names in the model
+- `derivatives`: List of derivative names (must match variables with `_dt` suffix)
+- `constants`: List of parameter names
+- `residuals`: List of residual equations that define the system
+
+### Optional Sections
+
+- `time_dependent`: List of parameters that depend on time (e.g., activation functions)
+- `helper_functions`: Python code defining helper functions used in the residuals
+
+### Example
+
+```yaml
+variables:
+  - Pin
+  - Qin
+  - Pout
+  - Qout
+
+derivatives:
+  - dPin_dt
+  - dQin_dt
+  - dPout_dt
+  - dQout_dt
+
+constants:
+  - R
+  - C
+  - L
+  - S
+
+residuals:
+  - Pin - Pout - (R + S * abs(Qin)) * Qin - L * dQout_dt
+  - Qin - Qout - C * dPin_dt + C * (R + 2 * S * abs(Qin)) * dQin_dt
+```
+
+## Output
+
+The script generates three C++ function implementations:
+
+1. `update_constant` - Sets up constant matrix coefficients for the system
+2. `update_time` - Updates time-dependent parameters
+3. `update_solution` - Computes solution-dependent terms and Jacobians
+
+## Workflow
+
+1. Create a YAML file defining your block's mathematical model
+2. Run `jacobian.py` on this file to generate C++ code
+3. Copy the generated code to your block implementation file
+4. Complete the implementation with necessary boilerplate code
+
+## Tips
+
+- Ensure the number of variables equals the number of derivatives
+- The number of residuals should be equal to the number of variables minus 2
+- Use helper functions for complex expressions to improve readability
+- Define time-dependent constants separately
+
+## Examples
+
+See the provided YAML examples:
+- `ChamberSphere.yaml` - Spherical heart chamber model
+- `BloodVessel.yaml` - Blood vessel model with optional stenosis
+- `ClosedLoopCoronaryBC.yaml` - Coronary boundary condition model

docs/references.bib

@@ -111,3 +111,17 @@ @article{kerckhoffs2007coupling
   year={2007},
   publisher={Springer}
 }
+
+@article{caruel13,
+ author = {M. Caruel and R. Chabiniok and P. Moireau and Y. Lecarpentier and D. Chapelle},
+ doi = {10.1007/s10237-013-0544-6},
+ journal = {Biomechanics and Modeling in Mechanobiology},
+ month = {dec},
+ number = {4},
+ pages = {897--914},
+ publisher = {Springer Science and Business Media {LLC}},
+ timestamp = {2022-02-28},
+ title = {Dimensional reductions of a cardiac model for effective validation and calibration},
+ volume = {13},
+ year = {2013}
+}

scripts/BloodVessel.yaml

@@ -0,0 +1,21 @@
+variables:
+  - Pin
+  - Qin
+  - Pout
+  - Qout
+
+derivatives:
+  - dPin_dt
+  - dQin_dt
+  - dPout_dt
+  - dQout_dt
+
+constants:
+  - R
+  - C
+  - L
+  - S
+
+residuals:
+  - Pin - Pout - (R + S * abs(Qin)) * Qin - L * dQout_dt
+  - Qin - Qout - C * dPin_dt + C * (R + 2 * S * abs(Qin)) * dQin_dt

scripts/ChamberSphere.yaml

@@ -0,0 +1,52 @@
+variables:
+  - Pin
+  - Qin
+  - Pout
+  - Qout
+  - radius
+  - velo
+  - stress
+  - tau
+  - volume
+
+derivatives:
+  - dPin_dt
+  - dQin_dt
+  - dPout_dt
+  - dQout_dt
+  - dradius_dt
+  - dvelo_dt
+  - dstress_dt
+  - dtau_dt
+  - dvolume_dt
+
+constants:
+  - rho
+  - thick0
+  - radius0
+  - W1
+  - W2
+  - eta
+  - act
+  - act_plus
+  - sigma_max
+
+time_dependent:
+  - act
+  - act_plus
+
+helper_functions: |
+  def CG(radius):
+      return (1 + (radius / radius0)) ** 2
+
+  def dCG(radius, dradius_dt):
+      return 2 * (1 + (radius / radius0)) * (1 / radius0) * dradius_dt
+
+residuals:
+  - rho * thick0 * dvelo_dt + (thick0 / radius0) * (1 + (radius / radius0)) * stress - Pout * CG(radius)
+  - -stress + tau + 4 * (1 - CG(radius) ** -3) * (W1 + CG(radius) * W2) + 2 * eta * dCG(radius, dradius_dt) * (1 - 2 * CG(radius) ** -6)
+  - 4 * pi * radius0 ** 2 * CG(radius) * velo - dvolume_dt
+  - dtau_dt + act * tau - sigma_max * act_plus
+  - dradius_dt - velo
+  - Qin - Qout - dvolume_dt
+  - Pin - Pout

scripts/ClosedLoopCoronaryBC.yaml

@@ -0,0 +1,32 @@
+variables:
+  - Pin
+  - Qin
+  - Pout
+  - Qout
+  - Vim
+  - Pa
+
+derivatives:
+  - dPin_dt
+  - dQin_dt
+  - dPout_dt
+  - dQout_dt
+  - dVim_dt
+  - dPa_dt
+
+constants:
+  - Ra
+  - Ram
+  - Rv
+  - Ca
+  - Cim
+  - Pim
+
+time_dependent:
+  - Pim
+
+residuals:
+  - Pout - Pin + (Ram + Ra) * Qin + Rv * Qout + Ram * Ca * dPa_dt - Ram * Ca * dPin_dt + Ram * Ra * Ca * dQin_dt
+  - Qin - Qout + Ca * dPa_dt - Ca * dPin_dt + Ca * Ra * dQin_dt - dVim_dt
+  - Cim * Pout + Cim * Rv * Qout - Cim * Pim - Vim
+  - Pa - Pin - Ra * Qin

scripts/README.md

@@ -0,0 +1,3 @@
+# svZeroDSolver scripts
+
+The [Jacobian Generator](../docs/pages/jacobian.md) tool generates C++ code for new block implementations using symbolic math.