Convert svZeroDSolver from Python to C++

[JonathanPham]

We want to convert the 0D solver from Python to C++. The C++ version will replace the GenBC code and remove the need to perform Python-to-C++ coupling in SimVascular.

I will take a stab at this. I am considering using the Eigen library to perform the matrix computations.

[ktbolt]

@JonathanPham I see that you are writing some C++ code in your cpp branch. Please describe what you are doing in this Issue.

Unit tests, excellent!

I realize you are just experimenting but keep in mind

• Naming convention: PressureVariable is good for class names but use LPNVariable::get_value() rather than LPNVariable::GetValue() for methods or functions.
• Indent two spaces, no tabs
• Except for primitive data types (int, double, etc.) pass values by reference and const them if they are not changed Wire(const std::string name&, const std::array<LPNBlock *, 2>& connecting_block_list)
• Don't use raw pointers unless you really need to because it is easy to lose track of who owns the pointer and they are not safely copied. If you need to use raw pointers then use smart pointers.
• Think carefully about whether to use inheritance or composition, or maybe a template would work.
• Getter and Setter methods are not really needed for simple classes like LPNVariable. You can just make member data public or have name() and value() methods to retrieve member data.

[JonathanPham]

@ktbolt Thanks for the feedback! I finally got around to implementing the changes you recommended.

Regarding the use of smart pointers, I understand that they are very useful when pointers are created using the new keyword. Thus far, in my conversion of blocks.py (to blocks.hpp and blocks.cpp), I have not used the new keyword to create any pointers and I don't have any plans to do so yet. I am planning to declare pointers by using the address-of operator, &, instead. As such, for now, should I continue just using normal pointers (not smart pointers)?

Current work: Currently performing a one-to-one conversion of the Python 0D solver to C++. Will make necessary changes to the code architecture, etc. as needed.

[ktbolt]

@JonathanPham It is fine to use whatever you feel comfortable with. The solver code is not very large and is not an interactive app so leaking memory might not be a problem. The main thing is to implement a code that is clean, documented, understandable, and maintainable (i.e. new functionality to be easily added).

Note that the & operator creates a reference to something. Pointers and references are very different.

[ktbolt]

@richterjakob I just saw that you are doing something with the C++ conversion. Why has this work not been proposed here and no plan discussed?

[mrp089]

@ktbolt, @richterjakob had a divine inspiration last week to rewrite the 0D solver in C++ and just pulled through with it. Even though the core of the 0D solver is already running, we can still discuss it here! @richterjakob can include a link to your branch and briefly outline the steps you have taken?

Some things we have to address before switching over from Python to C++:

• Calling the solver from SimVascular / Python API
• Integrating 0D blocks that use TensorFlow (@natalia-rubio). Maybe we can have a C++ / Python interface for these where all of the block matrices are computed in a Python file?
• Post-processing the results as Python dictionaries saved in NumPy. This could also be done in an external Python file that the solver calls if such an output is desired.

I'm excited that we now have a super-fast version of the 0D solver and by the new possibilities this opens up for us (UQ with >10^6 model evaluations, coupling to 3D solvers, generating synthetic data for machine learning models, performing large-scale parametric studies, ...)!

[ktbolt]

@mrp089 Thanks for the update!

I don't of course wish to dampen inspiration, divine or human, but I do need developers to communicate what they plan to do with the rest of the team before they do it. Other team members may have already spent time thinking about this and may have actual requirements defined.

Increasing solver performance is indeed an important goal but so is implementing a code that is documented, maintainable, and extensible. A hacked C++ version is just a short-term solution that will probably need to be reimplemented in future.

The primary requirements for the 0D solver was to implement it in C++ so it could be loaded as a shared library by another application for BCs (e.g. svSolver and svFSI) to replace GenBC. An application would interact with it in the same process address space using an API that is called just like a function.

The C++ version could have a Python interface, perhaps generated by SWIG.

[richterjakob]

@ktbolt I understand that this might have been a bit uncommunicated. Let me explain, why I made the decision last week. After several discussions about UQ with 0D models, it had become clear that the performance of the 0D solver is crucial for the decision of which methods I can use for my thesis. I estimated whether the python version is sufficient for my needs and unfortunately it was not, or at least it would have limited my UQ possibilities. With my mid-term thesis presentation approaching quickly and therefore only a few days I would be able to spend on this, it was more of a now-or-never decision. Moreover, from what I heard in the lab I was quite sure that nobody is working on this right now.

Now regarding my changes. Although I wrote everything rather quickly, I tried to keep a high quality of code so my changes don’t have to be completely rewritten for future directions. With respect to that, I considered the following points:

