A Fast Semi-implicit Method for Anisotropic Diffusion

A Fast Semi-implicit Method for Anisotropic Diffusion

Prateek Sharma111Chandra Fellow psharma@astro.berkeley.edu Gregory W. Hammett Theoretical Astrophysics Center and Astronomy Department, University of California, Berkeley, CA 94720 Princeton Plasma Physics Laboratory, Princeton, NJ 08543

Simple finite differencing of the anisotropic diffusion equation, where diffusion is only along a given direction, does not ensure that the numerically calculated heat fluxes are in the correct direction. This can lead to negative temperatures for the anisotropic thermal diffusion equation. In a previous paper we proposed a monotonicity-preserving explicit method which uses limiters (analogous to those used in the solution of hyperbolic equations) to interpolate the temperature gradients at cell faces. However, being explicit, this method was limited by a restrictive Courant-Friedrichs-Lewy (CFL) stability timestep. Here we propose a fast, conservative, directionally-split, semi-implicit method which is second order accurate in space, is stable for large timesteps, and is easy to implement in parallel. Although not strictly monotonicity-preserving, our method gives only small amplitude temperature oscillations at large temperature gradients, and the oscillations are damped in time. With numerical experiments we show that our semi-implicit method can achieve large speed-ups compared to the explicit method, without seriously violating the monotonicity constraint. This method can also be applied to isotropic diffusion, both on regular and distorted meshes.

implicit methods, finite differencing, monotonicity, anisotropic diffusion

1 Introduction

Anisotropic diffusion equation arises frequently in diverse applications: microscopic transport in magnetized plasmas bra65 (); image processing per90 (); diffusion-tensor magnetic resonance imaging bas02 (); thermal properties of crystals dia91 (); transport in geological formations ber02 (), etc. In sha07 () we showed that simple finite-differencing of the anisotropic diffusion equation resulted in unphysical numerical heat fluxes, which lead to negative temperatures at large temperature gradients. Negative temperatures, in addition to being unphysical, result in an imaginary sound speed and associated numerical instabilities.

Anisotropic diffusion equation satisfies important mathematical properties such as monotonicity preservation along the direction of diffusion, so that no new extrema are created in that direction and any existing extrema are not accentuated, i.e., maxima must drop or be unchanged, and minima must increase or be unchanged (e.g. see ker81 (); nor07 () and references therein). For simplicity we will refer to this as an extrema reducing property. In two or three dimensions the anisotropic diffusion equation of the form we are considering here can always be transformed to the form , where is a coordinate that varies along field lines and are field-line label coordinates that are constant along a magnetic field line, so the monotonicity properties of 1-D diffusion should be preserved.

In sha07 () we proved that the use of slope limiters (e.g., see lev02 ()) to interpolate temperature gradients at cell faces guarantees that temperature extrema are not accentuated, as required physically. The extrema-reducing property ensures that the temperature is positive for a CFL stable timestep. Since the temperature is positive, numerical instabilities that plague simple finite differencing of anisotropic diffusion (because of an imaginary sound speed!), do not arise with the use of limiters. Because of this desirable property our method has been used in astrophysical magnetohydrodynamic (MHD) simulations where thermal conduction is anisotropic, and large temperature gradients can arise sha10 (); par09 (); rus10 (). In addition to thermal conduction, limiters have proved useful for anisotropic viscosity with large gradients don09 (), and may help in various problems with anisotropic transport; e.g., cosmic ray streaming sha09 (). There is another explicit method that has been applied to astrophysical problems and avoids the problem of negative temperature in presence of large temperature gradients, though it is non-conservative and somewhat more complex ras08 (). We also discuss some other recent work below.

A key limitation of all explicit methods is that the timestep is limited by the usual CFL condition, , where is the grid spacing and is the anisotropic diffusion coefficient. For some applications, this timestep constraint is rather severe and the conduction timestep can be much smaller than the MHD CFL time limit. In such cases an implicit method, where there is no stability limit on the diffusion timestep, is desirable. Although it is straightforward to difference a linear anisotropic diffusion equation implicitly, the resulting scheme is still not monotonicity-preserving. E.g., see Table 1 in bal08 (), which shows that temperature oscillations remain till late times even with an implicit method, just as with the explicit schemes without limiters. One can try to solve the anisotropic diffusion equation fully implicitly, using limiters to prevent temperature oscillations, but the nonlinearities in the limiters will require a careful iterative treatment (some studies have found that these kinds of nonlinear limiters make iterative solvers more difficult).

