SimVascular / svFSIplus

svFSIplus is an open-source, parallel, finite element multi-physics solver.
https://simvascular.github.io/documentation/svfsi.html
Other
8 stars 22 forks source link

Arbitrary function & ODE parser #218

Open zasexton opened 2 months ago

zasexton commented 2 months ago

Use Case

Often for perfusion and metabolic modeling, different systems of differential algebraic equations are included as reaction terms to represent kinetic components of overall metabolism. Often hypotheses about metabolic functionality require the form of these equations to be changed substantially (e.g. the inclusion of fatty acid oxidation in a cardiomyocyte model versus glucose oxidation alone). Therefore, users interested in simulating multiple models for tissue perfusion may require a flexible framework to easily change governing reaction equations.

Problem

Currently, any changes to governing equations require recompiling svFSIplus to integrate. Futhermore, the lack of modularity hinders incorporation of complex reaction terms in more complex advection-diffusion-reaction systems.

Solution

We shall introduce flexible string parser to allow for inclusion of arbitrary algebraic and ordinary differential equations into on-demand c++ functions. Such an approach is attractive for general purpose uses and computational efficiency (as opposed to symbolic solvers which may be large third party packages or too slow for integration in a governing PDE equation).

Alternatives considered

Symbolic approaches may be considered if direct construction of c++ functions cannot easily be accomplished (although this is not desirable).

Additional context

Specifically this will be important for perfusion and metabolic modeling although there are many more general use cases.

Code of Conduct

zasexton commented 1 month ago

Arbitrary ODE Solver Test Cases

This document provides test cases for a new ODE solver using the Van der Pol oscillator and the Rossler attractor. Each test case includes the problem definition, parameters, initial conditions, and expected behavior.

In each case, the only provided inputs are a string value of the ode system, initial conditions, known constant values, and the time span for the initial value problem. No c++ functions have been compiled a priori as these numerical functions (and, if needed, their Jacobians) should be constructed by the arbitrary solver.

Explanation

These test cases are designed to assess initial functionality and performance of the new ODE solver on both periodic and chaotic systems.

Test Case 1: Van der Pol Oscillator

Problem Definition

The Van der Pol oscillator is a non-conservative oscillator with non-linear damping. It is described by the following second-order differential equation:

$$ \frac{d^2x}{dt^2} - \mu (1 - x^2) \frac{dx}{dt} + x = 0 $$

This can be converted into a system of first-order ODEs:

$$ \begin{cases} \frac{dy_0}{dt} = y_1 \ \frac{dy1}{dt} = \mu (1 - y{0}^2) y_1 - y_0 \end{cases} $$

Parameters

Initial Conditions

Time Span

Expected Behavior

The Van der Pol oscillator should exhibit a limit cycle behavior for $(\mu = 1.0)$. The solution should converge to a periodic orbit.

Input Parameters for Solver

Obtained Solution

Combined Results

cpp_van_der_pol_result_combined

$y_0$ Results

cpp_van_der_pol_result_y0

$y_1$ Results

cpp_van_der_pol_result_y1

Test Case 2: Rossler Attractor

The Rossler Attractor is a system of three coupled non-linear differential equations. It is described by:

$$ \begin{cases} \frac{dy_0}{dt} = -y_1 - y_2 \ \frac{dy_1}{dt} = y_0 + ay_1 \ \frac{dy_2}{dt} = b + y_2(y_0 - c) \end{cases} $$

Where:

Parameters

Initial Conditions

Time Span

Expected Behavior

The Rossler Attractor should exhibit chaotic behavior with a fractal structure in the state space. The trajectory should be sensitive to initial conditions and parameters.

Input Parameters for Solver

Obtained Solution

Combined Results

rossler_attractor_result_combined

$y_0$ Results

rossler_attractor_result_y0

$y_1$ Results

rossler_attractor_result_y1

$y_2$ Results

rossler_attractor_result_y2

Summary & Importance

Initial testing demonstrates that strings representing ODE systems can be effectively implemented and yield solutions to some classic problems. Both the solver and the ODE constructor leverage base classes already defined in svFSIplus. This added functionality may help to extend the solver package to be more modular and avoid hard-coding of ODE physics. This added functionality requires one additional third party package (exprtk.h) and, optionally, a symbolic engine library (symengine) which is used in a preprocessing stage. This work is crucial for the subsequent implementation of a systems biology metabolic model.

Next Steps

The obtained solutions for these two test cases will be validated against benchmark ODE solver packages to test the accuracy of solutions. If satisfactory, efficiency testing and code clean-up will be started to speed up numeric evaluations (if needed) and provide performance metrics for solutions.

ktbolt commented 1 month ago

Very nice @zasexton !

Should these tests be added to the svFSIplus Unit Tests?

Note that the svZeroDSolver is also using exprtk.hpp.

All external software should be placed in svFSI/Code/ThirdParty.

zasexton commented 1 month ago

Reviewing memory efficiency and checking for potential memory leaks using valgrind.

==27238== HEAP SUMMARY:
==27238==     in use at exit: 0 bytes in 0 blocks
==27238==   total heap usage: 16,954,566 allocs, 16,954,566 frees, 720,433,544 bytes allocated
==27238==
==27238== All heap blocks were freed -- no leaks are possible
==27238==
==27238== For lists of detected and suppressed errors, rerun with: -s
==27238== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)

No leaks detected on test cases.

zasexton commented 1 month ago

@ktbolt these might make good test cases in the future but I'm still not sure how exactly I will pull this solver into the svFSIplus code structure yet... after I double check the obtained solutions I will look to clean up this code to test inclusion within svFSIplus

zasexton commented 1 month ago

The following ODE Solver types have been added:

These solvers will be directly validated against their analogous performant solvers available in Scipy over a set of test cases beginning with the Van der Pol oscillator and Rossler Attractor. This framework of general ODE Solvers will make the architecture of svFSIplus more modular

zasexton commented 1 month ago

Attached is a comparison between the rk45 scheme used in svFSIplus and the Scipy for the Rossler Attractor ODE test case. We observe that at long time spans the solutions diverge. These solutions are carried out at equal step sizes. There will have to be a bit of investigation here to determine why this difference emerges. Scipyrk45_cpprk45_comparison_rossler_attractor

zasexton commented 1 month ago

Observed differences in the solutions seem to be from the Dormand-Prince-Shampine (DPS) triple for the interpolation/extrapolation obtained during the Runge-Kutta 45 scheme. The values obtained for the implementation of the code are given here: https://www.ams.org/journals/mcom/1986-46-173/S0025-5718-1986-0815836-3/S0025-5718-1986-0815836-3.pdf

The svFSIplus scheme does not seem to employ this type of DPS triple for the specific implementation of RK

zasexton commented 1 month ago

We are reaching better agreement by including the DPS scheme for efficient error estimation. Still need to finish the full verification of the code but this getting better. Scipyrk45_cpprk45_comparison_rossler_attractor_partial_agreement

zasexton commented 1 month ago

We have achieved verified agreement of the RungeKutta 4th-order with 5th-order error estimate between Scipy and svFSIplus for the rossler attractor test case. Scipyrk45_cpprk45_comparison_rossler_attractor_verified_agreement

zasexton commented 1 month ago

Similar verified agreement observed for the remaining Van der Pol Oscillator using the respective RK45 schemes. Scipyrk45_cpprk45_comparison_van_der_pol_verified_agreement