EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 615 (July 2018) A&A, 615 (2018) A67 Full HTML
Free Access
Issue
A&A
Volume 615, July 2018
Article Number A67
Number of page(s) 14
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201731916
Published online 16 July 2018

© ESO 2018

1 Introduction

Magnetic fields play an important role in the dynamics of many physical systems. The ideal magnetohydrodynamic (MHD) theory has been an excellent tool with which to understand the physical interaction between the magnetic field and the plasma in the solar atmosphere, with a long list of examples of successful modelling of a number of phenomena (e.g., Asplund et al. 2000; Schaffenberger et al. 2005; Vögler et al. 2005; Stein & Nordlund 2006; Khomenko & Collados 2006; Nordlund et al. 2009). Nevertheless, the assumptions underlying ideal MHD theory are not always valid. In the conditions of solar photosphere and low chromosphere, the presence of neutrals together with the weakening of collisional coupling due to the density fall off with height leads to the invalidation of the assumptions underpinning the MHD approach and produces a series of non-ideal effects. These non-ideal effects have been actively studied during recent years (e.g., Kumar & Roberts 2003; Khodachenko et al. 2004; Forteza et al. 2007; Pandey & Wardle 2008; Vranjes et al. 2008; Soler et al. 2009, 2010; Zaqarashvili et al. 2011; Khomenko & Collados 2012; Martínez-Sykora et al. 2012). The importance of the non-ideal effects can be evaluated by examining the ratio between the cyclotron frequency and the collisional frequencies. In the photosphere, the cyclotron frequency can exceed the collisional one in strongly magnetised areas such as flux tubes and sunspots, while in the chromosphere this happens over most of its volume even in regions with a weaker magnetic field (Khomenko et al. 2014). It can be expected that in such conditions the neutral and charged particles behave in a different way, with direct repercussions on the plasma dynamics and energy exchange. A recent review of the behaviour of partially ionised plasmas in different astrophysical conditions is provided by Ballester et al. (2017).

Mathematically, the treatment of partially ionised plasmas depends on the degree of collisional coupling. The most widely used formalisms are either a two-fluid or single-fluid approach, see for example the recent discussion in Khomenko et al. (2014) and a review by Ballester et al. (2017). The two-fluid formalism, including the treatment of ionisation-recombination terms was recently provided by Meier & Shumlak (2012). While the single-fluid formalism yields a less detailed description, it is numerically more efficient for strongly stratified solar plasma, and therefore it has been adopted in several numerical codes that aim at realistic simulations of the solar atmosphere such as MURaM (Vögler et al. 2005), Bifrost (Gudiksen et al. 2011), and MANCHA3D (Felipe et al. 2010), the latter being the one used in the current study. The single-fluid formalism has an advantage since numerically problematic collisional terms (both elastic collisions and ionisation-recombination) do not appear explicitly in the equations, greatly simplifying the modelling. For the range of collisional frequencies typical for the solar photosphere and chromosphere, these collisional terms require extremely short time steps, unsuitable for most frequently used explicit numerical codes. In the solar case, the single-fluid formalism can safely be used for the modelling of the dense and collisional photosphere and low chromosphere, for the typical time scales of interest (Zaqarashvili et al. 2011). For both the single and two-fluid formalisms the generalised Ohm’s law (Cowling 1957) has to be used to describe the behaviour of electric currents. The generalised Ohm’s law has a series of additional terms derived from the presence of neutrals, while other standard terms have slightly modified expressions (Khomenko et al. 2014). As a consequence, there appear additional terms in the induction and in the energy conservation equations. The leading non-ideal terms in the single-fluid formalism are: the Ohmic diffusion, the ambipolar diffusion, Hall term, and Biermann battery term.

In an astrophysical context, ambipolar diffusion is envisaged as a process of diffusion suffered by the magnetic field due to collisions between neutral and charged particles, the latter being frozen into the field. Ambipolar diffusion helps converting magnetic energy into thermal energy on scales that are typically much faster than classical Ohmic diffusion. For the Sun, it has been demonstrated that the ambipolar diffusion is orders of magnitude larger than Ohmic diffusion in the chromospheric layers, therefore greatly contributing to chromospheric heating via Joule dissipation (Khomenko & Collados Vera 2012; Priest 2014). In the two-fluid treatment, the energy equation for the charged fluid also has a Joule heating term, similar to the single-fluid case, and the expression for the electric field comes from the generalised Ohm’s law (see e.g., Khomenko et al. 2014; Ballester et al. 2017). It is interesting to note, however, that ambipolar diffusion term is not explicitly present in the two-fluid Ohm’s law, and in the energy conservation equation in the two-fluid treatment, and the corresponding effect is expressed via explicit ion-neutral collisional terms.

The Hall effect appears when ions decouple from the magnetic field while electrons remain attached to it. In a fully ionised plasma ions decouple from the magnetic field due to the difference between the inertia of ions and electrons, and this decoupling typically happens for processes at high frequencies. In weakly ionised plasmas, however, the ions can become detached from the magnetic field due to collisions with neutrals, leading to a similar effect. The spatial and temporal scales of the Hall effect in partially ionised plasmas are different from those in fully ionised plasmas, but the dynamics are similar in both cases. The Hall term is not dissipative, so it does not contribute to the heating of the environment, however it does redistribute magnetic energy in the system. Therefore, the Hall effect has an important role at many scales, from laboratory to astrophysical plasmas, while ambipolar diffusion is typically only important at low ionisation fractions (Parker 2007; Pandey & Wardle 2008; Priest 2014; Ballester et al. 2017).

The Biermann battery effect is somewhat different to the two non-ideal effects discussed above. Its magnitude in the solar case is rather small and it is usually ignored as it is not deemed important. But recently Khomenko et al. (2017) have shown that this term is able to seed the magnetic field for the local solar dynamo. In the present paper we will consider this term for completeness.

From a numerical point of view, the battery term only introduces a small source term into the hyperbolic system of equations of single-fluid MHD. So long as this source term is small in some sense it is not necessary for either accuracy or stability reasons to introducea time step restriction based on it. However, the ambipolar term is diffusive (i.e., parabolic), and the Hall term is dispersive. These two latter terms introduce important restrictions into the integration time step and into the stabilityof explicit numerical schemes, which are frequently used in hydrodynamic codes solving hyperbolic equations.

This paper presents a new numerical scheme aiming to overcome these restrictions. This scheme is implemented as separate modules in the existing single-fluid radiative MHD code MANCHA3D (Khomenko & Collados 2006, 2012; Felipe et al. 2010; Khomenko et al. 2017). The proposed scheme is relatively easy to add into an existing explicit numerical scheme, which allows one to keep all the benefits and flexibility of the explicit approach. We also provide a suite of numerical tests to verify the newly implemented modules and to check their order of accuracy.

This paper is organised as follows. In Sect. 2 we show the set of equations corresponding to a single fluid description with a generalised Ohm’s law. In Sect. 3 we present very briefly our single-fluid code MANCHA3D and the details of the new numerical scheme introduced here. In Sect. 4 we show tests demonstrating the feasibility and reliability of the new scheme introduced, while in Sect. 5 a brief summary is presented.

2 Equations

As discussed in the Introduction, the peculiarity of solar plasma in the photosphere and in the low chromosphere is that, while being only weakly ionised, it is still very dense and, therefore, collisions play an important role in coupling all plasma and neutral components. Since typically the scales of interest are significantly larger than collisional scales (see Khomenko et al. 2014), the most efficient approach for the mathematical treatment of such plasma is a single-fluid quasi-MHD formalism. The derivation of conservation equations and the generalised induction equation for such a case can be consulted elsewhere (see e.g., Ballester et al. 2017). Following the derivation of the equations from Bittencourt (2004), and using a specific definition for the macroscopic variables, see Appendix D of Khomenko et al. (2014), the non-ideal effects are encoded into several terms in the induction and total energy conservation equations. The set of equations is formed by the continuity equation (1)