We have experimented with a Jacobian-free nonlinear iterative implicit method (a two-stage Richardson iteration extension of the LGMRES(1,1) version of Loose GMRES Baker05 (), a variant of the Generalized Minimal RESidual methodsaa86 () with restarting) but found that it requires a fairly large number of iterations per time step, because the Jacobian matrix is not strongly diagonally dominant for large timesteps and has a large condition number.222Even for a linear 2-D Poisson problem on an grid, Loose GMRES or conjugate-gradient methods require iterations by themselves, or iterations if combined with a sufficiently good preconditioner like Modified ILUGustafsson78 (). On the other hand, the splitting method employed in this paper uses tri-dagonal implicit solvers that are equivalent to only iterations and so are quite fast by comparison. The fast method proposed here might be able to serve as an effective preconditioner to further accelerate an unsplit iterative method.

We have also experimented with an explicit method which is stable for timesteps longer than the CFL limit, and where the internal iteration time-steps are chosen based on the properties of Chebyshev polynomials ale96 (). We were not able to obtain a speed-up of more than 10 compared to the CFL-limited scheme, irrespective of resolution, for any of the parameters that we varied. Moreover, the parameters for which maximum speed-up is obtained, without becoming numerically unstable, are difficult to choose (this is true even for isotropic diffusion!).

Here we present a conservative, directionally-split, semi-implicit method which is numerically stable for any choice of timestep, and is easy to implement. The method is based on directional splitting where the heat fluxes in each direction are updated sequentially. The heat flux in each direction, e.g., (see Eq. 2), consists of two terms: , the ‘normal’ term where temperature gradient need not be interpolated, and , the ‘transverse’ term which involves temperature interpolation with limiters. In our directionally-split method the ‘normal’ terms are treated implicitly, and the transverse terms are treated explicitly. The directional splitting of ‘normal’ implicit terms results in a tridiagonal matrix which can be solved very quickly. The explicit treatment of ‘transverse’ terms with limiters ensures that extrema are not accentuated. The resulting scheme, while not strictly monotonicity-preserving for large timesteps, results in only small amplitude temperature oscillations which are damped in time. Speed-up of order 100-1000, compared to the explicit scheme, is easily achieved for our test problems.

There has been some interesting recent work on another approach to the problem of preserving positivity in presence of anisotropic diffusion tensors (or diffusion on distorted meshes), based on expressing the flux at cell faces in terms of the advected quantity at the cell centers on either side of the face (a “two-point flux expression”), but where the coefficients of this flux depend on the transverse gradients and so is nonlinear (for example see Lipnikov10 (); Sheng09 ()). Future work could compare the nonlinear limiters and implicit solvers that are used in our algorithm with these other algorithms on the types of test problems considered here and in sha07 ().

The paper is organized as follows: Section 2 presents the method in detail and shows that the scheme is linearly stable for large timesteps. Section 3 presents results from three test problems which show the practical utility of our method. We conclude and discuss applications of our method in Section 4.

2 The Method

The anisotropic diffusion equation in its simplest form is given by


where is the temperature, is the heat flux along magnetic field lines, is the thermal diffusion coefficient (with dimensions ), and is the magnetic field unit vector. We will assume as a given function of space and time, and as being a constant for simplicity (in sha07 () we showed that a harmonic average should be used for interpolation of the diffusivity for numerical stability). We assume two-dimensions with a uniform Cartesian grid and a constant diffusion coefficient. The generalization to a nonlinear diffusivity, general coordinate system, and three dimensions is straightforward.

The semi-implicit method is obtained by directional splitting of the heat flux updates in each direction, given by


where magnetic field unit vectors are interpolated at the appropriate cell faces (simple averaging is fine for magnetic field unit vectors), and


are the temperature differences centered at appropriate faces, and stands for a limiter. See sha07 () for a discussion of limiters in this context; here we will use the slope limiter of van Leer van77 (); lev02 (); is symmetric in its arguments, where


