Scalable explicit implementation of anisotropic diffusion with Runge-Kutta-Legendre super-time-stepping

TL;DR

This paper presents a second-order Runge-Kutta-Legendre super-time-stepping scheme for efficiently and accurately simulating anisotropic thermal conduction in astrophysical plasmas, outperforming previous methods in speed and precision.

Contribution

It introduces a robust, fast, and second-order accurate explicit scheme for anisotropic diffusion, improving over existing super-time-stepping methods in stability and accuracy.

Findings

RKL scheme is more accurate than first-order STS methods.

Explicit conduction with RKL scales efficiently on large parallel architectures.

The method effectively handles saturated thermal conduction in simulations.

Abstract

An important ingredient in numerical modelling of high temperature magnetised astrophysical plasmas is the anisotropic transport of heat along magnetic field lines from higher to lower temperatures.Magnetohydrodynamics (MHD) typically involves solving the hyperbolic set of conservation equations along with the induction equation. Incorporating anisotropic thermal conduction requires to also treat parabolic terms arising from the diffusion operator. An explicit treatment of parabolic terms will considerably reduce the simulation time step due to its dependence on the square of the grid resolution ($\Delta x$) for stability. Although an implicit scheme relaxes the constraint on stability, it is difficult to distribute efficiently on a parallel architecture. Treating parabolic terms with accelerated super-time stepping (STS) methods has been discussed in literature but these methods suffer from poor accuracy (first order in time) and also have difficult-to-choose tuneable stability parameters. In this work we highlight a second order (in time) Runge Kutta Legendre (RKL) scheme (first described by Meyer et. al. 2012) that is robust, fast and accurate in treating parabolic terms alongside the hyperbolic conversation laws. We demonstrate its superiority over the first order super time stepping schemes with standard tests and astrophysical applications. We also show that explicit conduction is particularly robust in handling saturated thermal conduction. Parallel scaling of explicit conduction using RKL scheme is demonstrated up to more than $10^4$ processors.

Tables

Table 1. Table 1: Expected Flops relative to MHD evolution

MHD	explicit	sub-cycling	STS	implicit†
1	$\left(1+r^{-1}\right)t^{-1}$	$1+t^{-1}{r}^{-1}$	$1+t^{-1/2}{r}^{-1}$	$1+r^{-1}$

Table 2. Table 2: A comparison of the execution time for the blast wave problem

$N$	Wall-clock time $\left[s\right]$
	explicit	RKL-2
$192$	$292$	$112$
$384$	$4448$	$1276$
$768$	$62700$	$11828$
$1536$	$1.11\times {10}^6$	$1.45\times {10}^5$

Table 3. Table 3: Wall-time for the local TI test run till 0.87 Gyr

explicit	sub-cycling	AAG-STS†	RKL-2
6 h 1 m 29 s	43 m 45 s	9 m	16 m 6 s

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})$$

$$\frac{\partial }{\partial t}\left(\rho \mathbf{v}\right)+\nabla \cdot \left(\rho \mathbf{v}\mathbf{v}-\mathbf{BB}+p_t\mathsf{I}\right)$$

$$\frac{\partial E}{\partial t}+\nabla \cdot \left[\left(E+p_t\right)\mathbf{v}-(\mathbf{B}\cdot \mathbf{v})\mathbf{B}\right]$$

$$\frac{\partial \mathbf{B}}{\partial t}-\nabla \times (\mathbf{v}\times \mathbf{B})$$

$${\mathbf{F}}_{\mathrm{c}}=\frac{q}{q+F_{\mathrm{class}}}{\mathbf{F}}_{\mathrm{class}},$$

$${\mathbf{F}}_{\mathrm{class}}=-{\kappa }_{\parallel }\hat{\mathbf{b}}(\hat{\mathbf{b}}\cdot \nabla T)$$

$$q=5\varphi \rho c_{\mathrm{iso}}^3.$$

$${\mathbf{F}}_{\mathrm{class}}=-\kappa \nabla T.$$

$$\frac{\partial (\rho ϵ)}{\partial t}=-\nabla \cdot {\mathbf{F}}_c,$$

$$-\nabla \cdot {\mathbf{F}}_c\approx -\sum _d\frac{({\hat{\mathbf{e}}}_d\cdot {\mathbf{F}}_c{)}_{i+\frac{1}{2}}-({\hat{\mathbf{e}}}_d\cdot {\mathbf{F}}_c{)}_{i-\frac{1}{2}}}{\Delta h_d}$$

$$\Delta t=\min\Big{(}\Delta t_{\rm h},\,\Delta t_{\rm p}\Big{)}\,,$$

$$\Delta t_{\mathrm{h}}=\frac{C_h}{\max \left(\frac{{\lambda }_x}{\Delta x}+\frac{{\lambda }_y}{\Delta y}+\frac{{\lambda }_z}{\Delta z}\right)}$$

$$\Delta t_{\mathrm{p}}=\frac{C_p}{\max \left(\frac{{\chi }_x}{\Delta x^2}+\frac{{\chi }_y}{\Delta y^2}+\frac{{\chi }_z}{\Delta z^2}\right)}$$