• Modularity: the code is completely object-oriented. Similar to the python version it contains classes for blocks and other components, and namespaces for grouping different scopes just to name a few examples.
• Templates: The classes are completely templated to enable easy switching between different scalar types (float, double) or between sparse and dense systems.
• Documentation: The code is completely documented and an HTML documentation is generated via Doxygen. You can see the first version here. On this page, you can also find more information about the architecture of the code. I can also recommend the pages of the element classes (like MODEL::BloodVessel), where I included the circuit drawing for each element.
• Testing: The code is tested by the same tests as the python version (the ones that are based on the currently supported elements) and achieves the same results.
• CI/CD: Automatic building and testing of the C++ version is implemented in a GitHub action (for ubuntu and macOS)

Please have a look at the HTML documentation. It hopefully provides a good overview of everything. I hope that this C++ version meets the requirements you had in mind in the first place or is at least somehow useful to go in that direction. For my thesis, it meant a huge improvement in my possibilities.

Please let me know if you have any comments on the implementation or on how to move forward with this. I would be very happy to hear your opinion on everything.

[richterjakob]

Ah and here is the link to my branch: https://github.com/richterjakob/svZeroDSolver/tree/cpp

[ktbolt]

@richterjakob Well, the code is open source so you can do whatever you want with it when it is convenient for you.

[ktbolt]

@richterjakob I had a look at your C++ code. It is definitely cleaner than the Python code and demonstrates nicely how to use Eigen for the linear algebra but we will not replace the current svZeroDSolver with it.

The header-only implementation does not fit in with the goal of making the solver callable as a plugin from other solvers. I also think that using templates for all of the classes does not make sense. There are other issues with the code and documentation but no need to discuss here.

The code is definitely useable by others so I suggest making a repository for it under your GithHub account.

[mrp089]

I think we should use @richterjakob's (running, tested, and documented) C++ implementation and fix/extend what we need to turn it into a version that can replace the Python one (and can be coupled to 3D solvers).

@ktbolt If you already have extensions/improvements in mind, can you put them in super brief issues? To tag the issues, I've created a milestone.

Thank you for your support!! Even if we got off to an uncoordinated start, I think this is a great chance to build a great 0D solver for us.

[ktbolt]

@mrp089 You are probably right that we should use @richterjakob implementation as the basis for the new C++ code. I would like to clarify the following questions I have (maybe @JonathanPham can help here)

1. Does the code support all of the Python functionality
2. Should the Model::Block class be pure virtual (i.e. just define an interface)
3. Are the circuit equations always nonlinear, seems that simple RCR circuits will be linear
4. Why are we using generalized alpha for the time stepping and not something simpler (e.g. forward Euler) or even explicit

[JonathanPham]

@ktbolt In response to questions 3 and 4:

3. Some circuit equations are linear, like you said. But sometimes, the equations will be nonlinear. For uniformity across the 0D solver, I think it makes sense to create a single solver that can handle both linear and nonlinear equations, rather than developing two separate solvers.
4. Generalized-alpha is a 2nd order implicit time-stepping method. If we want to do something simpler, Implicit Euler can work too (it is only 1st order accurate, however). Explicit solvers (e.g., Forward Euler) have the issue of numerical driftoff when solving algebraic equations using Newton iteration (which is needed by the implicit solver when solving nonlinear equations). Aekaansh looked at this driftoff issue before and it causes the solution to effectively diverge from or be offset from the true solution when solving algebraic systems (e.g., resistor-only system).

[mrp089]

2. Yes

[ktbolt]

@JonathanPham Thanks for reply! It will be interesting to see how this all works for large systems, may want to implement separate linear and nonlinear solvers.

[richterjakob]

@ktbolt

1. Apart from two small details the C++ version has feature parity with the python version. The only features that are missing are: Periodic cubic spline interpolation (linear interpolation is currently used) #70 ; DOF based output writing (only element based output is supported right now)
2. In my opinion yes. The MODEL::Block class defines common members between elements without having an own mechanical behavior.
4. The generalized-alpha acts as a linear solver when the system is linear. It will terminate the non-linear loop after one iteration. Though there was a bug in the python version that resulted in one iteration too many but that is fixed now.

Before writing the C++ version I tried my architecture plans for C++ in the Python version. So there is a Python version with similar architecture than the C++ version, a few bug fixes and more documentation. All changes are on the master branch of my fork. Let me know if we want to keep the Python changes too. The functionality hasn't changed except that a few visualization functions have been removed that were not crucial for the C++ version. If we want to keep them, I will have to re-add them.

[ktbolt]

@menon-karthik needs an interface to the svZeroDSolver from svSolver. I will take @richterjakob C++ code from the https://github.com/richterjakob/svZeroDSolver master branch and start modifying it. I will later merge this into https://github.com/SimVascular/svZeroDSolver.

[ktbolt]

I've added a SolverInterface class to the C++ code that prototypes an interface to the 0D solver. It is compiled into a shared library that can be loaded by a C/C++ program and the interface functions called just like any other function.