We have experimented with other slope limiters (e.g., minmod and monotonized central [MC] limiters). We observed that monotonicity properties are much better with diffusive limiters (such as minmod and van-Leer) as compared to sharper limiters such as the MC limiter (see lev02 () for properties of different limiters) for our semi-implicit scheme with large timesteps. However, more diffusive limiters result in larger perpendicular diffusion.

Defining the components of the diffusion operator on the right hand side of Eq. (1) as


(there is no implied summation on the right hand side of this definition), then the method in Eqs. (3) & (4) can be expressed as


Our formulation treats the ‘normal’ temperature derivative terms, which are guaranteed to result in a heat flux in the correct direction, implicitly. The ‘transverse’ terms are treated explicitly, and employ limiters that ensure that the temperature extrema are not accentuated. Directional splitting results in a quickly solvable tridiagonal matrix for each directional update. Instead of updating followed by we can also update the heat fluxes in the reverse order. Results do not depend substantially on the order of updates. Since our split operators are individually only first order accurate in time, Strang-splitting (e.g., see lev02 ()) will not improve the accuracy of our method (and high accuracy is not a priority for components of the solution that are strongly damped anyway). Moreover, Strang-splitting applied to Eqs. (3) and (4) is numerically unstable for large timesteps.

2.1 Linear Stability Analysis

One can perform the von Neumann linear stability analysis on Eqs. (3) and (4). Let us assume a single temperature mode , where is the amplification factor in time, and , are the wavenumbers in the and directions. The amplification factor can be written as , where and are amplification factors for the substages in Eqs. (3) and (4), respectively. On substituting the discretized temperature eigenmode in Eq. (3), and using trigonometric identities, one obtains


where, for simplicity, we have assumed that , are constant in space. We use arithmetic averaging instead of the nonlinear limited averaging for the transverse temperature gradient. Similarly, for update in the direction,


Now the amplification factor after a full timestep is


where and . From Eq. (12) we get

which is guaranteed to be since for all real , . Thus our scheme (Eqs. 3 & 4) is unconditionally stable like the usual implicit methods (see Fig. 1). The unconditional stability of our scheme also holds for the case of a general symmetric, positive diffusion tensor

with non-negative eigenvalues (i.e., ). Although large timesteps are stable, we cannot use very large timesteps because of loss of accuracy.

Figure 1: Contour plots of the amplification factor (see Eqs. 10, 11) for and different Courant factors (ncfl, see Eq. 16). Also shown is the analytic amplification factor (, assuming , where ) for ncfl=10; there is no damping along because temperature gradient (which is along ) is perpendicular to .

Fig. 1 shows the amplification factor () for a timestep longer than the Courant time step (characterized by ncfl; see Eq. 16) for a fixed magnetic field unit vector . As expected, the amplification factor is close to the analytic expectation (also shown in the same figure) for smaller ncfl. The damping rate is not sufficiently large (compared to the analytic solution) for unit vectors corresponding to the first and the third quadrant in -space. The maximum growth rate for the points in the first and third quadrants occurs at small scales (compared to the box size), but this scale and the amplification factor becomes larger for a larger ncfl. Similar low damping rate arises for the Crank-Nicolson method for isotropic diffusion for large time steps. The slow damping of modes parallel to at small scales can be avoided if we modify our scheme to a 4-step scheme, with the implicit terms in Eqs. (3) & (4) only applied for and two extra fully-implicit steps applied for . Schematically this 4-step scheme is given by (see Eq. 8 for notation)


However, for the test problems discussed in §3 our 2-step method behaves satisfactorily even for ncfl as large as 1000. And moreover, perpendicular diffusion, non-monotonicity, and computational cost are slightly worse for the 4-step method (one can try different orderings of the operators in Eq. 13 but this, and its analog where and are interchanged, worked best for our test problems). Thus we do not discuss the 4-step method in detail. Another feature in Fig. 1, which is seen for all ncfl, is that the smallest scale modes in the direction perpendicular to field lines are damped; this corresponds to cross-field diffusion at small scales.

Expanding the amplification factor (; Eqs. 10 & 11) in the limit , (i.e., large length scales), one gets