the momentum conservation equation (2)

and the induction equation (3)

where we retained the convective, Ohmic, Hall, Ambipolar, and Biermann battery terms on the right hand side. The coefficients multiplying these terms are, respectively, the Ohmic (η), Hall (ηH), ambipolar (ηA), and battery (ηB) coefficients, and all have the same units of diffusivity (l2t−1, i.e. m2s−1 in SI). Those coefficients are defined later on, in Eqs. (9)–(11). The total energy conservation equation (4)

is written in terms of the total energy density per volume unit e, which is the sum of the kinetic, internal and magnetic energies, (5)

The electric currents are defined as (6) (7)

and Gauss’s law for magnetism is (8)

To close the system, the equation of state is used.

Equations (2) and (4) have two extra terms: an external time dependent force S(t) describing an arbitrary external perturbation, and a source term Qrad to take into account radiative energy exchange. These terms will be neglected in this work since the implementation of the numerical scheme introduced below is independent of the presence of these terms. The coefficients multiplying the non-ideal terms in the induction and energy equations are defined as (9) (10)

and (11)

The expressions for the collisional frequencies between ions and neutrals (νin) and electrons with neutrals (νen) are taken from Spitzer (1956). For the collisions between electrons and ions (νei) we use the expressions from Braginskii (1965) (12) (13) (14)

where min = mimn∕(mi + mn) and men = memn∕(me + mn) are the reduced masses. The cross sections for a weakly ionised plasma assuming elastic collisions between solid spheres are σin = 5 × 10−19 m2 and σen = 10−19 m2 (Vranjes & Krstic 2013; Huba 2016). The Coulomb logarithm, ln Λ, is defined as (15) (16)

The equations above assume charge neutrality (ne = ni) and negligible electron mass (me∕mi ≈ 0). A complete derivation of these equations and a discussion of the relative importance of the non-ideal terms along the solar atmosphere can be found in Khomenko et al. (2014).

3 Description and implementation of the new numerical scheme

3.1 MANCHA3D code

Before describing the MHDSTS1 scheme, we first briefly introduce our single-fluid MANCHA3D code. The current implementation of this numerical code is a parallel tridimensional extension of the code developed by Khomenko & Collados (2006) and solves the non-linear 3D MHD equations for the perturbations written in conservative form (Felipe et al. 2010).

To avoid high frequency numerical noise and improve the numerical stability some artificial diffusion terms are added instead of their physical counterparts following the work of Vögler et al. (2005), so each physical diffusivity has its own contribution formed by a constant contribution, a hyperdiffusivity part and a shock-resolving term. On the other hand, because in some simulations a high diffusivity can modify the perturbations’ amplitude (especially when solving waves) the code can performan additional filtering to dampen small wavelengths using the method described by Parchevsky & Kosovichev (2007). In this way, the artificial diffusivity remains low and the numerical noise is avoided. For spatial integration MANCHA3D uses a centred 6th order accurate derivative, following Nordlund & Galsgaard (1995). To advance in time it uses an explicit multi-step Runge-Kutta (RK) scheme, accurate to 4th order (Vögler et al. 2005; Khomenko & Collados 2006; Felipe et al. 2010). To prevent wave reflection at the boundaries, especially when working with complex magnetic field configurations, the code has implemented a perfectly matched layer (PML; Berenger 1994, 1996; Hu 1996, 2001; Hesthaven 1998; Parchevsky & Kosovichev 2009). It also has the possibility of using some adaptive mesh refinement (AMR) method (Berger & Colella 1989) to follow the evolution of fine structures.

More recently, a new radiative transfer module has been added to solve the radiative transfer equations using the short characteristic method. It is also possible to choose between an ideal equation of state or a realistic one in thermodynamical equilibrium based on precalculated tables, which include the ionisation equilibrium solution and the electron pressure needed to calculate the non-ideal terms, see Khomenko & Collados (2012) and Khomenko et al. (2017) for details.

The purely MHD version of the code has been verified against standard numerical tests by Felipe et al. (2010) and the code has been used so far for a wide range of numerical experiments (both ideal and realistic), such as modelling of wave propagation, instabilities, and magneto-convection (e.g., Felipe et al. 2010, 2013, 2016; Khomenko & Collados 2012; Luna et al. 2016; Santamaria et al. 2016; Khomenko et al. 2017). At this point, it is convenient to briefly explain the time integration method used by MANCHA3D as we shall refer to it in the following sections. In the single-fluid approximation, the system of equations written in conservative form can be schematically represented as (17)

where represents an operator formed by the sum of the spatial derivatives of fluxes F(u), the source terms, S(u), and the non-ideal terms (Ohmic, Ambipolar, Hall, and Battery), G(u). The vector u is the set of variables in the equations: [ρ, ρv, B, e](r, t).

The next time step is calculated using an explicit RK scheme that can be written in a compact form as

where u(n+1) corresponds to the solution at tn+1 = tn + Δt, with Δt given by the Courant, Friedrichs and Lewy (CFL) condition. The order of the RK scheme is given by the chosen m ≤ 4, and the coefficients αk can be found in van Leer et al. (1992, Table 6) or easily computed using the expression (19). (19)

MANCHA3D uses a fourth order RK scheme to solve the non-ideal single-fluid MHD equations, which using this notation corresponds to m = 4 in Eq. (18).

3.2 New MHDSTS scheme

As mentioned in Sect. 1, the non-ideal MHD terms are numerically problematic. In the solar atmosphere the value of the ambipolar diffusion coefficient (ηA) strongly increases with height and becomes larger than the classical Ohmic diffusion (η) by several orders of magnitude (Khomenko et al. 2014). When the ambipolar term dominates the equations, the system turns from hyperbolic to parabolic, and very small integration time steps are required in order to fulfil the CFL condition. To avoid such small time steps, the proposed MHDSTS scheme uses the technique known as super time-stepping (STS). This numerical technique was proposed by Gentzsch (1980) to speed up explicit time-stepping schemes for parabolic problems. Later on, Alexiades (1995) and Alexiades et al. (1996) presented a variant of this technique, which is used by O’Sullivan & Downes (2006); O’Sullivan & Downes (2007) to overcome the limitations caused by a large ambipolar term. As mentioned earlier, this numerical technique allows us to speed up, in a easy and efficient way, explicit time-stepping schemes for parabolic problems.

Besides, the Hall coefficient, ηH dominates over the other two coefficients (η and ηA) over most ofthe photosphere and then, the Hall effect can dominate the system (Khomenko et al. 2014). In this case, the system becomes skew-dominant, and the integration can not be done using the STS approach because it is unstable. In this situation, the system’s eigenvalues become complex and to fulfil the CFL condition we would require that our time-step gets close to zero to maintain stability (O’Sullivan & Downes 2006; O’Sullivan & Downes 2007; Gurski & O’Sullivan 2011). Thus, the proposed MHDSTS scheme uses yet another scheme, Hall diffusion scheme (HDS), which was proposed by O’Sullivan & Downes (2006) to overcome this problem.