The solver is initialized using a JSON file passed to the interface initialization function and returns a problem_id.

    const char* file_name = "steadyFlow_RC_R.json";
    int problem_id;
    initialize(file_name, problem_id);

The solution is incremented in time using the increment_time interface function

  increment_time(problem_id, time);

[ktbolt]

The C++ code does not build on Ubuntu 18 and Ubuntu 20, cmake fails with

Could not find a package configuration file provided by "simdjson" with any
of the following names:

    simdjsonConfig.cmake
    simdjson-config.cmake

The C++ code does build on Ubuntu 22 but none of our SV applications are currently supported for that platform. We can probably build the solvers on Ubuntu 22 OK but I don't think SV will build without new externals.

[richterjakob]

I encountered the same problem when using the package manager dependencies. However, the second way of dependency management via git clone (see cluster build instructions in html documentation) should work for all ubuntu versions.

[ktbolt]

@richterjakob Good idea to just clone the code. I am also going to have a look at using CMakeFetchContent.

[ktbolt]

On MacOS Catalina the 0D solver crashes if I use a MODEL::Model<T> object as a class data member, for example

template <typename T>
class MainTest {
  public:
    MODEL::Model<T> model;
};

MainTest<T> mtest;
mtest.model = config.get_model();

ALGEBRA::Integrator<T, S> integrator(mtest.model, time_step_size, 0.1, absolute_tolerance, max_nliter);

It seems that data is being corrupted in Block objects. This works on Ubuntu 22.

I've changed ConfigReader<T>::get_model() to return MODEL::Model<T>* and things now work, I can run a 0D simulation from an external program.

[ktbolt]

I've done some prototyping from both C++ and Fortran, looks good. People can work off of my branch to interface to 3D solver codes.

There is a lot more work to do to get this code ready to replace the current Python code. I will continue the conversion in a couple of months or so.

[mrp089]

Awesome, @ktbolt, happy to hear!!

Can you list what more work would be necessary from your view to replace the current Python code? Then we already have some steps in place if anybody in the lab is overcome with a spontaneous burst of coding motivation.

[mrp089]

TLDR: C++ development will be moved to a separate repository because @ktbolt does not think the current version is good enough to be in this repository under SimVascular.

@ktbolt, @richterjakob, and I discussed how to proceed with the development of the C++ version. Here's a summary.

Interested parties

• current developers @richterjakob, @ktbolt, @menon-karthik, @natalia-rubio
• potential future developers (e.g. @aabrown100-git and @zinanhu0810 who are interested in 0D-3D coupling)

Background @richterjakob developed a C++ version of svZeroDSolver. It provides the following advantages:

• On average 35 times faster for all N=72 models in this paper (1.8 min to run all models, compared to 61 min before)
• Same (plus more) functionality than the Python implementation
• Tested extensively
• Documented more than the Python version
• Deployed automatically

Not only does a fast solver allow quicker development and many-query applications (UQ, optimization), but the C++ nature also allows coupling to 3D solvers (svSolver, svFSI). Instead of hand-typing massively coupled ODEs in GenBC, we could use the modular svZeroDSolver and just swap out 0D and 3D as desired for fast prototyping.

Development so far C++ (and some Python) development moved between forks and multiple branches and has not been merged into this repository's master:

• @richterjakob base solver: https://github.com/richterjakob/svZeroDSolver
• @ktbolt interface: https://github.com/ktbolt/svZeroDSolver/tree/convert_cpp
• @menon-karthik closed-loop coronary boundary condition: https://github.com/menon-karthik/svZeroDSolver/tree/closed_loop

Future development @ktbolt @ktbolt does not want the merge the current C++ implementation into this repository because it does not meet SimVascular software standards. Thus, C++ development should happen only in people's forks as their private software projects and not be affiliated with SimVascular. He mentioned changes that he wants to implement himself in a couple of months:

• no templates
• more error handling
• implement as a library instead of header-only for 3D solver coupling

He is concerned that developers will be confused if we merge the C++ version now and the implementation changes in a couple of months.

Future development @mrp089 I want to have both C++ and Python versions in this repository now for the following reasons:

• The current C++ version runs 35 times faster with identical results and is well documented
• Anybody can easily use and develop the latest C++ version
• We keep track of development in a central, visible place
• No duplicate code and documentation

Compromise As a compromise, @ktbolt suggested creating a separate, internal repository for all C++ development, splitting the repositories:

• Python https://github.com/SimVascular/svZeroDSolver
• C++ either forked from master or copied from @richterjakob's fork

I will work with @richterjakob and @menon-karthik to merge all development in the new C++ repository.

[mrp089]

I did not want to discount @richterjakob's efforts :wink: If we only export the last time step, we get a speedup factor of 51 in C++ over the Python version (averaged for all N=72 models in the publication).

[mrp089]

Closing this. The C++ implementation has been moved to https://github.com/StanfordCBCL/svZeroDSolver