Simulate an open-loop 0D model until it converges

[JonathanPham]

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]

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

[menon-karthik]

This looks really great @JonathanPham, thanks for doing this!

I have one minor comment (apart from the one about the print message) - would it be possible to for you to check what the computational cost of checking the convergence is? Maybe take the test cases you added, run them ~100 times with the convergence check on and ~100 times with the convergence check off (but the same number of cardiac cycles). Just wondering what the impact of this check is on the performance, and if it is measurable we can add a note to the documentation saying so. What do you think?

[menon-karthik]

I think we should print this message even in non-debug builds so the user has this information when it occurs. Similar to the messages below showing the percentage error in flow/pressure.

[JonathanPham]

That's a good idea. I can check on the computational cost and report back. I can also print the above message as well.

[JonathanPham]

In the following tests, I set use_cycle_to_cycle_error to be either true or false. In both cases, number_of_cardiac_cycles is set to be the same value. The times reported below correspond to the total time needed to run each case 100 times.

Note that for both tests, I set number_of_cardiac_cycles to be the total number of cardiac cycles found when use_cycle_to_cycle_error is true.

1. Test pulsatileFlow_R_coronary_cycle_error.json: true takes 12.9 seconds and false takes 1.7 seconds --> 7.5x slowdown when checking cycle-to-cycle errors
2. Test pulsatileFlow_bifurcationR_RCR_cycle_error.json: true takes 0.093 seconds and false takes 0.083 seconds --> 1.12x slowdown to compute new number of cardiac cycles

The above timing results suggest that with the current code, when we check for cycle-to-cycle errors, the computational cost increases almost an order of magnitude for models without RCR boundary conditions. This is quite expensive. However, I have an idea of how to update the code to make this faster. Will work on adding these updates next.

Script used for timings:
test_times.txt

[JonathanPham]

Just reran the above timing tests with the newest commit 19d7595 and the timings are pretty much equivalent now (meaning there is no significant slowdown from checking the cycle-to-cycle errors).

@menon-karthik I think we are good to merge now.

[menon-karthik]

Great, glad you checked this!