The STS technique has been used in different fields successfully, for example: Choi et al. (2009) described the implementation of the ambipolar term into a multi-dimensional MHD code based in a total variation diminishing (TVD) scheme showing its use to follow complex MHD flows such as molecular cloud cores or protostellar discs; Grewenig et al. (2010) introduced STS into a new class of numerical schemes called fast explicit diffusion (FED) schemes to speed up the application of anisotropic diffusion filters for image enhancement or image compression tasks; Tsukamoto et al. (2013) presented an explicit scheme for Ohmic dissipation with smoothed particle magnetohydrodynamics (SPMHD) and used the STS technique to solve the Ohmic part of the induction equation; Tomida et al. (2013) also adopted the technique to solve the Ohmic part of the induction equation in problems of protostellar collapse; Wurster et al. (2014) provided an implementation of the STS technique into the PHANTOM-SPMHD code (Price & Federrath 2010a,b; Lodato & Price 2010; Price et al. 2017) for solving the ambipolar part of the induction equation. At this point we also mention that there is another family of STS scheme based on the Legendre polynomials and its properties which have been applied in problems involving thermal conduction in astrophysical plasmas as prominences or coronal rain successfully (Meyer et al. 2012, 2014; Xia & Keppens 2016; Xia et al. 2017).

Equation (17) can be split into different terms (ideal MHD, Ohmic diffusion, ambipolar diffusion, etc.), and these can be grouped together in a number of operators. There are a few options to group those terms into operators, but for our new MHDSTS scheme we have grouped them into three operators (MHD, STS, and HDS) as described in the following sections. The splitting technique we have chosen to solve the system is the Strang splitting. This splitting is written in a way that keeps asecond order accuracy along the simulation by permuting the order in which we apply the operators (Strang 1968; LeVeque 2002).

3.2.1 MHD operator

This operator is used to integrate the ideal MHD equations together with the battery term. We decided to include the non-ideal contribution of the Battery term into this solver because it does not limit the time step in solar-like problems (Khomenko et al. 2014). This operator can be defined as (20)

where the vector u is, as before, the set of conserved variables, LiMHD(u) = −∇F(u) +S(u) is (as per Eq. (17)) the standard operator for ideal MHD, and LBatt (u) is the operator for the battery term.

The temporal integration of this operator follows the same equations as Eq. (18), but in this case the operator being instead of and the time-step being ΔtMHD, calculated as (21)

where dtiMHD is obtained from the CFL condition applied to the operator LiMHD(u). The Batteryterm time step, dtBatt = min(dx2, dy2, dz2)∕ηB, does not play a roll into this operator due to its source term behaviour. It is important to point out that if there are no non-ideal terms this operator becomes a RK operator because there is no need to apply the Strang splitting. In this case, selecting m = 4, we recover the RK scheme. This standard part of the MANCHA3D code was throughly verified via either specific numerical tests, applied in 1D, 2D and 3D (Felipe et al. 2010), or through the scientific numerical simulations by the code published in the literature in the last ten years.

3.2.2 STS operator

This operator is used to treat only the contribution of the Ohmic and Ambipolar terms and it can be written as (22)

where the vector w corresponds to the conserved variables [B, e](r, t) or [B, eint](r, t) in the case where we prefer to use the internal energy equation, LOhmic(w) is the operator for Ohmic diffusion, and LAmbi(w) is the operator for ambipolar diffusion.

The main idea behind STS is to demand stability not at each time step, but at the end of a bigger step called super-step (Δ tSTS). The Δ tSTS is calculated as a sum of sub-steps τj. To maximise the size and stability of the super-step, the small sub-steps are obtained using the optimality properties of modified Chebyshev polynomials (Alexiades et al. 1996). (23)

in this expression, dtdif is the minimum time step given by the Ohmic term (dtOhmic = min(dx2, dy2, dz2)∕η) and by the ambipolar term (dtAmbi = min(dx2, dy2, dz2)∕ηA). NSTS is the number of sub-steps, and ν is a damping parameter, associated to the Chebyshev polynomials, to reduce higher frequencies and obtain a strong stability condition for the numerical method. This parameter is within the interval (0,λmin∕λmax], where λmin and λmax are the smallest and the bigger eigenvalues of the matrix respectively (Alexiades et al. 1996; Gurski & O’Sullivan 2011). When ν takes valuesclose to one, the system is highly damped, and the super-step becomes more accurate in the same way that any numerical explicit scheme has smaller errors resulting from small time steps. When ν takes values close to zero, the system becomes unstable.

Thus, the value of the parameters (NSTS, ν) control the stability, accuracy and speed of the method and depend on the spectral properties of the matrix . But, as pointed out by Alexiades et al. (1996), a precise knowledge of the eigenvalues is not needed to obtain a robust method, so the values can be arbitrarily chosen by the user considering certain limits (stability, speed, error, etc.). The time update by the STS operator can be written as (24)

where tn+1 = tn + ΔtSTS and the super-step ΔtSTS is given by (25)

with (26)

In the limit when ν tends to zero, the super-step formed by N sub-steps covers an N times bigger interval compared to N explicit time steps, namely: (27)

The STS scheme is a first order scheme and it belongs to the family of Runge-Kutta-Chebyshev (RKC) methods with N steps. Therefore, if used as described above, the overall precision of the MHDSTS scheme would be of first order, even if the precision of the other operators was bigger. In order to increase its order of accuracy one could perform a Richardson extrapolation (Richardson 1911; Richardson & Gaunt 1927) as suggested by O’Sullivan & Downes (2006). However, this implies many executions of the STS operator: three executions to reach second order, seven executions to reach third order and fifteen executions to reach fourth order.

For that reason, instead of performing a Richardson extrapolation, we consider the STS scheme as an ‘Eulerian’ first-order step of our multi-step RK scheme, so now we should reach fourth order just with four calls to the STS scheme. In other words, we use the STS scheme to obtain the solution at the stable points of the RK scheme and in this way, get a high order solution. We refer to this technique for increasing the order of a numerical scheme as RK-wrapper. Thus, the complete temporal scheme for the STS operator can be written similarly to Eq. (18) but, in this case, using the operator and its corresponding time step ΔtSTS:

Each step of this RK-wrapper has a size of αkΔtSTS, with αk given by Eq. (19). Each of these steps represents a whole STS super-step, where dtdif is scaled with αk and the corresponding τj values are calculated according to Eq. (23). With these, the conserved variables of the vector w(k−1) are calculated using the Expression (24).

3.2.3 HDS operator

The HDS operator only solves the Hall term in the induction equation. It is designed by O’Sullivan & Downes (2006) to overcome the problems originated by a skew-symmetric Hall term dominated system, and it can be written in conservative form as (29)

where the vector b is the conserved variable B(r, t), and LHall(b) is the operator for the Hall term. This operator works similarly to the MHD operator described above. However, unlike it, the update of the magnetic field components is done using all the information available at the moment. O’Sullivan & Downes (2006) applied the HDS scheme by using a first order RK scheme and then, increasing the accuracy order by the application of a Richardson extrapolation. Instead, we follow the same idea proposed in Sect. 3.2.2, treating the individual HDS operator as an “Eulerian” step to use the RK-wrapper to increase its accuracy order. By doing so, we keep the stability properties of this implicit-like scheme but reaching a higher accuracy order in time and with less computational effort. The proposed scheme would be:

In these equations dtHall is the time step imposed by the Hall term through the CFL condition given by O’Sullivan & Downes (2007), that is, . The HDS operator advances in time ΔtHDS = NHDSdtHall, by repeating the execution of the Eq. (30) NHDS times, which is the number of sub-cycles needed to cover one super-step Δ tSTS and/or one ΔtMHD; this is further elaborated in the following section.

