How to construct MX computation graph from existing SX functions wisely?

[程迪]

Hi everyone,

I am using casadi to implement a chemical kinetics to generate c code for reactor simulation. The reaction mechanism consists of N species and K reactions, N and K is about 100 to 1000. The state variable is temperature T and concentration c[N]. The reactor is controlled by an equation like:

dT/dt = Tdot(T,c)

dc/dt = cdot(T,c)

The most used constructing pieces are:

1.  Arrhenius formula:

kf(T;A,b,Ta) = A*b**T*exp(-Ta/T)

T is input

A,b,Ta is parameters per rate which can be tabulated as dense DM.

2. piecewise polynomial interpolation for thermodynamics

Cp(T;T_min,T_mid,T_max;cp_low_coeffs,cp_high_coeffs)=pw_const(T,[T_min,T_mid,T_max], [cp_min,cp(T;cp_low_coeffs),cp(T;cp_high_coeffs),cp_max])

cp(T;cp_low_coeffs) = # a polynomial of T, must use`((a1*T+a2)*T+a3)*T+a4)+logT*a5+a6/T` style to be stable. 

all parameters are dense can be tabulated

3. reaction rates

RR(T,c;F,R) = kf(T)*(c[1]**F1*c[2]**F2* ... )+kb(T)*(c[3]**R1*c[3]**R2)

F and R are highly sparse because the reaction is only related to only small set of species. Usually, F and R's elements are integers.

4. matrix multiplication

dc/dt = stoichimetric_matrix * RR

stoichimetric_matrix is a highly sparse matrix with almost 

5. some element-wise operations

There is small fraction of reaction rates must be modifed individually.

for example, for finite kinds of reactions called Falloff reactions, there are two Arrhenius rates which must be combined together using a complex formula related to almost all species concentration.

Using SX to implement it is straight forward, but the generated code is too big if I N and K is beyond 60. So I would like to re construct it using MX. 

My quesiton is is there some guides or best practice on how to use map() or other aggregation operations to construct the computational graph wisely, so the generated code would have a lot of for loops and can be effectively parallelized using SIMD-like ILP techniques rather than just OpenMP.

And here is some confusions I encountered:

1. How to use map function to generate F(T,coeffs[K]) from F(T,coeffs)? I only want to parallelize only a fraction of the input arguments.

2. How to construct a effective MX or SX function to modify only a fraction of elements in a vector.