Implicit and explicit (Super Time Stepping) solvers for diffusion equation in Energy Step

In order to solve the Energy Step corresponding to the radiative diffusion we use an implicit solver (sor). The reason is that, as it is well known, the numerical solution of the diffusion equation obtained via explicit algorithm is numerically unstable for time-step larger than:

\Delta t_p = {{ C_p }\over {max({{\Chi_x} \over {\Delta x^2}} + {{\Chi_y}\over {\Delta y^2}} + {{\Chi_z}\over {\Delta z^2}})}}

where C_p<0.5 is the parabolic courant number and \Chi_x, \Chi_y and \Chi_z are the thermal diffusivities in units [L^2]/[T].

We recall that in the code the hydro time-step computed via the CFL condition scales linearly with the grid resolution (\Delta x, \Delta y and \Delta z); therefore for large grid resolution and/or thermal diffusivities, \Delta t_p is much shorter than the hydro time step and the explicit method is very inefficient.

The implicit solution has the great advantage of being unconditionally stable (one may loose precision of the solution but not the stability for large time-steps). The great disadvantage is that implicit algorithm are based on iterative solutions and is not solvable in parallel in a scalable way.

There exists an alternative method, which is explicit and allows to use a time-step \tau larger than \Delta t_p iterating over a small number s of sub steps of length \Delta t_j such that:

\tau = \sum _{i=1}^{s} \Delta t_j.

This is obtained requiring that the solution is stable on \tau and not on each \Delta t_j. The method is called super-time stepping (Alexiades,Amiez,Gremaud 1995).

The solution at intermediate steps \Delta t_j is obtained either via Tchebychev or Legendre polynômes. An example for s=3 is given in Mignone et al. 2017

For our code fargOCA the interest of the super time stepping is twofold:

1. explicit method are easy to parallelize: we would really benefit of high degree of parallelization in this phase in which we are moving the code on Kokkos (practically finished by @alainm) in order to use GPUs and future accelerators. We stress the fact that new machines coming on the national panorama are based on GPUs.
2. we need to solve th diffusion equation faster for all the projects involving fully radiative model like for example the "inflow" case that we are studying with @morby