If the order in which the Eqs. (30a)–(30c) are solved is maintained for successive time steps, an artificial handedness is introduced into the scheme. This could be avoid by reversing the order between time step or, as it was pointed by O’Sullivan & Downes (2007), performing a random permutation of the order in which the magnetic field components (30a)–(30c) are solved over successive time steps. This is especially important under certain circumstances, as for example if a strong directional bias is introduced in the initial state or in 3D numerical experiments. Such permutation can result in a small loss of numerical stability.

3.3 All together: timing and accuracy order

This new scheme is only worthwhile with respect to the RK when the Ambipolar or Hall terms introduce a heavy constraint in the time step as we commented in Sects. 3.2.2 and 3.2.3. Then, one should select the time steps for the operators according to the following rule: (31)

When the requirements of dtdif ≤ ΔtMHD is not fulfilled, the evolution of the system is dominated by hyperbolic ideal MHD terms, and there is no need to apply this new MHDSTS scheme, then we are out of the STS (and/or HDS) regime.

To choose the time step for each iteration, one first has to choose the pair (NSTS, ν) to keep the accuracy and stability during the simulation and also to try to maximise the size of the super-step along the simulation. One option for choosing those values can be obtained using Eq. (27) and forcing the STS operator to jump in time as far as the MHD operator allows it, so in the limit when ν → 0, we can write that (32)

getting in this way an optimum approximation for NSTS which should be very close to the number of sub-steps needed to cover the time step imposed by the ideal MHD CFL. Selecting the parameter ν can be a bit tricky, and for that reason it is convenient to define the STS efficiency as: (33)

This expression reaches an asymptotic value of 0.5 when NSTS →∞ for a given ν, see dashed line in Fig. 2. It can be also approximated as , see “+" markers in Fig. 2. Using this approximation for the efficiency we can obtain a value of ν to achieve a high efficiency and to maintain the stability and accuracy of the STS scheme. To do so, we can impose a certain efficiency and calculate, for a chosen NSTS what value of ν it corresponds to. After several tests we have concluded that a good compromise between speed and stability is reached fixing the STS efficiency around a 76% of its maximum. This is equivalent to equating tanh to one, in which case we can write . Finally, replacing the value of NSTS obtained from Eq. (32) we get an estimation for ν. The choice of number of HDS sub-cycles NHDS is done in execution time as .

The upper limit in time for the integration time step at any moment of simulation is given by the CFL condition applied to the ideal MHD operator ΔtMHD. The goal then is to adjust dtiMHD, dtdif and NHDS to fulfil the identity ΔtMHD = ΔtSTS = ΔtHDS (see Fig. 1 for an example).

Eventually, as we have mentioned, due to the Strang splitting, the overall order of the MHDSTS scheme is limited to second order whenever each operator reaches second order accuracy or higher using the RK-wrapper method. This accuracy is expected to decrease due to the different combinations of values we can choose for the pair (NSTS, ν) and the physical characteristics of each simulation.

thumbnail Fig. 1

Example of time step adjusting. Upper panel: situation before the adjust, all the super steps are different. Lower panel: the red square shows in each case the quantity that was changed. In this case, the code starts removing one sub-cycle of the HDS, then reduces d tdif to match the HDS super-step and finally reduces the MHD step to match the other two.
thumbnail Fig. 2