which should be compared to the analytic amplification factor , where . In the limit , the difference between the numerical and analytical amplification factors is ; i.e., this method is first order in time.

It is instructive to do a similar stability analysis in three dimensions, with , , as the amplification factor in each directional update. In this case,

and analogous expressions are obtained for other directions. It is easy to show that the absolute value of the amplification factor is not guaranteed to be for a large . Thus, our scheme is not unconditionally stable in three dimensions. By numerically evaluating for different parameters333We calculate the maximum of on a grid with resolution upto in , for a given . Then we try to find the maximum of (assuming ) for which . we have verified that the stability condition for our scheme in three dimensions is (assuming )

The corresponding stability limit for the explicit scheme in three dimensions is . Thus, our scheme can attain a speed-up of relative to the explicit method. Although our 2-step scheme (Eq. 9) is not unconditionally stable in three dimensions, we have numerically verified that the 4-step method (Eq. 13) is unconditionally stable in three dimensions.

We have also experimented with a variant of the alternate direction implicit (ADI) schemes for Eq. (1), inspired by their application to isotropic diffusion pre92 (). Specifically, we tried (see Eq. 8 for notation)

However, presence of the transverse terms () makes the scheme unstable for timesteps larger than a few times the CFL timestep, unlike in the case of isotropic diffusion where it is unconditionally stable. Even for isotropic diffusion, ADI does not give strong damping in the large time step limit, i.e., it is A-stable but not L-stable. Therefore, we do not consider ADI further.

Notice that our fully implicit scheme is only first order accurate in time, but for dissipative processes this is often adequate, since one is most interested in well-resolved components of the solution which are only weakly damped in a single time step. Another variant of Eqs. (3) and (4), where the explicit ‘transverse’ () terms are symmetrized with respect to the and updates (we were trying this to see if this scheme has better monotonicity properties as compared to our method), is numerically unstable for large s even though the linear stability analysis predicts an unconditional stability. Thus, linear stability is only a necessary (and not sufficient) condition for numerical stability, especially since the limiters are nonlinear.

3 Numerical Tests

In this section we describe various tests for our semi-implicit scheme (Eqs. 3 and 4).

3.1 Diffusion in a ring

The ring diffusion test (see par05 (); sha07 ()) involves the diffusion of a hot patch in fixed circular magnetic field lines. This is a crucial test to check monotonicity properties of the anisotropic diffusion scheme because field lines make all possible angles with respect to the Cartesian grid. At late times (a few diffusion times across the ring), the temperature is expected to be uniform along each magnetic field line in the ring. The computational domain is a Cartesian box. The initial temperature distribution is


where and , and , . Reflective boundary condition ( at boundaries in the direction, at boundaries) is used for temperature; magnetic field and conduction vanishes outside . The parallel conduction coefficient ; there is no explicit perpendicular diffusion.

Figure 2: Temperature contour plots at for the ring diffusion test problem using a grid. Our semi-implicit method (Eqs. 3 and 4) is used with different CFL numbers (ncfl). Also shown is the temperature plot with the fully explicit method (using the van Leer limiter) for comparison.

We quantify the timestep by ncfl, where


Since our scheme is unconditionally stable in two dimensions, we experiment with the CFL number (ncfl) to quantify the speed-up relative to the explicit method that we can obtain, without degrading the solution. Fig. 2 shows the temperature contour plots at using different ncfl for a box. As expected from section 2.1, our scheme is numerically stable even for ncfl 1. However, the solution deteriorates for an extremely large ncfl; e.g., the temperature profile for ncfl=10000 looks quite different from rest of the others as there is considerable numerical diffusion out of the circular ring. This figure shows that large speed-ups ( 1000 in the case of Fig. 2) are possible as compared to the explicit method. Moreover, temperature oscillations at extrema are not as severe as schemes without limiters.

Extrema are not accentuated for our semi-implicit method because the transverse term, responsible for non-monotonicity, vanishes at temperature extrema because of limited averaging. However, it is not guaranteed that the temperature will be bound by the initial temperature extrema; this is because temperature oscillations may arise at non-extremal locations for a large CFL factor (ncfl). These newly created extrema will not be accentuated though. We never encountered temperature oscillations for the fully-explicit method with limiters which uses a CFL-stable timestep.

