Choosing the thermal conduction equation solution method in SPH

The set of explicit iterative schemes for parabolic problems (heat conduction as example) is reviewed. The iterative process is constructed based on algebraic polynomials (Chebyshev, Lanczos, or Legendre) properties. They combine the simplicity of explicit time-integration with the increased allowed timestep. This approach may be applied at smoothed particle hydrodynamics (SPH) simulations as a replacement for the popular implicit time-integration schemes, whose drawbacks might get stronger especially on SPH-like contrary to mesh-style space description.


Introduction
In SPH method the space differential approximation is rather wide (∼30-60 neighbours in 1D, ∼50-90 in 2D, and ∼80-200 in 3D). That's why the well-known implicit time-integration approach to the thermal conduction equation solution might become more complicated, comparing to the mesh methods with only cell boundary neighbours stencil.
On the other side, the simple and robust explicit time-integration method for parabolic equations suffer from the stability criteria for timestep τ : C = τ K h 2 ≤ 1 2 , where h is the typical space discretization resolution, K is the so-called kinematic diffusion coefficient with dimension L −2 t (L stands for unit of distance, t -unit of time).
That's why the idea to overcome the timestep limit while retaining the explicit scheme may bring extra efficiency on massively parallel architectures, where overheads of solvers at implicit algorithms are considerable.

Heat conduction notation in SPH Parabolic equation for heat conduction is
with ρ -is for density, e -specific (i.e. per unit mass) internal energy, κ -heat conduction coefficient, u -for temperature.
The simplest case will become the one-dimensional (1D) diffusion equation for u In SPH the 1D function g(x) space-differentiation operator at coordinate x i acts in another manner h eff means the effective size of "particle" (also thought as discretization point or material point) at coordinate x j , positive h i is the characteristic width of the dimensionless non-negative bell-shaped smoothing function w q i (x j ) is even function, therefore the support area is presented by the support radius (cut-off Application of space differentiation operators in a sequential manner to obtain the higherorder operator seems bad practice [1,2], due to necessity to take summation over second stage of neighbourhood (i.e. through neighbours k of neighbours j of particle i), and moreover the additional losses of accuracy on each summation.
Fortunately, some techniques of obtaining higher-order approximation of sequential differentiation operators in SPH were developed [3,4]. We utilize the universal approach to higher-order approximation of SPH-summation operators representation, where we can chose the desired order of approximation easily, transforming the sequential first-order operators to immediate (pure) versions of higher-order differential operators As the SPH space-differentiation operator stencil is rather wider, then at typical mesh operator (i.e. {x + i · h} in 1D, where i ∈ {−N, . . . , N }, N ∼ 1 in mesh case, and ∼ 2.5 ÷ 3 -in SPH), the parabolic equations implicit solution methods, widely used in mesh discretizations, seem less efficient applied to SPH. That's why we try to modify the standart explicit methods for use in SPH approach.

Methods developed by Kozyrev&Litvinov (KL) [5]
For linear parabolic equation ∂ ∂t u(x, t) = ∂ 2 ∂x 2 u(x, t) the numerical scheme with two time-layers is where superscripts n + 1 is for "t + τ ", n -for "t" time layers, G is the differential timestep operator.
For explicit scheme G → E + τ A, where E is the unitary operator, A -space discretization of parabolic operator ∂ 2 ∂x 2 : For numerical scheme stability max 0≤λ≤4C |P Q (λ)| ≤ 1 is necessary and sufficient.

based on Chebyshev polynomials
Guess the linear representation To write the scheme as predictor-corrector, reform numerical scheme from [5] based on Chebyshev polynomials (KL-C):

based on Lanczos polynomials
Guess the linear representation λ = Q(Q+2) To write the scheme as predictor-corrector, reform Iterative (5) numerical scheme from [5] based on Laczos polynomials (KL-L): 4. "Locally iterative" (LI) numerical schemes [6] LI and LI-M (i.e. "modified") schemes for (5):   To obtain the LI-type scheme of minimal variation, the sequence of Chebyshev polynomial roots should be reordered in a special manner. This ensures stability for non-linear problems.
The O(t 2 ) LI-2 scheme can be constructed in a predictor-corrector manner from the LI-M one, we just apply the LI-M for the τ 2 as a predictor step, and then perform the corrector step:

Runge-Kutta-Legendre (RKL) methods [7]
Super-time-stepping methods fall into a category of stabilized Runge-Kutta methods for which additional stages are used to increase the stability of the method, thereby increasing the maximum permissible time-step. This should be contrasted with traditional Runge-Kutta methods, where the additional stages are used to increase the accuracy of the method.
. The internal Q-stage update of (3) can be expressed in terms of the stability polynomial , RKL methods are based on using shifted Legendre polynomials as the stability polynomial of the scheme: Stability is ensured if ∀λ : |P Q (λ)| ≤ 1.

O(t) scheme
We require that every internal stage also has a corresponding stability polynomial Thus, each stage in the scheme can be thought of as a first-order accurate approximation to the solution at "inner" time t = ν(ν+1) Q(Q+1) τ . O(t) Runge-Kutta numerical scheme based on Legendre polinomials (RKL1), (Q = 1 → explicit): Chebyshev T ν (ξ) = cos (ν arccos ξ) Runge-Kutta numerical scheme based on Legendre polynomials (RKL2), 1 < Q : 6. Travelling heat wave problem simulation Simulation setup is:  scheme (eq.) steps quantity Q ≥ 1 steps time-layers final iteration time consistency C 1 C = 0.5 ν n order on iter (10) We chose x L = 0, x R = 1, D = 4 and run simulations till t end = 0.2 in two cases: with s = 4 and s = 8. Timestep τ was defined from the relation τ D h = α, where h = 10 −2 is the mesh resolution. And we vary the amount of mesh cells the thermal wave propagates per τ choosing α as 0.5 or 1.
The obtained results of simulations are presented on figures 3-4 (horizontal axis -x, vertical