Solving PDE implicitly or explicitly depending on stiffness

Question

I've got a system of several PDEs for a multitude of parts which represent real hydraulic parts like pipes or thermal energy storages. Each of these parts may have an arbitrary number of nodes and/or or dimensions.
The PDEs of most parts are non-stiff, while a few parts tend to have quite stiff PDEs.
Thus filling one sparse matrix with all of these parts and solving this matrix using one and the same solver, like embedded RK(4)5/Dormand-Prince or embedded Heun methods is not easily feasible.

Currently I solve the PDEs of each part for itself. For this example let's say using an embedded Heun method with adaptive stepsize control. Therefore each part has a function (in the simulation program) which returns its derivative for the given set of variables. This derivative is then used to calculate the intermediate steps of the Heun method for each part separately. Values shared between parts, like temperatures of contact areas, are updated in between stages AND steps. This results in the following programm control flow with the stepsize $h$:

1. Determine predictor step input variables for all parts
2. Predictor step:
   1. Get derivative of part 1 and calculate the intermediate result of it with $y_{interm,i+1} = y_i + h*deriv_{pred}(y_i)$
   2. Do the same for part 2
   3. ... until part n.
3. Update value array which stores the variables of contact areas
4. Corrector step:
   1. Get derivative of part 1 and calculate the final step result of it with $y_{i+1} = y_i + h/2*(deriv_{pred}(y_i) + deriv_{corr}(y_{interm,i+1}))$
   2. ... until part n.
5. Check if von Neumann stability was violated and if yes, reduce the stepsize $h$, reset value array for connection variables and go back to 2.
6. Get error estimate and if error too big, reduce stepsize $h$, reset value array for connection variables and go back to 2.

This approach works really well for most parts. But as soon as there is only one stiff part out of 1000 non-stiff-parts, the stepsize $h$ may easily be reduced to about $1/100$ of the "old" stepsize.
To deal with this, I thought that solving only the stiff parts with an implicit scheme may be a good approach, since there will be several hundred parts of which only 2 or three will be stiff. This would retain the fast calculation speed of the non-stiff parts as well as the "big" stepsize, at the cost of a higher computational cost for the stiff parts and an increased truncation error in these parts.
The new control flow would be:

1. Determine predictor step input variables for all parts
2. Predictor step:
   1. Get derivative of part 1 and calculate the intermediate result of it with $y_{interm,i+1} = y_i + h*deriv_{pred}(y_i)$
   2. Do the same for part 2
   3. ... until part n.
3. Get implicit PDE solution for stiff parts with half stepsize with $y_{imp,i+0.5} = y_i + h/2*deriv(y_{imp,i+0.5})$
4. Update value array which stores the variables of contact areas
5. Corrector step:
   1. Get derivative of part 1 and calculate the final step result of it with $y_{i+1} = y_i + h/2*(deriv_{pred}(y_i) + deriv_{corr}(y_{interm,i+1}))$
   2. ... until part n.
6. Check if von Neumann stability was violated and if yes, reduce the stepsize $h$, reset value array for connection variables and go back to 2.
7. Get error estimate and if error too big, reduce stepsize $h$, reset value array for connection variables and go back to 2.
8. Get final implicit PDE solution for stiff parts with half stepsize with $y_{imp,i+1} = y_{imp,i+0.5} + h/2*deriv(y_{imp,i+1})$

In this example I used Euler implicit method, but I'd preferably integrate something like Crank-Nicolson.

---

Will this be ok from a numerical/mathematical point of view? Is this a good idea? Or is there any other way to do it? Any concerns?
Thanks for your help in advance!

Answer 1

If you just slap together an implicit and an explicit method you will likely have order loss. You can do so with low order methods though, and Crank-Nicholson mixed with some other integrator is an easy way to get a decent second order integrator. Higher order IMEX integrators like the Kennedy and Carpenter Additive Runge-Kutta methods or the SBDF schemes are specifically designed to solve $u' = f_{stiff} + f_{nonstiff}$ at higher order (with adaptivity and all of that jazz). These are the kinds of methods you're looking for if you want efficiency.

There are two open source implementations that I know of. The Sundials ARKODE library and the DifferentialEquations.jl Split ODE Solvers both have these same higher order additive Runge-Kutta methods, and DifferentialEquations.jl also has the lower order ones plus the SBDF methods in development. I would recommend the DifferentialEquations.jl ones not just because I am the lead developer of the JuliaDiffEq organization (disclosure), but also because the ARKODE implementations tend to diverge without coefficient tweaks on any heavily stiff equation according to the benchmarks (for example, here's a test on OREGO, and there's a lot more tests here. These same tweaks are shown in the ARKODE C++ example code in order to make ROBER not diverge, so even the developers know), so while it's a part of Sundials it's not as robust as the more well-known CVODE integrator.