$$\chi_{d}=\left\{\begin{array}[]{ll}\displaystyle\kappa\frac{\gamma-1}{\rho}&\qquad\mathrm{(HD)}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\kappa_{\parallel}b_{d}\cos\theta\frac{\gamma-1}{\rho}&\qquad\mathrm{(MHD)}\end{array}\right.$$

$$\frac{du}{dt}=\mathsf{M}u,$$

$$u(t+\tau )=R_s(\tau \mathsf{M})u(t).$$

$$R_s(\tau \lambda )\equiv \prod _{i=1}^s(1+\lambda \Delta t_i)$$

$$u(t+\tau )=\left(1+\tau \mathsf{M}+\frac{1}{2}(\tau \mathsf{M}{)}^2+...\right)u(t).$$

$$R_s(\tau \lambda )=a_s+b_s{\mathcal{P}}_s\left(w_0+w_1(\tau \lambda )\right)$$

$${\mathcal{P}}_j\left(x\right)=\left(\frac{2j-1}{j}\right)x{\mathcal{P}}_{j-1}\left(x\right)-\left(\frac{j-1}{j}\right){\mathcal{P}}_{j-2}\left(x\right).$$

$$\begin{array}[]{lcll}\Delta t_{\rm h}&=&\displaystyle\Delta t_{\rm p}\frac{(s^{2}+s)}{2}&\text{for RKL-1,}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\Delta t_{\rm h}&=&\displaystyle\Delta t_{\rm p}\frac{(s^{2}+s-2)}{4}&\text{for RKL-2}.\end{array}$$

$$\Delta t_{\mathrm{h}}=\Delta t_{\mathrm{p}}\frac{s}{2\sqrt{\nu }}\left[\frac{(1+\sqrt{\nu }{)}^{2s}-(1-\sqrt{\nu }{)}^{2s}}{(1+\sqrt{\nu }{)}^{2s}+(1-\sqrt{\nu }{)}^{2s}}\right].$$

$$C_p\to \frac{s^2}{4}(\mathrm{RKL}-1),C_p\to \frac{s^2}{8}(\mathrm{RKL}-2).$$

$$C_{p}\to\left\{\begin{array}[]{lcl}\displaystyle\frac{s^{2}}{2}&\quad\mathrm{for}&s\lesssim 1/(2\sqrt{\nu}).\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\frac{s}{4\sqrt{\nu}}&\quad\mathrm{for}&s\gtrsim 1/(2\sqrt{\nu}).\end{array}\right.$$

$$\rho ϵ=\frac{p}{\gamma -1}$$

$$\frac{\partial (\rho ϵ)}{\partial t}=-\nabla \cdot {\mathbf{F}}_c,$$

$$\frac{\partial }{\partial t}T(x,t)=\frac{\partial }{\partial x}\left(\kappa \frac{\partial T(x,t)}{\partial x}\right),$$

$$T^{\mathrm{ref}}(x,t)=\frac{1}{\sqrt{1+\frac{2\kappa t}{{\sigma }^2}}}\exp \left[-\frac{x^2}{2{\sigma }^2\left(1+\frac{2\kappa t}{{\sigma }^2}\right)}\right].$$

$${ϵ}_{L1}=\frac{1}{N_x}\sum _i\left|T_i^n-T^{\mathrm{ref}}(x_i,t^n)\right|.$$

$$\frac{\Delta t_{\mathrm{h}}}{\Delta x}=\frac{C_p\Delta x}{\kappa }=\frac{C_{p0}}{N_{x0}\kappa }=\mathrm{const},$$

$$T(x,0)=10+20\left[\frac{x}{10}-\mathrm{floor}\left(\frac{x}{10}+\frac{1}{2}\right)\right]$$

$$N_{\mathrm{step}}=\frac{t_{\mathrm{stop}}\kappa }{C_p}{\left(\frac{N_x}{L_x}\right)}^2=\frac{1600}{C_p},$$

Full text

1 Introduction

Since most baryons in the universe are in a magnetised plasma state, magnetic fields play a crucial role in the dynamics and thermodynamics of astrophysical objects — ranging from stars and interstellar medium to the intra-cluster and intergalactic medium. Magnetohydrodynamic (MHD) simulations have matured (Evans & Hawley 1988; Tóth 2000; Balsara 2001) and have contributed to several breakthroughs in our understanding, from accretion to the interstellar medium (Hawley et al. 1995; Korpi et al. 1999). Magnetic fields not only produce forces and stresses in a plasma, they also affect transport properties by predominantly allowing diffusion of heat and momentum along field lines and suppressing transport across them. Anisotropic transport affects fundamental properties such as convection/buoyancy in stratified plasmas (Balbus 2000; Quataert 2008), and thermal instability and condensation of cold gas out of the hot phase (Field 1965; Sharma et al. 2010b).

While the numerical solution of the ideal MHD equations can be carried out very accurately in highly nonlinear regimes across thousands of processors (Beresnyak 2011; Federrath et al. 2011), simulations with anisotropic thermal conduction (and, likewise, similar diffusive processes) have not yet reached the same level of fidelity and scalability.

Taking into account diffusion processes changes the mathematical structure of the underlying system of conservation laws from purely hyperbolic to mixed hyperbolic-parabolic type. The numerical discretisation of such systems can then pose a more restrictive limitation on the choice of the time step which, for a standard explicit scheme, must scale with the square of the grid size rather than with $\Delta x$ alone. In addition, when saturation effects are considered, thermal conduction itself becomes a mixed hyperbolic/parabolic operator. Two approaches are commonly employed to circumvent the stability time step constraint.

In the first one, the parabolic part of the equations is solved using fully- or semi-implicit methods (Press et al. 1986), allowing for a time step much longer than the explicit one (see e.g., Sharma & Hammett 2011). However, scalable implementation of implicit schemes on massively parallel clusters is technically challenging (Kannan et al. 2015) as it requires solving large sparse matrices.

In the second approach, a particular stability polynomial can be employed to construct an explicit multi-stage time stepping scheme that has extended stability properties. If $s$ is the number of stages, the length of the super-step that one can take (asymptotically) can be shown to scale as $\sim s^{2}\Delta t_{\rm p}$, where $\Delta t_{\rm p}$ is the standard explicit parabolic time step. This technique in general is classified as super-time-stepping (STS). This technique offers therefore an overall speedup $\propto s$ and, being explicit, it can be easily parallelized. van Der Houwen & Sommeijer 1980; Verwer et al. 1990 introduced Runge-Kutta Chebyshev (RKC-STS) methods. A variant of this method was popularized by Alexiades et al. 1996 (AAG-STS). More recently, super-time-stepping based on Legendre polynomials, known as Runge-Kutta Legendre STS (RKL-STS), was proposed by Meyer et al. (2012, 2014).

In this paper we review and compare the various STS schemes in the context of anisotropic thermal conduction in a MHD plasma. We show that RKL-STS method is an attractive option for simulating anisotropic diffusion on massively-parallel supercomputers and offers enhanced stability with respect to its predecessors, in particular, AAG-STS. We note that RKL-STS can also be used to model other isotropic and anisotropic transport processes such as plasma viscosity (e.g., Dong & Stone 2009), non-ideal terms in Ohm’s law (e.g.,O’Sullivan & Downes 2007) , and cosmic ray diffusion (e.g., Sharma et al. 2009; Pakmor et al. 2016). There are also applications in other diverse areas such as image processing (e.g., Weickert 1998).

The paper is structured as follows. In section 2 we review the fundamental equations of MHD with thermal conduction. In Section 3 we review the numerical approaches to anisotropic thermal conduction and also discuss the monotonicity properties in multiple dimensions. Selected numerical benchmarks are presented in 4. The parallel scalability of explicit conduction is demonstrated in 5 and conclusions are drawn in section 6.

2 Governing Equations

Here we consider the magnetohydrodynamic (MHD) equations in presence of conduction:

[TABLE]

where $\rho$ is mass density, $\bm{v}$ is fluid velocity, $\bm{B}$ is magnetic field, $p_{t}=p+B^{2}/2$ is the total (gas + magnetic) pressure, $E=\rho\epsilon+\rho v^{2}/2+B^{2}/2$ is the total energy density while $\mathsf{I}$ is the identity tensor.

Thermal conduction effects, important for dilute gases and plasmas, are accounted for by the additional term on the right hand side of the total energy equation (1c). The thermal conduction flux $\bm{F}_{c}$ smoothly varies between classical and saturated regimes,

[TABLE]

where $F_{\rm class}$ and $q$ are, respectively, the magnitudes of the classical thermal conduction flux (Cowie & McKee, 1977; Balsara et al., 2008)

[TABLE]

and the saturated flux

[TABLE]

In Eq. (3), $\kappa_{\parallel}$ is the thermal conductivity parallel to the direction of the magnetic field (a function of temperature; for expressions see pp. 38 of Huba 2004) and $\bm{\hat{b}}\equiv\bm{B}/B$ is the magnetic field unit vector. We do not consider thermal conduction perpendicular to magnetic field lines because it is typically much smaller, but including it is straightforward. In absence of magnetic fields, the classical heat flux is simply given by (with $\kappa=\kappa_{\parallel}$)

[TABLE]

Saturation of heat flux (Eq. 4) comes into play when the mean free path of the electrons is large compared to the temperature scale-height $L_{T}=T/|\nabla T|$ (Cowie & McKee, 1977). The saturated heat flux is given by $q\bm{\hat{F}}_{\rm class}$ with $q$ specified in Eq. (4), while $\bm{\hat{F}}_{\rm class}=\bm{F}_{\rm class}/F_{\rm class}=-\bm{\hat{b}}~{}{\rm sgn}(\bm{\hat{b}}\cdot\nabla T)$ is a unit vector along the local magnetic field line but down the temperature gradient. In Eq. (4) the parameter $\phi$ is an uncertainty factor of the order of unity (taken to be 0.3) while $c_{\rm iso}\equiv(p/\rho)^{1/2}$ is the isothermal sound speed. In the large temperature gradient limit, thermal conduction is described by a hyperbolic operator. On the contrary, for $F_{\rm class}/q\to 0$, thermal conduction is described by a parabolic (diffusion) operator. The flux given in Eq. (2) is a harmonic mean of the fluxes in the two regimes (weighted toward the smaller of the two) and reflects the mixed parabolic/hyperbolic mathematical nature of the underlying differential operators. Eqs. 2.10 & 2.10 in Balsara et al. (2008) gives some other ways of combining saturated and classical heat fluxes. For more details on our implementation of saturated conduction see Appendix A of Mignone et al. (2012).

Note that thermal conduction described by Eq. (2) is very similar to the cosmic ray streaming equation (see Eqs. 1.1, 1.2 & 6.1 in Sharma et al. 2010a) in the sense that at temperature extrema, where gradient is zero and heat flux is not saturated, the behaviour is diffusive but at other points (in the long mean free path regime) it is hyperbolic. The harmonic mean used in Eq. (2) is analogous to the regularisations of the cosmic ray streaming equation at extrema proposed in Sharma et al. (2010a).

3 Numerical Approaches for Anisotropic Thermal Conduction

In this section we review the various methods for implementing anisotropic thermal conduction.

We shall use the PLUTO code (Mignone et al., 2007, 2012) to solve Eqs. (1a-1d) with the thermal conduction flux given by Eq. (2). In PLUTO, STS schemes have been implemented in an operator split fashion so that we first advance Eqns. (1a-1d) without the conduction flux using standard Godunov-type methods, followed by the solution of

[TABLE]

by means of the STS/explicit approach.

A conservative approach is used during both steps so that the same flux is used for adjacent grid cells sharing a face. During the parabolic step, Eq. (6) is discretised using standard second-order finite differences:

[TABLE]

where $d=x,y,z$ spans across dimensions, $\Delta h_{d}=(\Delta x,\,\Delta y,\,\Delta z)$ and $\bm{\hat{e}}_{d}$ are unit vectors in the three directions. The interface flux is computed using central differences for the temperature gradient and an upwind scheme for the saturated heat term, as shown in appendix A of Mignone et al. (2012).

3.1 Explicit methods

In the standard explicit approach (forward in time centred in space) the time step required to advanced the full system (1a)-(1d) is given by the shortest timescale:

[TABLE]

where

[TABLE]

is the MHD/hydro (hyperbolic) time step whereas

[TABLE]

is the diffusion time scale. In Eqs. 9 and 10, the maximum is taken over the whole computational domain, $C_{h}\leq 1$ is the hyperbolic Courant number, $C_{p}\leq 1/2$ is the parabolic Courant number, $\lambda_{d}=|v_{d}|+c_{d,f}$, $c_{d,f}$ is the fast magneto-sonic speed in the $d$ direction, $\Delta x,\Delta y$ and $\Delta z$ are the mesh spacings in the three directions while $\chi_{x},\chi_{y}$ and $\chi_{z}$ are the thermal diffusivities with the dimensions of $[{\rm length}]^{2}/[{\rm time}]$. Adopting an ideal EoS ($p=[\gamma-1]\rho\epsilon$) one obtains, from Eq. (6) together with Eq. (5) or Eq. (3):

[TABLE]

where $b_{d}=B_{d}/B$ is the component of the magnetic field unit vector in the direction $x_{d}$ and $\theta$ is the angle between $\bm{B}$ and $\nabla T$ while $d=x,y,z$ (note that, owing to operator splitting, density can be considered constant during the diffusion step).

Various explicit methods are discussed in detail in section 2 of Sharma & Hammett (2007). A notable symmetric explicit scheme with a very small numerical diffusion perpendicular to the local field lines is described in Günter et al. (2005). This method has been used widely in modelling fusion plasmas, in which the temperature and field gradients are relatively gentle. However, for highly nonlinear problems with large temperature gradients, this method gives unphysical temperature oscillations (see top right panel in Fig. 6 of Sharma & Hammett 2007).

It must emphasised that, for large grid resolutions and/or thermal conductivities, the diffusion timescale can become much shorter than the MHD/hydro time step (e.g., in galaxy clusters; Wagh et al. 2014; Yang & Reynolds 2015), that is $\Delta t_{\rm p}\ll\Delta t_{\rm h}$. A fully explicit method can thus become very inefficient. In these cases the ability to take larger steps than the stability limit (Eq. 8) is highly desirable.

A slight improvement over explicit time-stepping is obtained via sub-cycling, in which the MHD module is evolved with $\Delta t_{\rm h}$ but the conduction module is applied for a number of sub-steps ($N_{\rm sub}\approx\Delta t_{\rm h}/\Delta t_{\rm p}$; assuming $\Delta t_{\rm p}<\Delta t_{\rm h}$), each using the stable time step $\Delta t_{\rm p}$.While this approach is faster than the fully explicit approach, as we show in section 3.3, we can do much better than this using super-time-stepping.

3.2 Implicit methods

Implicit methods are very attractive because their time step is not limited by the stability limit (Eq. 8). A larger time step only leads to a slight loss in accuracy. Also, for very high resolution simulations (i.e., with very small $\Delta x$), the CFL time step (Eq. 9) becomes much longer than the stability time step with diffusion (Eq. 10). For this reason implicit and semi-implicit schemes for diffusion are very popular (e.g., Sharma & Hammett 2011; Dubois & Commerçon 2016; Kannan et al. 2015; Pakmor et al. 2016).

The key problem with implicit methods is that they involve solving a sparse matrix equation, which is difficult to solve in parallel in a scalable way. Even a tridiagonal matrix equation which results from the semi-implicit method of Sharma & Hammett (2011) is non-trivial to solve in parallel. The MHD codes using implicit diffusion typically use sparse matrix libraries such as Hypre111http://acts.nersc.gov/hypre/ or PETSc222http://www.mcs.anl.gov/petsc/ to solve the matrix equation iteratively. This global solution approach is fundamentally different from the explicit update in which information only propagates across nearest grid cells in a single time step. Most importantly, even the scalable implicit methods for diffusion with MHD do not scale on more than $\sim$1000 processors on massive distributed memory supercomputers (e.g., see Fig. 11 in Caplan et al. 2017). In contrast, the STS methods that we advocate (in particular, RKL-STS) in this paper, being explicit, are trivial to parallelise and show good scaling up to tens of thousands of processors (see section 5).

3.3 Super-time-stepping: RKC, RKL, AAG

Super-time-stepping (STS) is a way of choosing explicit time steps that are on-average much longer than the explicit stability time step (Eq. 10). The anisotropic diffusion equation has non-positive eigenvalues; i.e., none of the Fourier modes grow with time. Following Meyer et al. 2014 (see their section 2), we can write the diffusion equation in a semi-discrete ODE form

[TABLE]

where, for a linear problem, $\mathsf{M}$ is a symmetric constant coefficient matrix which results from the discretization of the parabolic operator. The eigenvalues $\lambda$ of such an operator are non-positive and real.

The STS method can indeed be thought of as a multi-stage Runge-Kutta method in which the intermediate stages are chosen for stability, rather than for a higher order accuracy. The update of an $s-$stage scheme is quantified in terms of an amplification factor, $R_{s}$, defined as

[TABLE]

We now define a stability polynomial for an $s$-stage STS scheme,

[TABLE]

with $\sum_{j=1}^{s}\Delta t_{j}=\tau$, and we impose $|R_{s}(\tau\lambda)|\leq 1$ for all values of $\lambda$ between the largest negative eigenvalue of $\mathsf{M}$ and 0 in order to ensure stability. Temporal accuracy is achieved by matching terms in the stability polynomial with the expansion of the analytic solution of Eq. (12),

[TABLE]

For first (second) order STS the stability polynomial only matches the above analytic expression up to the linear (quadratic) term.

The sub-step sequence in STS was originally chosen by writing the stability polynomial, $R_{s}$ (a polynomial of degree $s$ in $\lambda$), in form of a shifted Chebyshev polynomial (CP; e.g., Verwer et al. 1990); we refer to this scheme as RKC (Runge-Kutta Chebyshev). The sub-steps in a super-step are chosen to exploit the recursion properties of CPs. Since the absolute value of CPs is always $\leq 1$ if their argument lies in $[-1,1]$, the magnitude of amplification factor after a super-time step is $\leq 1$. Since CPs attain unity within $[-1,1]$, RKC-STS requires damping for the method to work for practical diffusion problems.

Recently, Meyer et al. (2012, 2014) have proposed Legendre polynomials (LPs) as the basis for the construction of robust STS schemes.

For a general $s$-stage RKL scheme, the stability polynomial is chosen to be of the form

[TABLE]

where $w_{0}=1$ is chosen for all RKL schemes. The LPs, $\mathcal{P}_{j}\left(x\right)$ obey a three point recursion relationship given by

[TABLE]

RKL-2 method enforces the stability polynomial at sub-stage $j$ to be $a_{j}+b_{j}\mathcal{P}_{j}\left(\tau\mathsf{M}\right)$ and uses the property of LPs that their absolute value is bounded by unity if their argument lies in (-1,1). Following Meyer et al. (2014), in Appendix A we illustrate the second order temporal accuracy of the RKL-2 scheme along with explicit formulae for an $s=3$ stage scheme.

The key advantage of RKL schemes over RKC is that they are more robust because LPs, unlike CPs, are always smaller than unity in magnitude if their argument lies in (-1,1). Thus, there is no need of an explicit damping parameter ($\nu$ in Eqs. 9 & 10 of Alexiades et al. 1996). Moreover, Meyer et al. (2014) show that RKL has superior linear stability and monotonicity properties compared to RKC.

For stability, the argument of the LP in Eq. (15) should be $\geq-1$ (note that for a parabolic operator the eigenvalues are non-positive); i.e., $1-w_{1}\tau|\lambda|_{\rm max}\geq-1$ or $\tau\leq 2/(|\lambda|_{\rm max}w_{1})$, where $|\lambda|_{\rm max}$ is the eigenvalue of $\mathsf{M}$ with largest absolute value. Note that $w_{1}=2/(s^{2}+s)$ and $4/(s^{2}+s-2)$ for RKL-1 and RKL-2, respectively. Given a super-time-step that can be conveniently chosen as the hyperbolic CFL time step $\Delta t_{\rm h}$ (Eq. 9), the number of sub-stages in RKL-STS is given by (see Eq. 20, 21 in Meyer et al. 2014)

[TABLE]

where, abbreviations RKL-1 and RKL-2 stand for the first and second order Runge-Kutta-Legendre schemes.We assume here that $s>1$; $s=1$ corresponds to a standard forward Euler update. In this paper, RKL without a qualifier refers to RKL-2. We explicitly mention if we are using the first order method. It is possible to make the scheme even faster by taking multiple $\Delta t_{\rm h}$ as the size of a super-step but of course one needs to worry about accuracy.

Likewise, for AAG-STS, the largest number of first order Euler sub-steps $s$ satisfies (see Eq. 10 in Alexiades et al. 1996)

[TABLE]

The maximum $s$ for STS schemes gives a speed-up of the conduction module scaling as $(\Delta t_{\rm h}/\Delta t_{\rm p})^{1/2}$. The speed-up is marginally better for RKC (compare Eqs. 18, 19 in Meyer et al. 2012) although RKL has better stability (as shown in section 4), and hence is much more attractive.

It is instructive to compute the effective parabolic CFL number $C_{p}=\chi\Delta t/\Delta x^{2}=\Delta t_{\rm h}/2\Delta t_{\rm p}$ as a function of the number of sub-steps $s$. In the limit of large $s$ ($s\gg 1$), the RKL schemes (17) yield

[TABLE]

However, in the same limit, AAG-STS yields an effective gain which depends on $\nu$:

[TABLE]

This equation sets an upper limit on the gain that can be reached using the traditional AAG-STS scheme: no substantial gain is obtained when $s\gtrsim 1/(2\sqrt{\nu})$.

As mentioned earlier, the key advantage of STS methods is their simplicity, parallel scaling on massively distributed clusters, and their large speed-up compared to the explicit method. Table 1 shows the expected scaling of Flops (equivalently, run-time) using various conduction schemes relative to pure MHD evolution. Here, $t=\Delta t_{\rm p}/\Delta t_{\rm h}$ (assumed to be $<1$) and $r$ is the ratio of Flops (floating point operations) for MHD and conduction modules. The typical value of $r$ is a few (say 5; this is relatively independent of the physical set-up and dimensionality) but $t$ can be smaller than 0.01 (this depends on the physical set-up). For these typical values the ratios of Flops required in Table 1 (for MHD, explicit, sub-cycling, STS, implicit) become: 1, 120, 21, 3, 1.2. Therefore, STS promises a good speed-up compared to the explicit method, and is competitive relative to implicit methods (especially given the better parallel scaling of the former).

3.4 Monotonicity in Multi-Dimensions

A serious problem for anisotropic diffusion, in presence of large gradients, is that the temperature can behave non-monotonically. The use of limiters for interpolating transverse temperature gradients has proved useful in preventing negative temperatures in such cases (Sharma & Hammett 2007; see Appendix B for details).

The use of limiters to calculate the transverse temperature gradients (Eq. 51) has shown to improve the robustness of all schemes (e.g., see Sharma & Hammett 2007, 2011). However, we show in section 4.1.3 that limiters reduce the accuracy of STS schemes for a large number of sub-steps. Given this, we keep limiters as an option to be used only when temperature gradients are large and the wrong sign of the heat flux gives a negative temperature in some grid cells.

4 Numerical Tests: effectiveness of RKL-STS

In the following we present selected numerical benchmarks. We adopt an ideal equation of state so that

[TABLE]

where $\gamma$ is the ratio of specific heats which equals $5/3$ unless otherwise stated.

4.1 Scalar Diffusion Equation

In this section we consider numerical tests based on the solution of the scalar diffusion equation

[TABLE]

where $\bm{F}_{c}$ is given by Eq (2). Other fluid variables are not evolved in time.

In all computations, we assign the super-step $\Delta t_{\rm h}$ and compute the number of sub-steps by solving either Eq. 17 (for RKL) or Eq. 18 (for AAG-STS) for $s$ and then rounding the solution to the next larger integer, i.e., $s\,\to\,1+{\rm floor}(s)$.

4.1.1 Gaussian diffusion in 1-D

A standard one dimensional diffusion test is useful to compare the stability and accuracy of different methods described in section 3.3. For this purpose, we solve the standard heat equation (without saturation) for temperature:

[TABLE]

using a constant diffusion coefficient $\kappa=1$. In the PLUTO code, this is equivalent to solving Eq. 6 by setting $\rho=1$ and $\gamma=2$, so that $\rho\epsilon=p=T$ in code units.

The initial condition consists of a Gaussian temperature profile $T(x,0)=\exp(-x^{2}/2\sigma^{2})$ inside the domain $x\in[-2,2]$. Then Eq. 23 admits the well-known analytical solution

[TABLE]

Boundary conditions are specified using the exact solution.

Figure 1 shows the computed temperature profile (symbols) versus the analytical solution (Eq. 24; solid line) at different times using RKL-2 on $N_{x}=128$ grid points and a parabolic CFL number $C_{p}~{}(\equiv\Delta t_{\rm h}/2\Delta t_{\rm p})=10$. As expected, the initial peak of the Gaussian temperature profile spreads symmetrically about $x=0$ with time.

In order to verify the order of accuracy of the various methods considered here, we first perform a resolution study by doubling the grid resolution from $N_{x}=64$ up to $N_{x}=16384$. We quantify the accuracy for different integration methods, AAG-STS (we use $\nu=10^{-3}$ unless otherwise stated), RKL-1, RKL-2, by estimating the error in $L1$ norm as

[TABLE]

Since we use standard centred differencing for the second derivative in Eq. 23 with local error $O(\Delta x^{2})$, we must ensure that the error in our simulation is not dominated by spatial discretisation. To this end, we adjust the time step such that the ratio of the time step and the grid spacing is fixed; i.e.,

[TABLE]

where $N_{x0}=64$ and $C_{p0}=0.45$ are, respectively, the number of zones and the parabolic CFL number at the lowest resolution. This yields a local truncation error (for $T$) of $O(\Delta t^{m})+O(\Delta t^{3})$ where $m=2,3$ for a first or second order scheme, respectively. Notice also that the parabolic CFL number increases linearly with the spatial resolution.

The left panel in Fig. 2 shows the errors obtained for AAG-STS, RKL-1 and RKL-2 methods together with the expected first and second-order convergence rates. From the figure one can verify that the order of accuracy for AAG-STS and RKL-1 is that of a $1^{\rm st}$ order method whereas RKL-2 converges as a $2^{\rm nd}$ order method, as expected. We remark that integration with AAG-STS shows the occurrence of negative temperatures during sub-steps for $N_{x}\geq 8192$ (corresponding to $s\gtrsim 30$), although the solution remains positive at the end of the super-step. This can pose serious difficulties in more complex applications in which the diffusion coefficient is a nonlinear function of the temperature (e.g., Spitzer conductivity $\kappa\propto T^{5/2}$) requiring $T$ to be non-negative at all times. On the contrary, RKL-1 and RKL-2 never exhibit such a behaviour (i.e., the solution remains positive for all sub-steps; this property is enforced by construction in the RKL scheme as described at the end of section 2.2 of Meyer et al. 2014).

In the middle panel of Figure 2 we plot the computational time (measured as the total number of steps and sub-steps) as a function of the number of grid points $N_{x}$. For an explicit scheme, let $N_{t}\propto N_{x}^{2}$ be the total number of steps required to reach some final time step. Then, from Eqs. 17 & 20, we expect the number of sub-steps using STS schemes to scale roughly as $s\propto\sqrt{\Delta t_{\rm h}/\Delta t_{\rm p}}\propto\sqrt{N_{x}}$. The total number of steps (integration steps $\times$ number of sub-steps) is therefore expected to be $N_{t}\propto N_{x}s\propto N_{x}\sqrt{N_{x}}$. This behaviour is verified in the middle panel of Figure 2 from which we conclude that STS techniques provide an effective asymptotic gain over standard explicit time-stepping proportional to the square root of the number of grid points.

Next, in the right panel of Figure 2, we plot the L1-norm errors by changing the parabolic CFL number $C_{p}=\kappa\Delta t_{\rm h}/\Delta x^{2}$ for a fixed grid resolution $N_{x}=4096$. Although comparable for $C_{p}\lesssim 1$, the errors grow linearly with the CFL number for the first order schemes like AAG-STS and RKL-1. However, for $C_{p}\gtrsim 200$, AAG-STS becomes unstable and integration is no longer possible, while RKL-1 remains stable without showing any significant undershoot even during each cycle sub-step. The second order scheme maintains approximately the same accuracy for $C_{p}\lesssim 10^{3}$ and the error starts to increase more rapidly for larger CFL although the solution remains well-behaved and positive at all times.

Finally we investigate the stability properties of the AAG-STS method alone, by varying both the parabolic CFL number and the $\nu$ parameter at the fixed grid resolution of $N_{x}=1024$ zones. Fig. 3 shows the logarithm of the error as a function of $C_{p}$ and $\nu$. The corresponding number of sub-steps $s$ is also over-plotted using white contour lines. For sufficiently small values of the damping parameter ($\nu\lesssim 10^{-4}$), we see that the error rapidly increases when $C_{p}$ exceeds $\sim 200$ ($s\sim 20$) and computations eventually become unstable for $C_{p}\gtrsim 10^{3}$ ($s\gtrsim 45$), irrespective of the value of $\nu$. For $s\lesssim 1/(2\sqrt{\nu})$, Eq. 20 shows that $C_{p}\approx s^{2}/2$ is independent of $\nu$. For larger values of the damping parameter ($\nu\gtrsim 10^{-4}$), however, the number of steps required to complete the calculation at a given CFL number increases and therefore we observe a loss of efficiency that reduces the limiting CFL number from $C_{p}\approx 10^{3}$ (at $\nu\approx 10^{-4}$) to $C_{p}\approx 200$ (for $\nu\approx 10^{-2}$). Of course, increasing the number of sub-steps at a given CFL number leads to larger stability at the cost of extra computational work. Finally, we point out that our temperature profiles are marginally affected by the numerical resolution, the effect of which is that of triggering (in case of numerical instability for $s\gtrsim$ 40) the growth of Nyquist ($k=\pi/\Delta x$) mode.

4.1.2 Saw-tooth profile with saturated heat flux

In the next example we compare the performance of the selected integration schemes by also including saturation of the heat flux. More specifically, we consider the initial sawtooth temperature profile

[TABLE]

and solve Eq. 6 by setting $\rho=1$ and $\gamma=2$ so that $\rho\epsilon=p=T$ in code units, as in the previous example. The thermal conduction flux is given by Eq. 2 with $\bm{F}_{\rm class}$ and $q$ defined in Eq. 5 and Eq. 4, respectively, with $\kappa=40$, $\phi=0.3$. Computations are carried out on the one-dimensional domain $x\in[-L_{x}/2,L_{x}/2]$ using $N_{x}=400$ zones and periodic boundary conditions. The total number of steps to reach some final time $t_{\rm stop}$ is therefore $N_{\rm step}=t_{\rm stop}/\Delta t_{\rm h}$, with $\Delta t_{\rm h}$ computed from the parabolic CFL number

[TABLE]

where we have used $t_{\rm stop}=0.1$.

We perform two sets of computations corresponding to $C_{p}=5$ and $C_{p}=50$ and compare the results with a reference solution obtained on a much finer grid. Results are shown in Fig. 4. In the case of AAG-STS we employ $\nu=10^{-3}$ (orange) and $\nu=10^{-2}$ (red).

Away from extrema, the evolution is initially dominated by the contribution of the saturated flux only (since $\partial^{2}T/\partial x^{2}=0$ everywhere except where the first derivative is discontinuous). The diffusive part of the flux acts in those region where a change in slope is present. For $C_{p}=5$, all methods yield well-behaved solutions with comparable errors whereas for $C_{p}=50$, high frequency spurious oscillations appear when using AAG-STS with a lower values of the damping parameter ($\nu=10^{-3}$) as shown by the orange line in Figure 4. Oscillations originate in proximity of the maxima where the second derivative of $T$ does not vanish and propagate downstream as the system evolves. Increasing the grid resolution or the parabolic CFL number tends to amplify this unstable behaviour. We note that oscillations disappear if saturation is not included.

On the other hand, computations remain stable when using RKL-1 (green) or RKL-2 (not shown) even for larger parabolic CFL numbers ($C_{p}\lesssim 200$). For $C_{p}\gtrsim 200$ we observe larger numerical errors in the solution even with RKL. The production of temperature oscillations at maxima is analogous to numerical oscillations produced at extrema when numerically solving the cosmic ray streaming equation (see Fig. 2.1 in Sharma et al. 2010a). Unlike in cosmic ray streaming, oscillations are not produced at temperature minima because of a smaller streaming speed ($\propto\sqrt{T}$). We have verified this dependence on streaming speed by running with a higher initial temperature (T = 30, instead of T = 10) and noticing oscillatory behavior both at temperature maxima and minima in the case of AAG-STS with $\nu=10^{-3}$, while other cases show stable solutions.

4.1.3 Ring diffusion

The 2-D Cartesian ring diffusion test problem, presented in Parrish & Stone (2005); Sharma & Hammett (2007), is useful to study the monotonicity properties of various numerical schemes for anisotropic diffusion in presence of temperature discontinuities. Temperature discontinuities are fairly common in astrophysical fluids and plasmas, and an ideal numerical scheme should not lead to negative temperatures in presence of large temperature gradients. For an explicit update with $\Delta t\leq\Delta t_{\rm p}$, Sharma & Hammett (2007) showed that temperature monotonicity is maintained if limiters are used to interpolate traverse temperature gradients (see Appendix B for details); using arithmetic mean for interpolation leads to non-monotonicity in general. The use of limiters with an implicit/semi-implicit update does not strictly maintain monotonicity but improves monotonicity substantially compared with arithmetic averaging (see Figs. 4 & 5 in Sharma & Hammett 2011).

We numerically solve the anisotropic thermal diffusion equation, Eq. 6 with $\kappa_{\parallel}=1$; saturation of heat flux is ignored for this test problem. We set $\rho=1$, $\gamma=2$, and $\rho\epsilon=p=T$ in code units. The magnetic field lines are circular with $B_{x}=y/(x^{2}+y^{2})^{1/2}$ and $B_{y}=-x/(x^{2}+y^{2})^{1/2}$, and the fluid is static. All variables except $T$ ($\rho$, $\bm{B}$, $\bm{v}$) are held fixed; $T$ evolves only because of anisotropic diffusion. The computational domain is $[-1,1]\times[-1,1]$, equal number of grid points ($N$) are used in the two directions, and periodic boundary condition is imposed on $T$. The initial condition on $T$ is

[TABLE]

where $r=(x^{2}+y^{2})^{1/2}$ and $\theta$ ($\tan\theta=y/x$) are polar coordinates. The magnetic field lines are not aligned with the Cartesian grid and the transverse heat flux (Eq. 51) is non-zero. The anisotropic diffusion equation is evolved until $t=1$ using different schemes.

Figure 5 shows the temperature at $t=1$ using RKL-2 scheme with different resolutions ($N=50,~{}100,~{}200$) and different number of STS stages ($s=20,~{}50,~{}100$); left (right) panel shows results without (with) limiters for interpolating transverse temperature gradients (see Appendix B). For all cases without limiters the minimum temperature at $t=1$ is less than 10, the initial minimum temperature.333For the present test problem with a maximum and minimum temperature of 12 and 10, the minimum temperature at later times is not negative but it will become negative if the initial minimum and maximum temperatures are say 0.1 and 10. A temperature ratio of 100 is commonplace in multiphase astrophysical flows, such as the interstellar medium. Once temperature becomes negative the MHD solver breaks down because of an imaginary sound speed. With limiters, only $s=50$ and $N=50$ run shows a temperature below 10. Therefore, monotonicity is better maintained with the use of limiters. A comparison of left and right panels shows that the use of limiters introduces large numerical diffusion perpendicular to field lines, especially for smaller resolution and larger number of stages (i.e., for smaller $N$ and larger $s$). A similar behaviour is observed for other schemes such as AAG-STS, RKC-STS and RKL-1. For very large $\Delta t_{\rm h}/\Delta t_{\rm p}$, corresponding to $s\sim\sqrt{\Delta t_{\rm h}/\Delta t_{\rm p}}\gtrsim 20$, it is better not to use the limiters because of excessive transverse diffusion. (A very large $s$, although numerically stable, also leads to a loss of accuracy.) However, limiters are necessary for preventing negative temperatures in presence of large gradients. For practical implementation we recommend the use of limiters only when temperature becomes smaller than a reasonable floor value in the anisotropic conduction step.

4.2 MHD Tests

Unlike section 4.1, the test problems described in this section solve the full set of MHD equations in presence of thermal conduction (Eqs. 1a-1d). These tests demonstrate the coupling of the MHD hyperbolic conservation laws with the parabolic update of thermal conduction.

4.2.1 Supernova blast-wave

We consider an MHD blast wave in cylindrical co-ordinates with initial parameters similar to the $L1$ model of Meyer et al. (2012) but without radiative cooling. For this 2D axisymmetric test, we solve the standard set of ideal MHD equations, taking into account anisotropy in thermal conduction flux along with saturated conduction with $\phi=0.3$ (see Eqs. 2, 3 and 4).

The initial magnetic field of $0.3~{}\mu$G is oriented along the $z-$ axis.Energy of $10^{51}\rm erg$, equivalent to a supernova explosion, is injected in a region of spherical radius ($[r^{2}+z^{2}]^{1/2}$) $7.8$ pc around the origin of a $300$ pc domain in cylindrical geometry. A high pressure is set inside this spherical region (resolved with 10 zones) such that one-third of the supernova energy goes as thermal energy. The radial velocity of this hot material is set by equating the kinetic energy to the remaining two-third of energy. The temperature in the ambient medium is set to be $8000$  K and the initial density is set to be $0.7m_{u}$ cm*-3* ($m_{u}$ is atomic mass unit) throughout the domain. Axisymmetric boundary conditions are imposed at $r=0.0$ and initial conditions are imposed at the vertical boundary $z=0.0$. The conditions at the outer boundaries in both $r-$ and $z-$ directions are imposed to be outflow.

Figure 6 shows the comparison of logarithmic value of temperature in Kelvin for runs with and without Thermal Conduction at time $t=0.1$ Myr. The RKL runs with thermal conduction use Spitzer conduction along field lines and no conduction across the field lines. The black lines shown in each panel depict the magnetic field lines. Since thermal conduction flux $\bm{F_{\rm c}}$ is only along the field lines, this anisotropy leads to an asymmetric expansion of the inner hot bubble along $z-$ axis. The evolution of temperature and magnetic field lines in Figure 6 evidently shows that the diffusion in the transverse direction is suppressed for the run with anisotropic thermal conduction. While in absence of thermal conduction the outer shock maintains its spherical structure.

Further, we compare the execution times for this test using the standard explicit time-stepping and the RKL-2 method. This comparison for varying grid resolutions is listed in table 2. The RKL-2 method is faster than the explicit method by one order of magnitude for our highest resolution run with a grid resolution of $1536^{2}$.

This observation is consistent with the fact that STS techniques like RKL-2 method provide an asymptotic gain over standard explicit time-stepping methods as demonstrated in the middle panel of Fig. 2 for the 1D Gaussian diffusion test (Sec. 4.1.1).

4.2.2 Local thermal instability

In this section we describe the results from 2-D MHD simulations of local thermal instability (TI) with anisotropic conduction (Sharma et al. 2010b). Our computational domain is a 40 kpc $\times$ 40 kpc periodic Cartesian box with a mean initial electron density of $n_{e}=0.1$ cm*-3* (mean mass per particle/electron, $\mu/\mu_{e}=0.62/1.17m_{p}$) and a uniform temperature of 0.7 keV (typical cool cluster core parameters). Only classical thermal conduction with the conductivity given by the Spitzer value (Eq. 11 in Sharma et al. 2010b) is included; i.e., $\bm{F}_{\rm c}=\bm{F}_{\rm class}$ instead of Eq. 2. We use a tabulated cooling function corresponding to the plasma of a third solar metallicity (Sutherland & Dopita 1993); the cooling function is set to zero for $T<10^{6}$ K. Following Sharma et al. (2010b), we include a spatially uniform heating rate density which is equal to the average cooling rate density over the computational domain. With this, the computational box is in global thermal balance, and this mimics the observed rough thermal balance inferred in cool core clusters. Moreover, this setup shows the exponential linear growth of the local thermal instability. With these parameters, the initial cooling time (which is approximately equal to the growth timescale of local TI; e.g., see Eq. 19c in McCourt et al. 2012) is 0.095 Gyr. A magnetic field of 5 $\mu$G is initialised at 45 degrees to the Cartesian box. Homogeneous, isotropic, isobaric random density perturbations are initialised to seed the local TI. The density perturbations ($\delta\equiv[\rho-\rho_{0}]/\rho_{0}$) are given by

[TABLE]

where $k,l$ are mode labels, $a_{k,l}=1.5\times 10^{-3}r(k^{2}+l^{2})^{-1/2}$ and $\phi_{k,l,m}=2\pi r$ ($r$ is a random number uniformly distributed between -0.5 to 0.5, which is different for the amplitude and phase and for different modes), and $L=40$ kpc is the box size. These choices give a maximum over-density amplitude ${\rm max}(\delta)\approx 0.003$. The setup is very similar to but not identical as Sharma et al. (2010b).

We run the local TI test problem using different methods for anisotropic thermal conduction: (a) fully explicit evolution in which both the MHD and conduction modules are evolved using a time step $\Delta t_{\rm p}$ (Eq. 10; this is typically shorter than the MHD CFL step $\Delta t_{\rm h}$; Eq. 9); (b) sub-cycling of conduction module in which MHD module uses $\Delta t_{\rm h}$ but the conduction module is sub-cycled and applied $\Delta t_{\rm h}/\Delta t_{\rm p}$ times using a time step of $\Delta t_{\rm p}$; (c) MHD module is evolved using $\Delta t_{\rm h}$ and conduction module is evolved using AAG-STS with $\nu=0.01$ and $s\sim(\Delta t_{\rm h}/\Delta t_{\rm p})^{1/2}$ stages (see Eq. 18); and (d) MHD module is evolved using $\Delta t_{\rm h}$ and conduction module is evolved using RKL-2 with $s\sim(\Delta t_{\rm h}/\Delta t_{\rm p})^{1/2}$ stages (see Eq. 17). The wall-time taken for different methods to run the $512\times 512$ TI test problem until 0.87 Gyr (9.16 cooling times) is listed in Table 3. All the methods use the monotonized-centered (MC) limiter to calculate the transverse terms in the anisotropic heat flux ($F_{x}^{T}$ in Eq. 51 and analogous expression for $F_{y}^{T}$). If we do not use limiters for interpolating transverse temperature gradients (and instead use simple averaging), all the different runs (a-d) blow up at some point in nonlinear evolution due to negative temperature somewhere in the computational domain. This test highlights the importance of using limiters for robustness in presence of large temperature gradients.

We note that the AAG-STS run crashes at 0.87 Gyr even with limiters because of a large $C_{p}=\Delta t_{\rm h}/2\Delta t_{\rm p}$ (see Fig. 3). Other runs could go for much longer without numerical problems. Both RKL-2 and AAG-STS methods clearly show a significant speed-up relative to sub-cycling and explicit methods. On comparing the STS methods, we see that AAG-STS is faster than RKL-2 due to a smaller number of computations per cycle. However, the first-order accurate AAG-STS scheme exhibits unstable behavior while RKL-2 maintains stability during the entire integration as described below.

Figure 7 shows temperature snapshots at different stages of TI evolution using RKL-2 (top panels) and AAG-STS (bottom panel) schemes. In the linear stage all schemes give a similar evolution. In the saturated nonlinear state AAG-STS and RKL-1 start to deviate quantitatively from each other. Unlike AAG-STS, the RKL-2 temperature snapshots at 0.86 Gyr are very similar to the ones obtained using explicit and sub-cycling methods (not shown in Fig. 7); this time is very close to the time when the AAG-STS run blows up (at 0.87 Gyr) due to numerical instability.

Figure 8 shows the time evolution of various quantities as a function of time using different methods for anisotropic thermal conduction. Top left panel of Figure 8 shows that the evolution of kinetic energy in the box is similar for all the runs. The top right panel shows the average magnetic energy evolution. Here, AAG-STS deviates from the other runs in the non-linear stage. Similar deviations are seen in AAG-STS for $T_{\rm min}$ and $T_{\rm max}$ (minimum and maximum temperature in the computational domain), but not for RKL-2. The numerical fragility of AAG-STS due to an imperfect stability parameter $\nu$ is evident from spikes in $T_{\rm min}$ and $T_{\rm max}$. The AAG-STS run blows up at 0.87 Gyr due to the numerical instability of AAG-STS for a large number of sub-stages.

5 Parallel scaling of RKL conduction

In this section, we demonstrate the parallel scaling of isotropic and anisotropic thermal conduction in PLUTO code. For this purpose, we use the blast wave test problem in three dimensions on a Cartesian grid with size $L=30$ pc. To set up a blast wave we initialise an ambient static medium with density $\rho_{0}$ and pressure $p_{0}$. A blast region with high density and pressure is set within a spherical radius $r_{0}$ around the origin. The density and pressure in this region is a factor 10 and 1000 times larger than the ambient medium respectively. We smooth the transition between the blast region and the ambient medium using a smoothing function as follows

[TABLE]

where, $\Lambda[r]=(1+\exp([r-r_{0}]/[0.1~{}{\rm pc}]))^{-1}$ is a smoothing function, $r=(x^{2}+y^{2}+z^{2})^{1/2}$ is the spherical radius, $p_{0}=2.295\times 10^{-9}$ dyn cm*-2*, $\rho_{0}=2.1\times 10^{-23}$ g cm*-3*, $r_{0}=1.0$ pc, and $\gamma=5/3$. Additionally, for MHD runs the initial magnetic field is along $y$ axis, $\bm{B}=B_{0}\hat{\bm{y}}$ where $B_{0}=\sqrt{2p/\beta}$ with $\beta=1$. Thus, the field strength inside the high pressure blast region is about 30 times larger than the ambient medium. We use Spitzer conduction along field lines and no conduction across the field lines for the MHD run. The HD run uses isotropic Spitzer conduction. For comparison, we also include a HD run without conduction.

We carried out strong scaling studies using a $512^{3}$ Cartesian grid on CPU-only nodes of IISc Cray XC40 cluster SahasraT444http://www.serc.iisc.in/facilities/cray-xc40-named-as-sahasrat/. Unlike the 2-D blast wave test problem in section 4.2.1, we choose a 3-D blast wave problem for testing parallel strong scaling on more than ten thousand processors. In strong scaling, the problem size remains fixed but the number of processors is increased progressively. A large enough problem size is needed (that is why 3-D instead of 2-D test problem is chosen) for communication not to dominate over computation even on largest number of processors.

We run the blast wave test problem for $t=0.5$ kyr on processors ranging from 256 to 16384. On average, the number of sub-steps $s$ ranges between 5-12 for the HD run and between 3-8 for MHD.

The left panel of Figure 9 shows the wall clock time as a function of processor count for our strong scaling studies. The wall clock time behaves as expected, with the MHD-RKL run taking longer than the HD-RKL run. Both runs with conduction show close to the inverse linear scaling with the processor count. The right panel of Figure 9 shows that the speed-up improves with increasing core count, but drops to 80% for more than $10^{4}$ processors. The same trend is also seen for pure HD run, indicating that the addition of thermal conduction has not resulted in any performance degradation. Our strong scaling tests show that RKL methods yield scalable algorithm on PetaFlop facilities.

This feature is usually not shared by implicit methods which, as already mentioned in Sec. 3.2, require inverting large sparse matrices, an operation difficult to achieve efficiently on multi-core systems (Botchev & van der Vorst, 2001). Nevertheless, the debate on which approach can be more efficient in massive parallel computations has yet to be settled. In the work by Pakmor et al. (2016), for instance, a semi-implicit scheme is employed to solve the cosmic-ray transport equation (coupled to the MHD equations) on unstructured moving mesh in the context of galaxy dynamics. Their semi-implicit solver requires solving a single linear system of equations per timestep and weak scaling tests claim good parallel efficiency, up to 480 cores. However, more complex scenarios including nonlinear systems of equations (e.g Ohmic diffusion, ambipolar diffusion, etc.) may introduce additional complexities making the applicability of implicit scheme less efficient. Although a comparison with implicit methods is outside the scope of this paper, we point out to the recent work by Caplan et al. (2017), where explicit STS schemes are compared with implicit Krylov solvers in the context of magnetised solar corona. The left panels of Fig. 11 in their paper shows the superiority of explicit schemes over implicit ones for their choice of problem and model. The comparison also shows that the strong parallel scaling of implicit methods saturates and deviates sharply from ideal behavior at large number of cores ($>1000$).

The explicit RKL schemes presented in this work confirm this prediction (the scaling performance with explicit RKL method shown in Fig. 9 is able to achieve a high efficiency of 80% even up to 1000 cores) and can be naturally extended to more complex systems of equations with with minimal modifications.

6 Conclusions

In this paper we have discussed various numerical methods implementing anisotropic thermal conduction coupled with the standard set of MHD equations. In particular, we have described the second order (in time) accurate Runge-Kutta Legendre (RKL-2) super-time stepping (STS) method implemented in PLUTO code. We have then compared these numerical methods on simple test problems, and also on astrophysical test problems like blast wave and thermal instability in which MHD evolution is coupled with conduction. The major conclusions of our paper are:

Using the 1-D Gaussian diffusion test, we show that the RKL-2 scheme is second order accurate in time, in comparison to the standard AAG-STS (Alexiades et al. 1996) scheme which is first order. 2. 2.

Super-time stepping schemes based on Chebyshev polynomials such as AAG-STS require an ad-hoc damping parameter ($\nu$), which has to be big enough for numerical stability. Figure 3 shows that AAG-STS becomes unstable for number of substages $s\gtrsim$ few 10s. Moreover, for $2s\sqrt{\nu}\gtrsim 1$ there is no real speed-up using AAG-STS (see Eq. 20). In absence of such a parameter and with sufficient inherent damping, RKL-STS schemes are more robust and can better exploit the super-time stepping strategy (see the right panel of Fig. 2). The robustness of RKL-STS is also useful for implementing saturated conduction (e.g., see Fig. 4). 3. 3.

The Cartesian ring diffusion test problem (Fig. 5) shows that all STS schemes break down if we use a large number of stages with low resolution. Moreover, the use of limiters to interpolate the transverse temperature gradient leads to somewhat larger perpendicular diffusion. Limiters do help prevent non-monotonicity of temperature in presence of temperature discontinuities. For practical implementation of anisotropic diffusion, limiters need not be used as a default option, but may be recommended at locations and times at which the temperature falls below a floor value. The local thermal instability test problem in section 4.2.2 demonstrates the utility of limiters in presence of temperature discontinuities. 4. 4.

As demonstrated in Tables 1, 2 & 3, and also indicated by the middle panel of Figure 2, the run-time with super-time stepping is substantially shorter than with explicit update or even sub-cycling. 5. 5.

Last, but perhaps most importantly as shown in Figure 9, the explicit schemes such as RKL-STS shows an excellent scaling (efficiency of $\sim 80\%$ up to 104 processors) on modern distributed PetaScale supercomputers.

Acknowledgments

The authors would sincerely like to thank the referee for the valuable comments which have played a significant role in improving the paper. BV would like to thank the support provided by University of Torino and also would like to acknowledge the hospitality of IISc during the research visit in 2016. This work is partly supported by the DST-India grant no. Sr/S2/HEP-048/2012 and an India-Israel joint research grant (6-10/2014[IC]). We thank the SERC-IISc staff for facilitating our use of SahasraT cluster for the parallel scaling runs; only with their intervention we could finish our runs in a reasonable time. PS acknowledges the hospitality of KITP where this paper was completed. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915.

Appendix A Illustration of RKL-2 scheme

This Appendix is closely based on Meyer et al. (2012, 2014) and is included here for completeness. Second order accuracy for the RKL-2 scheme can be achieved by matching the first three terms in Eqs. (14) & (15); i.e., by imposing $R_{s}(0)$ = 1, $R_{s}^{{}^{\prime}}(0)$ = 1 and $R_{s}^{{}^{\prime\prime}}(0)$ = 1. These three requirements can be used to estimate the values of three coefficients in Eq.(15) as,

[TABLE]

For $s<2$, we can choose $b_{0}=b_{1}=b_{2}=1/3$ (Meyer et al., 2014). Here we have used various properties of the Legendre polynomials and their derivatives, i.e.,

[TABLE]

Like Eq. (15) for the stability polynomial at the end of $s$ stages, the stability polynomial till $j$ stages is chosen to be $a_{j}+b_{j}\mathcal{P}_{j}\left(1+w_{1}(\tau\lambda)\right)$.

The RKL-2 scheme with $s$ stages can be expressed as,

[TABLE]

where the coefficients, $\mu_{j}$, $\tilde{\mu}_{1}$, $\tilde{\mu}_{j}$, $\nu_{j}$ and $\tilde{\gamma}_{j}$ can be obtained using the recursion relation for the Legendre polynomials (see Eq.(16)) and rearrangement of terms. The expressions for these coefficients are given by,

[TABLE]

To get a better sense of the scheme, we explicitly list the value of the above coefficients for a representative small value of $s=3$. For this three stage scheme we have,

[TABLE]

On substituting the above values of various coefficients for $s=3$ in Eq.(A) we get;

[TABLE]

where $\lambda$ is the eigenvalue of parabolic operator matrix $\mathsf{M}$. Note that the first three terms of final solution of $Y_{3}$ (Eq. A) match the Taylor expansion of the exponential function (Eq.(14). Thus this scheme is second order-accurate with a leading order in error given by the fourth term, i.e., $(\tau\lambda)^{3}/15$.

Finally to estimate the value of the super-step $\tau$, we should have $|R_{3}(\tau\lambda)|\leq 1$. The stability condition will be satisfied only if the argument of LP is $\geq-1$ (note that $\lambda$ is non-positive), i.e., $-5\leq\tau\lambda\leq 0$. While the similar stability condition for one step Euler-scheme requires $-2\leq\tau\lambda\leq 0$. This clearly shows the $s=3$ stage RKL-2 method allows to choose a large time-step as compared to the Euler scheme. Higher values of sub-steps $s$, will result in a wider range of stable time steps giving an obvious advantage over standard explicit schemes.

Appendix B Limiting transverse temperature gradients for monotonicity

Simple finite differencing of Eq. 6 accentuates temperature extrema in presence of large temperature gradients expected naturally in astrophysical plasmas (e.g., Sharma & Hammett 2007). The heat flux can be decomposed into normal and transverse components (Eqs. 50, 51). Expressing the heat flux , $\bm{F}_{\rm c}$ with components normal and transverse to the field lines we obtain the following explicit formulation of Eq. 6 in 2D Cartesian geometry (generalisation to 3-D and non-Cartesian coordinates is straightforward),

[TABLE]

For simplicity we consider the classical limit (i.e., ignore saturated conduction). On comparing Eq. 49 with Eq. 3, we get

[TABLE]

as the normal and transverse components of the heat flux (analogous expressions can be constructed for $F_{\rm c,y}^{N}$ and $F_{\rm c,y}^{T}$).

The normal component of the heat flux is naturally face-centred and always carries heat from higher to lower temperatures. However, the transverse temperature gradients are not naturally located at the faces, and need to be interpolated there. As the transverse flux can have any sign (Eq. 51), without special treatment, this component can lead to negative temperatures in regions with large temperature gradients (this also applies to saturated conduction in Eqs. 2 & 4; therefore, the transverse temperature gradient required to evaluate ${\rm sgn}[\bm{\hat{b}}\cdot\nabla T$] for saturated flux should also use limiters for robustness). Sharma & Hammett (2007) introduced limiters (similar to those used in the reconstruction step in finite volume methods; LeVeque 2002) to interpolate the transverse temperature gradients at the cell faces and showed that the resulting explicit scheme preserves temperature extrema.