Figure 3: Temperature contour plots at for different grid resolutions but a fixed CFL number (ncfl=1000), using our semi-implicit method.

Fig. 3 shows the temperature profiles at for different grid resolutions but with a fixed ncfl=1000. The figure shows that the temperature profiles are very similar for . The maximum speedup relative to the explicit method (i.e., maximum value of ncfl), without seriously affecting the solution, is achieved for the highest resolution simulations (where significant speed-up is in-fact desired), as seen from Figs. 2 & 3. While ncfl=1000 can be used for , the solution for a lower grid resolution with ncfl=1000 is quite diffusive in the perpendicular direction; also parallel diffusion seems to be suppressed (as is the case for ncfl=10000 and in Fig. 2).

Figure 4: Minimum temperature over the computational domain as a function of time (all timesteps are plotted) for our semi-implicit method using different Courant factors (ncfl). As in Fig. 2, grid resolution is . For comparison, the minimum temperature for the explicit method without limiters and with ncfl=1 is -0.41.

As mentioned earlier, our semi-implicit scheme is not guaranteed to be monotonicity-preserving. However, the oscillations originating at large temperature gradients are of small amplitude and are damped away quickly in time (see Fig. 4). In contrast, temperature oscillations for the explicit methods without limiters are large and persist till late times (see Fig. 7 in sha07 ()). Of course, the amplitude of temperature oscillations is proportional to the ratio of maximum to minimum temperature at the discontinuity, but Fig. 4 shows that the temperature is still maintained positive for ncfl as large as 10000! The temperature ratio of 100 (as in our test problem) is similar to the temperature range found in practical applications such as the transition of the solar chromosphere at K to the coronal temperature of K. Notice that the minimum temperature respects monotonicity constraint for the first timestep because all points in the initial condition are extrema and the transverse term vanishes. The temperature oscillations start from non-extremal points at later times.

Figure 5: Minimum temperature over the computational domain as a function of time (all timesteps are plotted) for different resolutions using our semi-implicit method. As in Fig. 3, the Courant factor (ncfl) is fixed to be 1000.

Fig. 5 shows the minimum temperature over the domain as a function of time for different grid resolutions and a fixed ncfl=1000. While Fig. 4 shows a clear trend of increased non-monotonicity as ncfl is increased, there is no systematic variation in the magnitude of temperature oscillations with the grid resolution for a fixed ncfl. This is expected because the factor is of the same order (ncfl) for all resolutions. However, Fig. 3 clearly shows that for a fixed ncfl a more accurate temperature profile is obtained for a higher grid resolution.

For a realistic problem ncfl should be , where is the scale on which we want temperature to be calculated accurately. Notice, that this factor increases with the grid resolution for a fixed , so higher resolution runs will more accurate for a fixed ncfl (see Fig. 3). Another constraint on ncfl comes from the positivity requirement. The magnitude of non-monotonicity is roughly independent of the grid resolution for a fixed ncfl (see Fig. 5), but depends on the ratio of maximum to minimum temperature at the discontinuity. E.g., for our test problem the initial temperature ratio is 100 and the maximum relative non-monotonicity (defined as , where is the initial minimum temperature and is the minimum temperature of all times) for ncfl=1000 (see Fig. 5) is . We have numerically verified that the relative non-monotonicity scales with the maximum to minimum temperature ratio (and of course non-monotonicity is larger for a larger ncfl). As mentioned before, non-monotonicity is worse for steeper limiters such as the MC limiter, but is less severe for diffusive limiters such as minmod.

3.2 Convergence & measuring

We perform a test problem with a smooth solution, described in sov04 (), to measure perpendicular numerical diffusion as a function of grid resolution and the Courant factor (ncfl). A two-dimensional Cartesian box ([-0.5,0.5][-0.5,0.5]) is initialized with a zero temperature. Temperature is fixed to be zero at the domain boundaries at all times. We solve the anisotropic diffusion equation (Eq. 1) with a source term


