Automatically compute periodicity error

[mrp089]

It would be a great addition if we printed the periodicity error of the last cardiac cycle to give users an idea of how converged the solution is. For every boundary condition, we would need to calculate a time constant τ=RC and estimate the periodic error ε_∞ as in here.

Once we have that, we could also leave the number of cardiac cycles open and let simulations run until reaching a user-defined periodic error ε_∞. I think this would generally increase our accuracy and also avoid unnecessarily many cardiac cycles.

[JonathanPham]

@mrp089 @menon-karthik

1. Do we want to print out periodicity error for only models with RCR boundary conditions at all outlets or for all models (including models with coronary boundary conditions, mixed, etc.)?
2. For models with just resistor boundary conditions, would we just not print anything?

[mrp089]

Happy that you're looking into this, @JonathanPham!

1. Usually, RCR outlets have very similar time constants, so they converge at the same time. However, that's user-dependent, so it could make sense to calculate the maximum time constant and use that to determine a periodic state to be on the safe side.
2. With only resistances, I think two cardiac cycles should be fine.

In addition/instead of printing out an error, we could also define an input parameter with the error metric and the simulation terminates once it reaches that (instead of specifying the number of cardiac cycles).

[mrp089]

If you introduce printouts, be careful to not slow the simulation down. We currently don't have any printouts in the Release build (only in Debug).

[JonathanPham]

I think introducing an input parameter with the periodicity error sounds good.

So, if both the periodicity error and the number of cardiac cycles are provided as inputs, should we terminate the simulation once the error level is achieved, even if the simulation has not run for the specified number of cycles yet?

What if the user wants to run their simulation for a specific number of cycles, instead of running until a certain error is achieved? Or should we just remove this capability?

[mrp089]

I'd say you can't input both, either one or the other, and throw a descriptive error if they're both given

[JonathanPham]

Sounds good!

[JonathanPham]

Just talked to @menon-karthik and we plan to incorporate this feature for all models, including models with non-RCR boundary conditions.

For a given model, we will compute the periodicity error based on the RCR boundary conditions present in the model, but we will print out the error between the last two cardiac cycles at the end of the simulation if the model has non-RCR boundary conditions, since the RCR-based periodicity error does not guarantee the model will converge.

[JonathanPham]

Summary of my thoughts on how to tackle this issue:

In the 0D solver input file, users will be able to specify either

• the number of cardiac cycles they want to run
• or a target periodicity error, $\epsilon_{\infty}$

If both are specified, then we would throw an error and not run the 0D simulation.

If $\epsilon_{\infty}$ is specified, then we would estimate the number of cardiac cycles, $n_{\infty}$, that we need to run for the 0D simulation to achieve $\epsilon_{\infty}$. We would conservatively estimate $n_{\infty}$ by evaluating equation 21 from this paper. Furthermore, $n_{\infty}$ will be based on the RCR boundary condition with the largest time constant, $\tau$, present in the 0D model. Finally, we will print the error between the last two cardiac cycles at the end of the 0D simulation, but we will only print in the Debug build (not in the Release build).