STS efficiency (ESTS) defined by Eq. (33) for different values of ν and itsapproximation (“+" markers). The dashed line indicates the asymptotic value of ESTS when NSTS →∞.

4 Numerical tests

As mentioned in Sect. 3.2.1, when all the non-ideal terms are switched off the MHDSTS scheme becomes a pure RK scheme because there is no need to apply the Strang splitting. This means that all the previous results obtained with MANCHA3D as well as the tests presented in Felipe et al. (2010) are also valid for this new version. Bearing this in mind, the sections below present specific numerical tests to only validate the newly introduced operators either individually, or working all together.

4.1 MHD operator test

For testing the MHD operator applied to the ideal part of the equations and the battery term, we performed the numerical experiment proposed by Tóth et al. (2012). In this test, a magnetic field is generated from scratch by means of the Biermann battery term in the induction equation. For this test, we have removed all the artificial diffusivities and considered only the contributionof the battery term.

The initial conditions for this test consist in a homogeneous isothermal background without magnetic field and no velocity fields (B =v = 0). To initiate the evolution, a small perturbations in electron density (Eq. (34)) and electron pressure (Eq. (35)) are applied.

The analytical solution (Eq. (36)) is obtained by introducing the initial condition into the MHD equations and integrating it one single time step. (36)

The numerical set-up consists in a 2.5D simulation box in the plane x − z. For this experiment we considered a double periodic domain in both x and z coordinates, covering a physical space of ± 10 meters, and using the constants, in S.I, n0 = 1∕e, p0 = 1, n1 = 0.1∕e, p1 = 0.1, kx = kz = 2π∕L, and L = 20.

The numerical solution is obtained using our MHD operator, Fig. 3, and it is compared with the analytical solution by calculating the absolute error between them, see Fig. 4. We only compare the numerical solution with the analytical one after one single time step because if the system keeps evolving the analytical solution (Eq. (36)) is not valid anymore, in other words, the convective term will contribute to the solution. It can be seen that most of the numerical domain does not exceed 10−9.3 G of absolute error and this is concentrated in the regions where the initial electron density and electron pressure perturbations are maximum.

thumbnail Fig. 3

Biermann battery test. The image shows the y component of the magnetic field vector at one time step (dt = 0.02 s) after the beginning of the simulation.
thumbnail Fig. 4

Logarithm of the absolute error in the By component of the magnetic field between the analytical solution and the numerical solution given by the MHD operator. The format of the figure is the same as Fig. 3. The error is minimum at locations with maximum electron density ne.

4.2 STS operator test

In order to test the implementation of the STS operator we use an experiment of Alfvén wave decay under the sole action of the ambipolar diffusion as it was proposed by Choi et al. (2009). This experiment is based on the analytical derivation of a dispersion relation characterising the propagation of an Alfvén wave in a homogeneous partially ionised plasma permeated by a horizontal and homogeneous magnetic field done by Balsara (1996) when a strong coupling approximation is taken into account. (37)

In this relation, is the Alfvén speed, k is a real wavenumberand ω = ωR + iωI is the complex angular frequency. It is known that the action of the ambipolar diffusion adds an imaginary term into the dispersion relation, which physically means that the Alfvén wave is unable to propagate when k ≥ 2cA∕ηA (Kulsrud & Pearce 1969). Under this strong coupling approximation and considering an ideal isothermal MHD system we can follow the evolution in time of a standing wave and compare the numerical results with the analytical solution of its first normal mode (Morse & Ingard 1968) (38)

where h0 is the initial wave amplitude.

For our experiments, we have considered an isothermal system in a 2.5D periodic domain in the plane x − z of side L with a constant ambipolar coefficient and a constant magnetic field B0 along the x coordinate, introducing a standing wave along the z coordinate with an initial velocity (39)

where vamp is a dimensionless initial peak amplitude. Taking k = kx in Eq. (37), it can be obtained that (40)

The analytical solution can be written as a superposition of two waves moving in opposite directions: (41)

with (42)

and tracking the evolution of the root-mean-square (RMS) magnetic field , the initial amplitude in Eq. (38) is .

We set the parameters vamp = 0.01, L = 1 m, T, p = 1 N/m, ρ = 1 kg/m3, and T = 1 K for the initial condition and background. We considered four cases with different constant ambipolar coefficients ηA = 0.01, 0.03, 0.1, and 0.3 m2 s−1.

First, we used NSTS = 1 and ν = 0 for the STS scheme to test how the error changes by increasing the order when using the operator splitting. In Fig. 5 we can see the results obtained in each case. In the upper panel of each subfigure we have the RMS evolution in time for the Bz component. In all cases the agreement between the analytical and numerical solutions is good, notice we only plot the numerical solution when the order is set to two. In the middle panel the relative error in time for the numerical solution obtained with the RK scheme and the STS scheme and three different orders is plotted. If we focus in the solution at order two, we can see how the shape and values of each schemes are very similar in time. The same can be seen for the other orders. Because the relative error goes to infinity when drops to zero a good tool to see the time evolution of the error can be the cumulative root-mean-square error (CRMSE) in time define as (43)

where n is the number of snapshots saved up to time tn, ti is the time of each output, and (Wurster et al. 2014). This quantity represents the sample standard deviation of the differences between the RMS value from our numerical experiment and the analytical solution and, as Hyndman & Koehler (2006) point, it is a good measure of the accuracy between models at the same scale. This function is plotted in the lower panel of the subfigures of the Fig. 5. When the rate at which the scheme adds error to the RMS is slower than the growth of the factor for increasing time, the CRMSE will decrease slowly with time to an asymptotic value. This is clearly seen when the wave is damped quickly due to the high ambipolar diffusion. Furthermore, it is clearly seen how both schemes have a similar behaviour for all the ambipolar diffusivities tested.

Now we study how the pair (NSTS, ν) affects to the error, computational time and accuracy order of the STS scheme in comparison with the RK. To do this, we made a couple of simulations using different values of them. We took three values of ν (0.1, 0.01 and 0.001) and for each of them, NSTS was taken between one and six.

The results are shown in Figs. 6 and 7 where in the upper panel of those figures we have plotted the CRMSE at a fixed time (t = 5τ, with τ = L∕cA) as function of (NSTS, ν). The time spent by each simulation to reach the same point in time (also at t = 5τ) is plotted at the middle panel. At the lower panel we show the accuracy order measured at t = 0.7τ. We chose this time to avoid values of the solutions close to zero (see vertical line marks in Fig. 5). All the tests where executed using a single node of the TeideHPC supercomputer (each TeideHPC node has 2 Intel Xeon E5-2670 processors, for a total of 16 cores per node).

In Fig. 6a, upper panel, we observe how the accumulated error for STS at order two gives results very similar to the RK one. However, for higher orders, the error is slightly smaller than the one corresponding to the RK schemes (horizontal dashed-lines). We see how increasing NSTS (i.e. increasing the number of times the STS operator will be called) introduces small extra errors. On the other hand we see how the CPU time increases linearly with NSTS due to these extra calls. Such behaviour is expected because our system is outside of the STS regime due to the low ambipolar diffusion and in this case, performing STS involves unnecessary computational effort. In Fig. 6b the system gets into the STS regime only for NSTS = 2, so for the rest of values we see the same behaviour as before. In the other cases shown in Fig. 7, ηA = 0.1 or 0.3, we have that the accumulated error increases over the error obtained for the RK1. However, by using the STS technique we obtain a speed-up of between two and eight times depending on the RK order chosen to compare. For NSTS > 4 the decreasing time tendency starts to invert, see Fig. 7c. This happens because in this particular test, the STS jump reaches the size of the MHD jump and the code starts to spend time with extra iterations of the STS scheme as in the previous cases.

Finally, following LeVeque (2002), to measure the accuracy order for conservation laws, we computed the L1 -norm against the analytical solution given by Eq. (41) as (44)

where the L1-norm is calculated at t0 = 0.7τ (see vertical dotted line in Fig. 5) and Np is the number of points into the numerical domain. The results are in general consistent with the selected order of the schemes and limitation imposed by the Strang splitting. We notice nevertheless, that for this experiment, the rounding errors are affecting strongly the determination of the order because we quickly reach a precise solution, even when working with low orders. Thus, an increase of the order does not improve the solution, and therefore the value of the error calculated according to Eq. (44) is meaningless. This problem of the error evaluation can be seen clearly for the RK schemes with high ambipolar diffusion.

In the upper panel of Fig. 7c we observe how the CRMSE gets values close to one in some cases. These cases correspond to the lower value of ν, and illustrate how the solution is entering into a unstable region, that is, the sub-steps τj are giving unstable values and the solution would not converge at the super-step. Checking the lower panel of the same figure, one can see how these points also correspond to the accuracy order values with no meaning. We conclude that one needs to be very careful in choosing the value of this damping parameter to avoid numerical instabilities. As we have seen, good values for this parameter can be obtained using Eq. (33). Another thing we notice is that for NSTS = 1 the behaviour of the CRMSE, CPU time or order matches the RK scheme very well. This is because when NSTS is one, the damping parameter ν is zero and then, the STS acts as a MHD operator, (i.e. a RK scheme), and then only the Strang splitting is affecting the results. In this case in particular, the results can be seen in Fig. 5 and verify that, despite the wrong values obtained for the accuracy order, see for example the lower panel of Fig. 7c (ηA = 0.1 m2 s−1), the numerical solution matches very well the analytic one.

For the STS we observe how by increasing NSTS or decreasing ν the CRMSE is increasing and the computation time and the accuracy order are decreasing. On the other hand, we see that by increasing the ambipolar diffusion coefficient the accuracy order tends to decrease slightly and both the CRMSE and the computational time tend to increase.

Despite the fact that STS is limited to second order due to the Strang splitting, we observe that setting the order of each operator individually at three or more, makes the whole system more stable. This way, we could capture better stiff problems and/or force the values of (NSTS, ν) to obtain a higher acceleration.

thumbnail Fig. 5

Upper panel: RMS value of the Bz componentof the magnetic field produced by the Alfvén wave, as a function of time. The solid black linecorresponds to the analytical solution and the crosses are the numerical solution obtained with the STS scheme at order two. It can be seen there is good agreement with the analytical solution for all the diffusivities tested. Middle panel: relative error in time for the RMS value of Bz component. The red, green and blue lines are for the numerical solution given by the RK scheme working in the second, third and fourth order, respectively. The dashed lines corresponds to the solution obtained with the STS operator working also in second, third and fourth order. Lower panel: CRMSE in time for both schemes defined by Eq. (43). As we expect, the behaviour of both schemes is similar as it is shown by the CRMSE. The time axis is in the units of the crossing time τ = L∕cA and the dotted vertical lines mark the time where the L1-norm is calculated, t0 = 0.7τ.
thumbnail Fig. 6

Upper panel: CRMSE values at t = 5τ s for different values of the pair (NSTS, ν). The horizontal dashed lines corresponds to the RK scheme values. The symbols correspond to the STS scheme, each of them indicate a different value of ν. The black, red, green and blue colours indicate first, second, third and fourth order accuracy respectively. The dotted lines linking the symbols are drawn just to make the plot more readable. Middle panel: CPU time spend by the schemes to reach the same point in the simulations. The colour and symbol code is the same as the upper panel. Lower panel: order computed using the L1 -norm at t0 = 0.7τ, see the dotted vertical line in Fig. 5; again, the colour and symbol code is the same as the upperpanel. In subfigure (a), we have a case with low ambipolar diffusion where the first order results unstable and the errors for the STS scheme are very close to the RK values. In this case, we have no acceleration because the system is out of the STS regime and the time increase linearly with NSTS as expected.
thumbnail Fig. 7

Same as Fig. 6, but showing a clear case inside the STS regime in subfigure (d). Here, the analytical solution is still well captured and the acceleration obtained is substantial. The errors, on the other hand, are limited, even thoughthey are higher than in the RK1 case.

4.3 HDS operator test

In order to test the HDS operator we have considered an experiment with a plane-polarised Alfvén wave propagating into a partially ionised plasma under the presence of a guiding magnetic field as it is proposed by Cheung & Cameron (2012). For this experiment we have neglected all the induction terms except the Hall term.

The numerical set up for the test consists in a 2.5D periodical box with a constant hydrodynamic background and a pure Alfvén wave propagating along the x-axis. The transverse components are considered small compare with the guide field, Bx.

Under these conditions, the Hall coefficient, ηH can be considered constant and the analytical solution given by Eqs. (45a)–(45c), can be obtained after linearised the single-fluid equations. This solution describes the precession of the polarisation plane of an Alfvén wave with a wavenumber k initially plane-polarised in the y-axis.

where ω is the angular frequency, B0 is a constant uniform magnetic field, b is a small perturbation and σ is the precession rate of the polarisation plane (46)

The relation between k and ω is obtained from the dispersion relation: (47)

For this experiment we select physical parameter to match those of the solar photosphere: ηH = 106 m2 s−1, B0 = 100 G, T = 6000 K, b = 0.1 G, pgas = 103 N m−2, ρ = 10−4 kg/m3, and L = 100 km so the domain covers one wavelength. For this set of parameters we have k = 2π × 10−5 m, σ =1.9 × 10−3 Hz, ω =5.6 × 10−2 Hz, and the rotational period τ = 2πσ−1 = 3183 s. For the initial conditions we set:

In the upper panel of Fig. (8) we can see the result of the simulations over the analytical solution given by Eq. (45). In the middle panel we can see how the relative error is higher near the nodes and, in the lower panel, how the CRMSEincrease in time as it is expected. We notice that simulations at first and second order are not plotted because they were not stable for the chosen CFL. However, at third and fourth order both schemes captures very well the analytical solution and they are very similar between them. Reducing the CFL slightly, the second order becomes stable.

thumbnail Fig. 8

Hall test, showing in the upper panel the time evolution for the numerical and analytical solutions of the Bz components of the magnetic field fluctuations due to the precession of the plane-polarised Alfvén wave in a fixed spacial point (x = 0). The solid line represents the analytical solution and the red crosses are the numerical solution. Middle panel: relative error of the numerical solution respect the analytical solution for both the RK and the HDS schemes, solid lines and dashed lines respectively. Lower panel: CRMSE in time for the Bz component. It can be seen how the two schemes have an almost identical behaviour in time and with the order. This it is expected because the HDS operator works as a RK operator. Also this shows us how well performs the Strang splitting.

4.4 MHDSTS scheme test

Finally, this section shows the results of testing simultaneously the MHD, HDS, and STS operators considering the effects of the Ohmic and ambipolar diffusions and the Hall term. For that, we used a shock-tube test with a magnetic precursor, or C-shocks (Draine 1980), under two different regimes depending of what term is dominating over the other non-ideal terms: a) Hall dominated and b) ambipolar dominated, see Falle (2003); O’Sullivan & Downes (2006); O’Sullivan & Downes (2007).

