# Solving PDE implicitly or explicitly depending on stiffness

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?

• What you're trying to do would normally be called an implicit-explicit (IMEX) solver, for when the ODE has the form $\dot y=f_{\text{non-stiff}}+f_{\text{stiff}}$? (Also makes this question easier to search for.) Is there a specific reason you are implementing all this yourself instead of using an off-the-shelf library? – Kirill Jul 27 '18 at 14:56
• Thanks for that hint, I'll look into IMEX solvers. Well, the reason for implementing this by myself is that I have so many different parts, all with different numbers of nodes and/or dimensions, that I did not find any solver which could handle this. Besides that the derivatives of my parts have exactly the same shape as the value arrays of my parts. For most parts this means 1D or 2D with the x/y-length the same as the x/y node numbers. Reshaping this to a tri-/penta-/octa-diagonal matrix seemed to be more effort to me, than writing my own solver. Or is there any solver which can solve this? – Scotty1- Jul 27 '18 at 15:04

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.
• Thanks again for your answer Chris. I'll try to sum up my thoughts as short as possible. The julia package looks great and I'd really like to use it, but I am restrained to python. Besides that my PDEs are shaped in a way, which is, looking at the tutorial examples (not fully understanding all parts of julia, thus I might be wrong), not supported by the package. For each part in my simulation environment, I have a value array, the temperature, which is shaped $M\times N\times K$, where $M$ is the the number of nodes in $x$-direction, $N$ the number of nodes in $y$-direction etc... – Scotty1- Jul 30 '18 at 11:01
• Each part has a diff-function, which returns the differentials of the part for a given set of input parameters, most importantly the temperature at the nodes and the massflow in the nodes. These differentials have exactly the same shape $M \times N \times K$ as the value arrays. The coefficients of the PDEs, like the diffusivity, are temperature and massflow dependent and partially also dependent on the grid position. This is approximated by empirical relations like Nusselt equations, since simulating these parameters with CFD or others is impracticable. – Scotty1- Jul 30 '18 at 11:09
• In the simulation environment I have several 100 parts of different shapes and number of dimensions (most are 1D, thus having vectors as arrays, some are 2D and 3D). Adding these to a matrix of linear equations would be a tremendous effort, especially since the empirical equations for things like the diffusivity would need to be updated in each step. Furthermore having only one 3D part in a system of 1000 1D parts would require the solver to solve a septa-diagonal matrix for all parts. So is there any feasible way to pass a system of diff-functions which share the boundaries to a solver? – Scotty1- Jul 30 '18 at 11:18