where . The fixed magnetic field is generated by a flux function , so that the magnetic field unit vectors are along the contours of constant . And since temperature is driven by the source term, temperature is always constant along field lines. If there is no diffusion across field lines, temperature should rise with time. However, because of finite numerical diffusion in the perpendicular direction, it reaches a steady state. The steady state solution for the temperature, if we assume a finite perpendicular diffusivity , is , independent of . We use the asymptotic value (in time) of the maximum temperature to calculate . This test is slightly modified from sov04 () in that we do not include an explicit perpendicular diffusivity; this is because for the problems of our interest perpendicular conduction is negligible.

Figure 6: The ratio of perpendicular (numerical) to parallel diffusivity as a function of grid size (triangles) for the smooth test problem in sov04 (). The Courant factor (ncfl) is fixed to be 1000. Dotted line shows a second order convergence.

Fig. 6 shows the ratio of the perpendicular numerical diffusivity and the parallel diffusivity () as a function of grid resolution for a fixed ncfl=1000. Perpendicular diffusion, which scales with , shows close to a second order convergence with the grid resolution, as expected. Numerical diffusion is not sensitive to the Courant factor (ncfl); e.g., we verified that is roughly independent of ncfl up to ncfl=10000 for . Close to second order convergence of perpendicular numerical diffusion (and the independence of ncfl) was also seen for the ring diffusion test problem in section 3.1.

3.3 Thermal Instability

Figure 7: temperature (keV) at Gyr ( 10 cooling times in the initial state) using our semi-implicit method (left) and the explicit method (right). The timestep for the conduction module in the semi-implicit scheme is chosen to be equal to the CFL timestep for the rest of the code () so that ncfl=, where , and is the maximum thermal diffusivity over the whole box. The conduction module is subcycled times for the explicit method.

We also tested our method for a realistic astrophysical application, namely thermal instability in the intracluster medium, the X-ray emitting hot plasma pervading the massive clusters of galaxies. For astrophysical motivation and details about numerical set up see sha10 (). We perform two-dimensional MHD simulations with anisotropic thermal conduction, using a periodic Cartesian box (40 kpc 40 kpc) with 1024 grid points in each direction. The initial temperature is 0.78 keV and initial electron number density is 0.1 cm. Identical pattern of small amplitude density/temperature perturbations are initialized (such that the pressure is uniform) to seed the thermal instability. The functional form of heating and cooling is such that cooling increases faster than heating for the cooler plasma, and vice versa for the hotter plasma. Thus, heating/cooling runs away and the plasma segregates into a two-phase medium. Net heating averaged over the whole box equals net cooling, so that the total thermal energy content of the box does not change with time. Thermal instability is in the isobaric limit; i.e., cooling time sound crossing time over all scales. Thermal conduction is primarily along field lines aligned initially at to the box; field lines roughly maintain their geometry even in the nonlinear stage. Small diffusion perpendicular to field lines is added for numerical convergence (see sha10 () for details). Cold filaments aligned along the direction of the local magnetic field arise nonlinearly (see Fig. 7) because thermal conduction along field lines suppresses growth of small scale modes.

Fig. 7 shows the temperature in the nonlinear state of thermal instability obtained by treating thermal conduction using our semi-implicit scheme (left) and using the explicit scheme (right) with the van Leer limiter. The temperature plots are almost identical, establishing the practical utility of our method. The conduction timestep (, where is the maximum thermal diffusivity over the whole box and is the grid size) in the initial state is 3 times the CFL timestep limit of rest of the code (). As thermal instability becomes nonlinear and the hottest plasma becomes hotter, decreases rapidly relative to , because of the sensitive dependence of conductivity on temperature (; see bra65 ()). The temperature-dependent conductivity is interpolated at the faces using the current temperature (see sha07 () for details about interpolation of conductivity). At 0.95 Gyr (the time corresponding to Fig. 7) ; thus the explicit method is subcycling the conduction module for 100 times, whereas the conduction module is applied only once for our semi-implicit scheme. Situation become worse with time because the hottest plasma in the box becomes hotter in time! Thus we are able to run much faster with our semi-implicit scheme, without affecting the solution and without violating temperature positivity. This example demonstrates the practical utility of our method. Also, recall that the stability limit for explicit diffusion scales as compared to for the hyperbolic terms, so our scheme will be even more useful at higher resolution and with mesh refinement.

4 Conclusions