To obtain the steady-state solution we consider our system isothermal with all the three diffusivity coefficients (η, ηA, and ηH) constants and no battery term. We also set all the time derivatives to zero, this is equivalent to transform the equations toa frame of reference where the shock is steady. With this, our MHD equations can be written as

with the induction equation as

where C0, Cx, Cy, and Cz are constants in time. C1 and C2 are constants of integration obtained from applying the pre-shock boundary condition.

The initial states for the pre-shock and post-shock plasma were obtained from the jump conditions of a front shock propagating into a magnetised fluid. The parameter for this initial state are given in Table 1.

Given the density and velocity at the post-shock side we can determine the magnetic field from the induction equation at each spatial position and then use the other four algebraical equations to obtain the rest of variables. To solve this ordinary differential system of equations with non-constant coefficients, we use a fourth order RK scheme implemented into a different code.

With these parameters, and avoiding the use of the artificial diffusion terms, the calculation is unstable for the RK scheme at all orders due to the high frequency oscillation (whistler waves) at the shock front introduced by the Hall term. The only way of solving this problem with the RK scheme without artificial diffusion is by forcing a reduced CFL condition and increasing the filtering cadence. Then, as a side effect we will spend more CPU time and tend to over-smooth the solution. The whistler waves can be seen clearly at the upper panel of Fig. 9, where we can also see how the analytical solution is very well captured by the MHDSTS scheme. In this case, if we choose order two without artificial diffusion, the experiment results unstable but choosing order three the experiment is stabilised due to the slightly higher numerical diffusivity introduced by the extra RK-step to reach order three. In the lower panel of the same figure we can see the relative error between both solutions. For this experiment, as we have obtained a semi-analytical solution, it is easy to compute the accuracy order by changing the mesh resolution. In this case the obtained order is ~ 2.1. In Fig. 10, we show a second experiment with a higher ambipolar diffusion to activate properly the STS operator but keeping also a high Hall term contribution. We see again the good agreement between both solutions and also how by increasing the ambipolar diffusion enough, the selected global order of the code needed to get the solution is lower. In this case, the calculated accuracy order is ~ 1.6, slightly smaller than the case (a) due to the action of the STS scheme.

Table 1

Parameters for the MHDSTS test.

thumbnail Fig. 9

Hall dominated shock tube test. Upper panel: steady solution for the z-component of the velocity field. The left side corresponds to the post-shock side. Lower panel: relative error. The order obtained for this experiment is ~ 2.1.
thumbnail Fig. 10

Shock tube test with high ambipolar diffusion and comparable with the Hall term. Upper panel: steady solution for the z-component of the velocity field. The left side corresponds to the post-shock side. Lower panel: relative error. The order obtain from this test is ~ 1.6 a slightly lower than the experiment (a) due to the STS contribution as mentioned in Sect. 4.2.

5 Conclusions

In this paper we have seen that STS and HDS are numerical techniques easy to implement and also that they can work together using the Strang operator splitting formalism forming what we have called MHDSTS scheme. We also show how it is possible to increase the temporal order of accuracy using a RK formalism as a wrapper of each operator. However, MHDSTS is limited to second order by the Strang splitting, but setting the order of each operator individually at three or more makes the whole system more stable.

So far, all the tests performed have shown a good agreement with the analytical solution. This make us confident enough to use this newnumerical scheme with MANCHA3D code in a production mode to investigate the effects of those diffusion terms with stiff problems.

We also check that the STS technique can speed up problematic simulations by overcoming the CFL condition imposed by the parabolic term corresponding to the ambipolar term. But as a note of alert, we have to be aware that we must choose carefully the pair ν and NSTS because, as Alexiades et al. (1996) mentions, and we show with a numerical test, the smaller ν and/or bigger NSTS means bigger the acceleration but also bigger errors and the increase of numerical instabilities, as well as a slight decrease of the accuracy order.

On the other hand, we also check that the HDS technique solves the problem with complex eigenvalues introduced by Hall term being more stable than our standard Runge-Kutta scheme and avoiding that our dtHall becomes very small, specially when the Hall term is dominating.

Working with the MHDSTS scheme, we have that the HDS time step is improved when the ambipolar or Ohmic diffusion are considered.

Finally, due to the whole scheme being explicit, even when HDS behaves as an implicit-like scheme, it is straightforward to use the other already implemented capabilities of the MANCHA3D code such as parallel scalability, the PML boundary condition, the AMR capability or the radiative transfer module.

Acknowledgements

