Arbitrary function & ODE parser

[zasexton]

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

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[zasexton]

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

• Van der Pol Oscillator:

  • The system is defined by converting the second-order ODE into two first-order ODEs.
  • Parameters and initial conditions are specified.
  • The expected behavior is a limit cycle.

• Rossler Attractor:

  • The system is defined by three first-order ODEs.
  • Parameters and initial conditions are specified.
  • The expected behavior is chaotic with a fractal structure.

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{dy_1}{dt} = \mu (1 - y_{0}^2) y_1 - y_0 \end{cases} $$

Parameters

• $(\mu = 1.0)$ (standard non-linearity parameter)

Initial Conditions

• $(y_0(0) = 2.0)$
• $(y_1(0) = 0.0)$

Time Span

• $(t \in [0, 20])$

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

• ODE System String: "dy0_dt = y1;
  dy1_dt = mu*(1-y0*y0)*y1-y0"
• Constant Values: {{"mu", 1.0}}
• Initial Conditions: {2.0,0.0}
• Time Span: [0, 20]

Obtained Solution

Combined Results

$y_0$ Results

$y_1$ Results

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:

• $( y_0(t) ), ( y_1(t) ), ( y_2(t) )$ are the state variables.
• $( a ), ( b ), ( c )$ are parameters that define the behavior of the system.

Parameters

• $(a = 0.2)$
• $(b = 0.2)$
• $(c = 5.7)$

Initial Conditions

• $(y_0(0) = 0.0)$
• $(y_1(0) = 1.0)$
• $(y_2(0) = 1.05)$

Time Span

• $(t \in [0, 100])$

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

• ODE System string: "dy0_dt = -y1-y2;\ndy1_dt = y0+ay1;\ndy2_dt = b+y2(y0-c)"
• Constant Values: {{"a", 0.2}, {"b",0.2}, {"c", 5.7}}
• Initial Conditions: {0.0, 1.0, 1.05}
• Time Span: [0,100]

Obtained Solution

Combined Results

$y_0$ Results

$y_1$ Results

$y_2$ Results

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]

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]

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]

@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]

The following ODE Solver types have been added:

• Forward Euler (as implemented in svFSIplus see example:

  svFSIplus/Code/Source/svFSI/CepModAp.cpp

  Lines 182 to 195 in b92add4

  	void CepModAp::integ_fe(const int nX, Vector<double>& X, const double Ts, const double Ti, const double Istim, const double Ksac)
  	{
  	double t = Ts / Tscale;
  	double dt = Ti / Tscale;
  	double Isac = Ksac * (Vrest - X(0));
  	double fext = (Istim + Isac) * Tscale / Vscale;
  	X(0) = (X(0) - Voffset) / Vscale;
  	Vector<double> f(nX);
  	getf(nX, X, f, fext);
  	//CALL AP_GETF(nX, X, f, fext)
  	X = X + dt*f;
  	X(0) = X(0)*Vscale + Voffset;
  	}

  )
• Backward Euler (with fixed-point iteration, no svFSIplus implementation yet)
• Runge Kutta (as implemented in svFSIplus, see example:

  svFSIplus/Code/Source/svFSI/CepModAp.cpp

  Lines 200 to 239 in b92add4

  	void CepModAp::integ_rk(const int nX, Vector<double>& X, const double Ts, const double Ti, const double Istim, const double Ksac)
  	{
  	double t = Ts / Tscale;
  	double dt = Ti / Tscale;
  	double dt6 = dt / 6.0;
  	double Isac = Ksac * (Vrest - X(0));
  	double fext = (Istim + Isac) * Tscale / Vscale;
  	X(0) = (X(0) - Voffset) / Vscale;
  	Array<double> frk(nX,4);
  	Vector<double> Xrk(nX);
  	Vector<double> frk1(nX), frk2(nX), frk3(nX), frk4(nX);
  	// RK4: 1st pass
  	Xrk = X;
  	getf(nX, Xrk, frk1, fext);
  	//CALL AP_GETF(nX, Xrk, frk(:,1), fext)
  	// RK4: 2nd pass
  	Xrk = X + 0.5 * dt * frk1;
  	getf(nX, Xrk, frk2, fext);
  	//CALL AP_GETF(nX, Xrk, frk(:,2), fext)
  	// RK4: 3rd pass
  	Xrk = X + 0.5 * dt * frk2;
  	getf(nX, Xrk, frk3, fext);
  	//CALL AP_GETF(nX, Xrk, frk(:,3), fext)
  	// RK4: 4th pass
  	Xrk = X + dt*frk3;
  	getf(nX, Xrk, frk4, fext);
  	//CALL AP_GETF(nX, Xrk, frk(:,4), fext)
  	X = X + dt6 * (frk1 + 2.0*(frk2 + frk3) + frk4);
  	//X = X + dt6*(frk(:,1) + 2._RKIND*(frk(:,2) + frk(:,3)) + frk(:,4))
  	X(0) = X(0)*Vscale + Voffset;
  	}

  )