We present a simple, directionally-split, conservative, semi-implicit method for anisotropic diffusion which is linearly stable for large timesteps. Directional splitting results in a tridiagonal matrix equation for each direction, which can be solved exactly and efficiently. For problem on a grid our scheme (Eqs. 3 & 4) requires two independent tridiagonal solves. In comparison, the fastest unsplit methods like LGMRES will require iterations to converge! Similarly, compared to the explicit method our scheme is ncfl ( 10-1000) times faster. Our method should be easily implemented in parallel using standard parallel linear algebra libraries like ScaLAPACK bla97 ().

Although our method is not monotonicity-preserving for arbitrarily large timesteps, the temperature oscillations are of quite small amplitude and are damped with time. Using test problems we show that large speedups (up to 100-1000 for our test problems) are achieved compared to the explicit method, without seriously violating the monotonicity constraint. A similar directional splitting may also prove effective for isotropic diffusion. Although ADI is numerically stable, fast, and second order accurate in time for isotropic diffusion, it does not give strong damping in the large timestep limit (just like the Crank Nicolson scheme; e.g., pre92 ()).

We also tried unsplit methods, both fully implicit using limiters for the ‘transverse’ terms, and semi-implicit where only ‘transverse’ terms with limiters are treated explicitly. The limiters lead to nonlinearities that require some care with an iterative solver, and these unsplit methods result in a large sparse matrix equation which is much more expensive to solve than the tri-diagonal systems of the split method, even with an iterative solver like conjugate gradients or Loose GMRES. Although both of these unsplit methods with limiters appear to be monotonicity-preserving for arbitrary , it takes many iterations to obtain a converged solution in both cases, and we generally find that the split algorithm is quite efficient by comparison. Our method might be able to serve as an effective preconditioner to further accelerate unsplit iterative methods, and works quite well as it is for our present purposes.

Thermal conduction is primarily along the magnetic field direction for hot plasmas. For astrophysical plasmas large temperature gradients exist, and it is important for the numerical scheme implementing anisotropic conduction to yield positive temperatures in regimes of interest. Another practical requirement is that the scheme be fast so that the conduction timestep is not much smaller than the MHD timestep. Here we present a simple, directionally-split method that gets close enough to monotonicity that the temperature remains positive in presence of relatively large temperature gradients, and results in a large speedup. Our method should find applications in modeling of hot astrophysical plasmas with large temperature gradients (e.g., the multiphase interstellar/intracluster medium, transition from chromosphere to corona in the Sun).

5 Acknowledgements