This work is supported by the Ministerio de Economía y Competitividad through project AYA2014-55078-P and the predoctoral training grants in Centres/Units of Excellence “Severo Ochoa” with reference SEV-2011-0187. It contributes to the deliverables identified in FP7 European Research Council grant agreement 277829, “Magnetic connectivity through the Solar Partially Ionised Atmosphere”. We also wish to acknowledge the contribution of Teide High-Performance Computing facilities to the results of this research. TeideHPC facilities are provided by the Instituto Tecnológico y de Energías Renovables (ITER, SA), http://teidehpc.iter.es.

References

  1. Alexiades, V. 1995, in 2nd International Conference on Dynamic Systems and Applications [Google Scholar]
  2. Alexiades, V., Amiez, G., & Gremaud, P. A. 1996, Com. Num. Meth. Eng., 12, 31 [Google Scholar]
  3. Asplund, M., Nordlund, Å., Trampedach, R., Allende Prieto, C., & Stein, R. F. 2000, A&A, 359, 729 [NASA ADS] [Google Scholar]
  4. Ballester, J. L., Alexeev, I., Collados, M., et al. 2018, Space Sci. Rev., 214, 58 [NASA ADS] [CrossRef] [Google Scholar]
  5. Balsara, D. S. 1996, ApJ, 465, 775 [NASA ADS] [CrossRef] [Google Scholar]
  6. Berenger, J.-P. 1994, J. Comput. Phys., 114, 185 [Google Scholar]
  7. Berenger, J.-P. 1996, J. Comput. Phys., 127, 363 [NASA ADS] [CrossRef] [Google Scholar]
  8. Berger, M. J.,& Colella, P. 1989, J. Comput. Phys., 82, 64 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bittencourt, J. A. 2004, Fundamentals of Plasma Physics, 3rd ed. (New York: Springer-Verlag) [CrossRef] [Google Scholar]
  10. Braginskii, S. I. 1965, Reviews of Plasma Physics, 1, 205 [Google Scholar]
  11. Cheung, M. C. M., & Cameron, R. H. 2012, ApJ, 750, 6 [NASA ADS] [CrossRef] [Google Scholar]
  12. Choi, E., Kim,J., & Wiita, P. J. 2009, ApJS, 181, 413 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cowling, T. G. 1957, Magnetohydrodynamics (New York: Interscience Publishers), 115 [Google Scholar]
  14. Draine, B. T. 1980, ApJ, 241, 1021 [NASA ADS] [CrossRef] [Google Scholar]
  15. Falle, S. A. E. G. 2003, MNRAS, 344, 1210 [NASA ADS] [CrossRef] [Google Scholar]
  16. Felipe, T., Khomenko, E., & Collados, M. 2010, ApJ, 719, 357 [Google Scholar]
  17. Felipe, T., Crouch, A., & Birch, A. 2013, ApJ, 775, 74 [NASA ADS] [CrossRef] [Google Scholar]
  18. Felipe, T., Braun, D. C., Crouch, A. D., & Birch, A. C. 2016, ApJ, 829, 67 [NASA ADS] [CrossRef] [Google Scholar]
  19. Forteza, P., Oliver, R., Ballester, J. L., & Khodachenko, M. L. 2007, A&A, 461, 731 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Gentzsch, W. 1980, in Proceedings of the Third GAMM–Conference on Numerical Methods in Fluid Mechanics (Wiesbaden: Vieweg+Teubner Verlag), 109 [CrossRef] [Google Scholar]
  21. Grewenig, S., Weickert, J., & Bruhn, A. 2010, in Pattern Recognition (Berlin, Heidelberg: Springer), 533 [Google Scholar]
  22. Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Gurski, K., & O’Sullivan, S. 2011, SIAM J. Numer. Anal., 49, 368 [CrossRef] [Google Scholar]
  24. Hesthaven, J. S. 1998, J. Comput. Phys., 142, 129 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  25. Hu, F. Q. 1996, J. Comput. Phys., 129, 201 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  26. Hu, F. Q. 2001, J. Comput. Phys., 173, 455 [NASA ADS] [CrossRef] [Google Scholar]
  27. Huba, J. L. 2016, Plasma Physics Division (Washington: Naval Research Laboratory) [Google Scholar]
  28. Hyndman, R. J., & Koehler, A. B. 2006, Int. J. Forecast., 22, 679 [CrossRef] [Google Scholar]
  29. Khodachenko, M. L., Arber, T. D., Rucker, H. O., & Hanslmeier, A. 2004, A&A, 422, 1073 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Khomenko, E., & Collados, M. 2006, ApJ, 653, 739 [Google Scholar]
  31. Khomenko, E., & Collados, M. 2012, ApJ, 747, 87 [NASA ADS] [CrossRef] [Google Scholar]
  32. Khomenko, E., & Collados Vera, M. 2012, ASP Conf. Ser., 463, 281 [NASA ADS] [Google Scholar]
  33. Khomenko, E., Collados, M., Díaz, A., & Vitas, N. 2014, Phys. Plasmas., 21, 092901 [NASA ADS] [CrossRef] [Google Scholar]
  34. Khomenko, E., Vitas, N., Collados, M., & de Vicente, A. 2017, A&A, 604, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kumar, N., & Roberts, B. 2003, Sol. Phys., 214, 241 [NASA ADS] [CrossRef] [Google Scholar]
  37. LeVeque, R. J. 2002, Finite Volume Methods for Hyperbolic Problems (Cambridge, UK: Cambridge University Press) [CrossRef] [Google Scholar]
  38. Lodato, G., & Price, D. 2010, MNRAS, 405, 1212 [NASA ADS] [Google Scholar]
  39. Luna, M., Terradas, J., Khomenko, E., Collados, M., & de Vicente, A. 2016, ApJ, 817, 157 [NASA ADS] [CrossRef] [Google Scholar]
  40. Martínez-Sykora, J., De Pontieu, B., & Hansteen, V. 2012, ApJ, 753, 161 [NASA ADS] [CrossRef] [Google Scholar]
  41. Meier, E. T., & Shumlak, U. 2012, Phys. Plasmas, 19, 072508 [NASA ADS] [CrossRef] [Google Scholar]
  42. Meyer, C. D., Balsara, D. S., & Aslam, T. D. 2012, MNRAS, 422, 2102 [NASA ADS] [CrossRef] [Google Scholar]
  43. Meyer, C. D., Balsara, D. S., & Aslam, T. D. 2014, J. Comput. Phys., 257, 594 [NASA ADS] [CrossRef] [Google Scholar]
  44. Morse, P., & Ingard, K. U. 1968, Theoretical Acoustics, 1st edn. (McGraw-Hill) [Google Scholar]
  45. Nordlund, Å., & Galsgaard, K. 1995, A 3D MHD code for Parallel Computers (Copenhagen: Astronomical Observatory) [Google Scholar]
  46. Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Sol. Phys., 6, 2 [Google Scholar]
  47. O’Sullivan, S.,& Downes, T. P. 2006, MNRAS, 366, 1329 [NASA ADS] [CrossRef] [Google Scholar]
  48. O’Sullivan, S.,& Downes, T. P. 2007, MNRAS, 376, 1648 [NASA ADS] [CrossRef] [Google Scholar]
  49. Pandey, B. P., & Wardle, M. 2008, MNRAS, 385, 2269 [NASA ADS] [CrossRef] [Google Scholar]
  50. Parchevsky, K. V., & Kosovichev, A. G. 2007, ApJ, 666, 547 [NASA ADS] [CrossRef] [Google Scholar]
  51. Parchevsky, K. V., & Kosovichev, A. G. 2009, ApJ, 694, 573 [NASA ADS] [CrossRef] [Google Scholar]
  52. Parker, E. N. 2007, Conversations on Electric and Magnetic Fields in the Cosmos (Princeton, USA: Princeton University Press) [Google Scholar]
  53. Price, D. J., & Federrath, C. 2010a, MNRAS, 406, 1659 [NASA ADS] [Google Scholar]
  54. Price, D. J., & Federrath, C. 2010b, Numerical Modeling of Space Plasma Flows, 429, 274 [NASA ADS] [Google Scholar]
  55. Price, D. J., Wurster, J., Nixon, C., et al. 2017 PASA, submitted, [arXiv:1702.03930] [Google Scholar]
  56. Priest, E. 2014, Magnetohydrodynamics of the Sun (Cambridge, UK: Cambridge University Press) [Google Scholar]
  57. Richardson, L. F. 1911, Phil. Trans. R. Soc. A, 210, 307 [NASA ADS] [CrossRef] [Google Scholar]
  58. Richardson, L. F., & Gaunt, J. A. 1927, Phil. Trans. R. Soc. A, 226, 299 [NASA ADS] [CrossRef] [Google Scholar]
  59. Santamaria, I. C., Khomenko, E., Collados, M., & de Vicente, A. 2016, A&A, 590, L3 [Google Scholar]
  60. Schaffenberger, W., Wedemeyer-Böhm, S., Steiner, O., & Freytag, B. 2005, ESA SP, 596, 65.1 [Google Scholar]
  61. Soler, R., Oliver, R., & Ballester, J. L. 2009, ApJ, 699, 1553 [NASA ADS] [CrossRef] [Google Scholar]
  62. Soler, R., Oliver, R., & Ballester, J. L. 2010, A&A, 512, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Spitzer, L. 1956, Physics of Fully Ionized Gases (New York: Interscience Publishers, Inc.) [Google Scholar]
  64. Stein, R. F., & Nordlund, Å. 2006, ApJ, 642, 1246 [NASA ADS] [CrossRef] [Google Scholar]
  65. Strang, G. 1968, SIAM J. Numer. Anal., 5, 506 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  66. Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2013, ApJ, 763, 6 [NASA ADS] [CrossRef] [Google Scholar]
  67. Tóth, G., van der Holst, B., Sokolov, I. V., et al. 2012, J. Comput. Phys., 231, 870 [Google Scholar]
  68. Tsukamoto, Y., Iwasaki, K., & Inutsuka, S.-i. 2013, MNRAS, 434, 2593 [NASA ADS] [CrossRef] [Google Scholar]
  69. van Leer, B., Lee, W. T., Roe, P. L., Powell, K. G., & Tai, C. H. 1992, Communications in Applied Numerical Methods, 8, 761 [NASA ADS] [CrossRef] [Google Scholar]
  70. Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Vranjes, J., & Krstic, P. S. 2013, A&A, 554, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Vranjes, J., Poedts, S., Pandey, B. P., & De Pontieu B. 2008, A&A, 478, 553 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Wurster, J., Price, D., & Ayliffe, B. 2014, MNRAS, 444, 1104 [NASA ADS] [CrossRef] [Google Scholar]
  74. Xia, C. & Keppens, R. 2016, ApJ, 823, 22 [Google Scholar]
  75. Xia, C., Keppens, R., & Fang, X. 2017, A&A, 603, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Zaqarashvili, T. V., Khodachenko, M. L., & Rucker, H. O. 2011, A&A, 529, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