• Crank Nicolson (as implemented in svFSIplus, see example:

  svFSIplus/Code/Source/svFSI/CepModAp.cpp

  Lines 105 to 177 in b92add4

  	void CepModAp::integ_cn2(const int nX, Vector<double>& Xn, const double Ts, const double Ti, const double Istim, const double Ksac,
  	Vector<int>& IPAR, Vector<double>& RPAR)
  	{
  	int itMax = IPAR(0);
  	double atol = RPAR(0);
  	double rtol = RPAR(1);
  	double t = Ts / Tscale;
  	double dt = Ti / Tscale;
  	double Isac = Ksac * (Vrest - Xn(0));
  	double fext = (Istim + Isac) * Tscale / Vscale;
  	Xn(0) = (Xn(0) - Voffset) / Vscale;
  	auto Im = mat_fun::mat_id(nX);
  	Vector<double> fn(nX);
  	getf(nX, Xn, fn, fext);
  	//CALL AP_GETF(nX, Xn, fn, fext)
  	int k = 0;
  	auto Xk = Xn;
  	bool l1 = false;
  	bool l2 = false;
  	bool l3 = false;
  	t = Ts + dt;
  	double eps = std::numeric_limits<double>::epsilon();
  	while (true) {
  	k = k + 1;
  	Vector<double> fk(nX);
  	getf(nX, Xk, fk, fext);
  	//CALL AP_GETF(nX, Xk, fk, fext)
  	auto rK = Xk - Xn - 0.5*dt*(fk + fn);
  	double rmsA = 0.0;
  	double rmsR = 0.0;
  	for (int i = 0; i < nX; i++) {
  	rmsA = rmsA + pow(rK(i),2.0);
  	rmsR = rmsR + pow(rK(i) / (Xk(i)+eps), 2.0);
  	}
  	rmsA = sqrt(rmsA / static_cast<double>(nX));
  	rmsR = sqrt(rmsR / static_cast<double>(nX));
  	l1 = (k > itMax);
  	l2 = (rmsA <= atol);
  	l3 = (rmsR <= rtol);
  	if (l1 || l2 || l3) {
  	break;
  	}
  	Array<double> JAC(nX,nX);
  	getj(nX, Xk, JAC, Ksac*Tscale);
  	//CALL AP_GETJ(nX, Xk, JAC, Ksac*Tscale)
  	JAC = Im - 0.50 * dt * JAC;
  	JAC = mat_fun::mat_inv(JAC, nX);
  	rK = mat_fun::mat_mul(JAC, rK);
  	Xk = Xk - rK;
  	}
  	Xn = Xk;
  	getf(nX, Xn, fn, fext);
  	//CALL AP_GETF(nX, Xn, fn, fext)
  	Xn(0) = Xn(0)*Vscale + Voffset;
  	if (!l2 && !l3) {
  	IPAR(1) = IPAR(1) + 1;
  	}
  	}

  )
• Backwards Differentiation Formulas (variable-order, adaptive sub-stepping, no svFSIplus implementation yet)

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]

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.

[zasexton]

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]

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.

[zasexton]

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.

[zasexton]

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