SimVascular / svZeroDSolver

A C++ lumped-parameter solver for blood flow and pressure in hemodynamic networks
https://simvascular.github.io/documentation/rom_simulation.html#0d-solver
Other
7 stars 18 forks source link

Simulate an open-loop 0D model until it converges #109

Closed JonathanPham closed 5 months ago

JonathanPham commented 5 months ago

Addresses issue #17.

These updates to give users the option to simulate an open-loop 0D model until it converges.

In the 0D solver input file, users must specify the number of cardiac cycles, $n$, they want to run (this parameter is number_of_cardiac_cycles in the 0D input file). They will also have the option (via parameter, use_cycle_to_cycle_error in the 0D input file) to specify a target cycle-to-cycle error, $\epsilon_{target}$ (parameter, sim_cycle_to_cycle_percent_error, in the 0D input file).

If only $n$ is specified, then we will just run the 0D simulation for $n$ cycles.

If $\epsilon_{target}$ is also prescribed, then we take separate steps depending on whether or not the 0D model includes a RCR boundary condition.

  1. if the model does not have any RCR boundary conditions: First, we first run the 0D simulation for $n$ cycles. Then, for each cap (inlet or outlet), we compute the cycle-to-cycle error, $\epsilon$, for both pressure and flow rate, where $\epsilon = \text{abs}\left(\frac{QOI{n} - QOI{n-1}}{QOI{n-1}}\right)$, where $QOI$ is a mean quantity of interest. If $\epsilon \leq \epsilon{target}$ for all caps (for both pressure and flow rate), then the model has converged periodically and we can stop. Otherwise, the model has not converged. In this case, we run another cardiac cycle, recompute $\epsilon$, and iterate until convergence.
  2. if the model has only one RCR boundary condition: For the outlet with the RCR boundary condition, we evaluate equation 21 from here to compute an estimated number of cardiac cycles, $n{\infty}$. (Note that $n{\infty}$ approximates that number of cardiac cycles that must be simulated in order to achieve a prescribed periodicity error, $\epsilon{\infty}$. Also, for simplicity, we assume $\epsilon{\infty}$ matches $\epsilon{target}$ in equation 21.) Run the 0D model for $n{\infty}$ cycles, instead of $n$ cycles. Then, for all caps, we compute $\epsilon$ for both pressure and flow rate and print them. No need to run for more cycles if the model has not converged, i.e., $\epsilon > \epsilon_{target}$. Note that this model is required to have only one RCR boundary condition, but it can have other boundary conditions as well.
  3. if the model has multiple RCR boundary conditions: We perform step 2 (above), but for the RCR boundary condition with the largest time constant, $\tau = R_{d}C$.

We also created a new test case for each case above.

  1. if the model does not have any RCR boundary conditions: Created a model based on pulsatileFlow_R_coronary.json, where we output results at all time steps. Then, we simulate the 0D model, compute the mean flow and pressure for the last two cardiac cycles simulated, and confirm that $\epsilon \leq \epsilon_{target}$.
  2. if the model has only one RCR boundary condition AND if the model has multiple RCR boundary conditions: Created a model based on steadyFlow_bifurcationR_R2.json, where we output results at all time steps and the two RCR boundary conditions have different time constants. Also, use a periodic inflow instead of a steady inflow. We compute $n_{\infty}$ per equation 21 of the paper, using the largest RCR time constant. Set number_of_cardiac_cycles to be a value that is different from $n{\infty}$. Then, we simulate the 0D model and confirm that the results have $n{\infty}$ cycles instead of the value specified for number_of_cardiac_cycles.
JonathanPham commented 5 months ago

@menon-karthik or @mrp089, can one of you guys review the code and merge the pull request? Thanks!