1

MHDSTS is the abbreviation for MHD, HDS and STS.

All Tables

Table 1

Parameters for the MHDSTS test.

In the text

All Figures

thumbnail Fig. 1

Example of time step adjusting. Upper panel: situation before the adjust, all the super steps are different. Lower panel: the red square shows in each case the quantity that was changed. In this case, the code starts removing one sub-cycle of the HDS, then reduces d tdif to match the HDS super-step and finally reduces the MHD step to match the other two.
In the text
thumbnail Fig. 2

STS efficiency (ESTS) defined by Eq. (33) for different values of ν and itsapproximation (“+" markers). The dashed line indicates the asymptotic value of ESTS when NSTS →∞.
In the text
thumbnail Fig. 3

Biermann battery test. The image shows the y component of the magnetic field vector at one time step (dt = 0.02 s) after the beginning of the simulation.
In the text
thumbnail Fig. 4

Logarithm of the absolute error in the By component of the magnetic field between the analytical solution and the numerical solution given by the MHD operator. The format of the figure is the same as Fig. 3. The error is minimum at locations with maximum electron density ne.
In the text
thumbnail Fig. 5

Upper panel: RMS value of the Bz componentof the magnetic field produced by the Alfvén wave, as a function of time. The solid black linecorresponds to the analytical solution and the crosses are the numerical solution obtained with the STS scheme at order two. It can be seen there is good agreement with the analytical solution for all the diffusivities tested. Middle panel: relative error in time for the RMS value of Bz component. The red, green and blue lines are for the numerical solution given by the RK scheme working in the second, third and fourth order, respectively. The dashed lines corresponds to the solution obtained with the STS operator working also in second, third and fourth order. Lower panel: CRMSE in time for both schemes defined by Eq. (43). As we expect, the behaviour of both schemes is similar as it is shown by the CRMSE. The time axis is in the units of the crossing time τ = L∕cA and the dotted vertical lines mark the time where the L1-norm is calculated, t0 = 0.7τ.
In the text
thumbnail Fig. 6

Upper panel: CRMSE values at t = 5τ s for different values of the pair (NSTS, ν). The horizontal dashed lines corresponds to the RK scheme values. The symbols correspond to the STS scheme, each of them indicate a different value of ν. The black, red, green and blue colours indicate first, second, third and fourth order accuracy respectively. The dotted lines linking the symbols are drawn just to make the plot more readable. Middle panel: CPU time spend by the schemes to reach the same point in the simulations. The colour and symbol code is the same as the upper panel. Lower panel: order computed using the L1 -norm at t0 = 0.7τ, see the dotted vertical line in Fig. 5; again, the colour and symbol code is the same as the upperpanel. In subfigure (a), we have a case with low ambipolar diffusion where the first order results unstable and the errors for the STS scheme are very close to the RK values. In this case, we have no acceleration because the system is out of the STS regime and the time increase linearly with NSTS as expected.
In the text
thumbnail Fig. 7

Same as Fig. 6, but showing a clear case inside the STS regime in subfigure (d). Here, the analytical solution is still well captured and the acceleration obtained is substantial. The errors, on the other hand, are limited, even thoughthey are higher than in the RK1 case.
In the text
thumbnail Fig. 8

Hall test, showing in the upper panel the time evolution for the numerical and analytical solutions of the Bz components of the magnetic field fluctuations due to the precession of the plane-polarised Alfvén wave in a fixed spacial point (x = 0). The solid line represents the analytical solution and the red crosses are the numerical solution. Middle panel: relative error of the numerical solution respect the analytical solution for both the RK and the HDS schemes, solid lines and dashed lines respectively. Lower panel: CRMSE in time for the Bz component. It can be seen how the two schemes have an almost identical behaviour in time and with the order. This it is expected because the HDS operator works as a RK operator. Also this shows us how well performs the Strang splitting.
In the text
thumbnail Fig. 9

Hall dominated shock tube test. Upper panel: steady solution for the z-component of the velocity field. The left side corresponds to the post-shock side. Lower panel: relative error. The order obtained for this experiment is ~ 2.1.
In the text
thumbnail Fig. 10

Shock tube test with high ambipolar diffusion and comparable with the Hall term. Upper panel: steady solution for the z-component of the velocity field. The left side corresponds to the post-shock side. Lower panel: relative error. The order obtain from this test is ~ 1.6 a slightly lower than the experiment (a) due to the STS contribution as mentioned in Sect. 4.2.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.