98 ODE Solvers

The process models are nothing but systems of ordinary differential equations, presented in the form of mass balance in unit time:

Mass_Accumulation_Rate = Mass_In_Rate - Mass_Out_Rate + Mass_Net_Gen_Rate

V * (dCi/dt) = Flow * (Ci_in - Ci_out) + r_Ci * V, where

V is the process unit's active volume; Ci is the ith model component as concentration; Flow is the total flow rate of Ci; r_Ci is the net change/reaction rate of Ci; and t is time.

Therefore,

dCi/dt = Ci_in - Ci_out + HRT * r_Ci

The ASM, as well as other process models, both biological and phys-chem, define the r_Ci terms.

Solving the systems of ODEs can be both easy and challenging. It is easy because one can say, just like in many false advertisements for particular programming languages, that he/she solves a system of ODEs with a single line of code in the XYZ language or with the ABC library.

That is pure bullshit.

Tens of thousands of lines, produced by math experts, tested by many more in decades, were wrapped into that "one line of code". For ODE solvers, most of the trusted packages were written in FORTRAN and/or C/C++.

I spent much time last year learning the development of ODE solvers, attempted making my own Euler, RK4, RKF45, and partitioning the ODE system into "fast" and "slow" rate sections as suggested by the IWA ASM report. The ASM models are stiff. As a result, they are very challenging to get both high convergence rate and accuracy. All my efforts were proven slow in convergence due to the stability requirement because the methods I tried were mainly good for non-stiff ODEs.

The partioning technique introduced in the IWA ASM report would work for getting a steady state solution by relaxation. I doubt it would work for real time integration for dynamic situations. Therefore, I did not consider it a good candidate in the long run.

I haven't had time to work up a solver based on BDF, but will attempt at least once for sure. For now, I resorted to using scipy.integrate.solve_ivp() to move the project forward, feeling comfortable choosing the right integration sub-routines for the problems. I am reading the works by the ODE experts like Gears and hope to build one that is efficient for stiff ODEs.