PS was supported by NASA through Chandra Postdoctoral Fellowship grant number PF8-90054 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060, and GWH was supported at the Princeton Plasma Physics Laboratory by DOE Contract No. DE-AC02-09CH11466. This research was supported in part by the National Science Foundation through TeraGrid resources provided by NCSA and Purdue University. Some of the runs were carried out on Henyey, the theoretical astrophysics computing cluster at the University of California, Berkeley. PS thanks Eliot Quataert for encouragement, and Jim Stone and Ben Chandran for useful discussions. We are grateful to Ian Parrish for discussions and for his comments on the paper. We thank the anonymous referees for very thorough referee reports that helped improve the quality of the paper substantially.


  • (1) V. Alexiades, G. Amiez and P. Gremaud, Super-time-stepping acceleration of explicit schemes for parabolic problems, Communications in Numerical Methods in Engineering 12 (1996), pp. 31-42.
  • (2) A. H. Baker, E. R. Jessup and T. Manteuffel, A technique for accelerating the convergence of restarted GMRES, SIAM Journal on Matrix Analysis and Applications 26 (2005), pp. 962-984.
  • (3) D. S. Balsara, D. A. Tilley and C. J. Howk, Simulating anisotropic thermal conduction in supernova remnants – I. Numerical methods, Monthly Notices of Royal Astronomical Society 386 (2008), pp. 627-641.
  • (4) P. J. Basser and D. K. Jones, Diffusion-tensor MRI: theory, experimental design and data analysis - a technical review, NMR in Biomedicine 15 (2002) pp. 456-467.
  • (5) Berkowitz, B., Characterizing flow and transport in fractured geological media: A review, Advances in Water Resources 25 (2002), pp. 861-884.
  • (6) L. S. Blackford, A. Cleary, J. Choi, E. D’Azevedo, J. Demmel, I. Dhillon, J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley, D. Walker and R. C. Whaley, ScaLAPACK Users’ Guide, SIAM, Philadelphia, PA (1997).
  • (7) S. I. Braginskii, Transport Processes in a Plasma, in: M. A. Leontovich (Ed.), Reviews of Plasma Physics, v. 1, Consultants Bureau, New York, (1965).
  • (8) Z. Dian-Lin, C. Shao-Chun, W. Yun-Ping, L. Li, W. Xue-Mei, X. L. Ma and K. H. Kuo, Anisotropic thermal conductivity of the 2D single quasicrystals: AlNiCo and AlSiCuCo, Physical Review Letters 66 (1991), pp. 2778-2781.
  • (9) R. Dong and J. M. Stone, Buoyant Bubbles in Intracluster Gas: Effects of Magnetic Fields and Anisotropic Viscosity, Astrophysical Journal 704 (2009), pp. 1309-1320.
  • (10) I. Gustafsson, A class of first order factorization methods, BIT 18 (1978), pp. 142-156.
  • (11) D. S. Kershaw, Differencing of the Diffusion Equation in Lagrangian Hydrodynamic Codes, Journal of Computational Physics, 39 (1981), pp. 375-395.
  • (12) R. J. Leveque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press (2002).
  • (13) K. Lipnikov, D. Svyatskiy and Y. Vassilevski, Interpolation-free monotone finite volume method for diffusion equations on polygonal meshes, J. Comp. Phys., 228 (2009), pp. 703.
  • (14) I. Nordbotten, I. Aavatsmark and G. T. Eigestad, Monotonicity of control volume methods, Numerische Mathematik, 106 (2007), pp. 255-288.
  • (15) I. J. Parrish and J. M. Stone, Nonlinear Evolution of the Magnetothermal Instability in Two Dimensions, Astrophysical Journal, 633 (2005), pp. 334-348.
  • (16) I. J. Parrish, E. Quataert and P. Sharma, Anisotropic Thermal Conduction and the Cooling Flow Problem in Galaxy Clusters, Astrophysical Journal 703 (2009), pp. 96-108.
  • (17) P. Perona and J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence 12 (1990), pp. 629-639.
  • (18) W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, Cambridge Univ. Press (1992).
  • (19) Y. Rasera and B. Chandran, Numerical Simulations of Buoyancy Instabilities in Galaxy Cluster Plasmas with Cosmic Rays and Anisotropic Thermal Conduction, Astrophysical Journal, 685 (2008), pp. 105-117.
  • (20) M. Ruszkowski and S. P. Oh, Shaken and stirred: conduction and turbulence in clusters of galaxies, Astrophysical Journal 713 (2010), pp. 1332-1342.
  • (21) Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM Journal on Scientific and Statistical Computing 7 (1986), pp. 856-869.
  • (22) P. Sharma and G. W. Hammett, Preserving monotonicity in anisotropic diffusion, Journal of Computational Physics 227 (2007), pp. 123-142.
  • (23) P. Sharma, P. Colella and D. F. Martin, Numerical Implementation of Streaming Down the Gradient: Application to Fluid Modeling of Cosmic Rays, submitted to SIAM Journal on Scientific Computing, arXiv0909.5426.
  • (24) P. Sharma, I. J. Parrish and E. Quataert, Thermal Instability with Anisotropic Thermal Conduction and Adiabatic Cosmic Rays: Implications for Cold Filaments in Galaxy Clusters, Astrophysical Journal, 720 (2010), pp. 652-665.
  • (25) Z. Sheng, J. Yue and G. Yuan, Monotone Finite Volume Schemes of Nonequilibrium Radiation Diffusion Equations on Distorted Meshes SIAM J. Sci. Comput., 31 (2009), pp. 2915.
  • (26) C. R. Sovinec et. al., Nonlinear magnetohydrodynamics simulation using high-order finite elements, Journal of Computational Physics 195 (2004) pp. 355-386.
  • (27) B. J. van Leer, Towards the Ultimate Conservative Difference Scheme. IV. A New Approach to Numerical Convection, Journal of Computational Physics, 23 (1977), pp. 276-299.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description