MAGNUS: A new resistive MHD code with heat flow terms

TL;DR

This paper introduces MAGNUS, a new 3D resistive MHD simulation code incorporating heat transfer, validated through various tests, and capable of accurately modeling wave propagation and magnetic reconnection in the solar atmosphere.

Contribution

The paper presents MAGNUS, a novel resistive MHD code with heat transfer, utilizing advanced numerical schemes and divergence control, specifically designed for solar atmospheric simulations.

Findings

The code accurately preserves magnetic divergence within machine precision.

Validated with 1D and 2D tests for ideal MHD, resistivity, and thermal conduction.

Successfully simulated wave propagation and magnetic reconnection in solar conditions.

Abstract

We present a new magnetohydrodynamic (MHD) code for the simulation of wave propagation in the solar atmosphere, under the effects of electrical resistivity, but not dominant, and heat transference in a uniform 3D grid. The code is based on the finite volume method combined with the HLLE and HLLC approximate Riemann solvers, which use different slope limiters like MINMOD, MC, and WENO5. In order to control the growth of the divergence of the magnetic field, due to numerical errors, we apply the Flux Constrained Transport method, which is described in detail to understand how the resistive terms are included in the algorithm. In our results, it is verified that this method preserves the divergence of the magnetic fields within the machine round-off error. For the validation of the accuracy and efficiency of the schemes implemented in the code, we present some numerical tests in 1D and 2D for the ideal MHD. Later, we show one test for the resistivity in a magnetic reconnection process and one for the thermal conduction, where the temperature is advected by the magnetic field lines. Moreover, we display two numerical problems associated with the MHD wave propagation. The first one corresponds to a 3D evolution of a vertical velocity pulse at the photosphere-transition-corona region, while the second one consists in a 2D simulation of a transverse velocity pulse in a coronal loop.

Tables

Table 1. Table 1 Information of the numerical parameters used in the simulations

Test	Domain	Grid points	$C_{\text{cfl}}$	$\Gamma$	Slope limiter
Brio Wu	$[-0.5,0.5]\times [-0.5,0.5]\times [-0.5,0.5]$	$800\times 4\times 4$	0.25	5/3	All
Current Sheet	$[-0.5,0.5]\times [-0.5,0.5]\times [-0.5,0.5]$	$128\times 128\times 4$	0.25	5/3	WENO5
MHD rotor	$[0.0,1.0]\times [0.0,1.0]\times [0.0,1.0]$	$200\times 200\times 4$	0.1	1.4	MINMOD
Cloud Shock Interaction	$[0.0,1.0]\times [0.0,1.0]\times [0.0,1.0]$	$400\times 400\times 4$	0.05	5/3	WENO5
Magnetic Reconnection	$[-0.5,0.5]\times [-2.0,2.0]\times [0.0,1.0]$	$200\times 800\times 4$	0.05	5/3	WENO5
Thermal Conduction	$[-0.5,0.5]\times [-0.5,0.5]\times [0.0,1.0]$	$200\times 200\times 4$	0.05	5/3	MINMOD
Transversal Oscillations	$[-25.0,25.0]\times [0.0,1.0]\times [0.0,50.0]$	$400\times 4\times 400$	0.25	5/3	MINMOD
MHD Gravity Waves	$[-0.75,0.75]\times [-0.75,0.75]\times [-0.25,5.75]$	$105\times 105\times 210$	0.25	1.4	MINMOD

Table 2. Table 2 Error L 1 superscript 𝐿 1 L^{1} for the density in the Brio Wu test

Method	N	$L^1$	$L_{\text{order}}^1$
	50	2.46$\times {10}^{-2}$	-
	100	1.46$\times {10}^{-2}$	0.75
HLLC/	200	9.33$\times {10}^{-3}$	0.65
MM	400	5.31$\times {10}^{-3}$	0.81
	800	3.22$\times {10}^{-3}$	0.72
	1600	1.88$\times {10}^{-3}$	0.78
	50	2.12$\times {10}^{-2}$	-
	100	1.27$\times {10}^{-2}$	0.74
HLLC/	200	6.92$\times {10}^{-3}$	0.87
MC	400	3.59$\times {10}^{-3}$	0.95
	800	2.19$\times {10}^{-3}$	0.72
	1600	1.27$\times {10}^{-3}$	0.79
	50	2.10$\times {10}^{-2}$	-
	100	1.22$\times {10}^{-2}$	0.79
	200	7.64$\times {10}^{-3}$	0.67
HLLC/	400	3.95$\times {10}^{-3}$	0.95
WENO5	800	2.28$\times {10}^{-3}$	0.80
	1600	1.27$\times {10}^{-3}$	0.85
	50	2.56$\times {10}^{-2}$	-
	100	1.68$\times {10}^{-2}$	0.60
HLLE/	200	1.11$\times {10}^{-2}$	0.60
MM	400	6.54$\times {10}^{-3}$	0.77
	800	4.018$\times {10}^{-3}$	0.70
	1600	2.37$\times {10}^{-3}$	0.76
	50	2.50$\times {10}^{-2}$	-
	100	1.44$\times {10}^{-2}$	0.80
HLLE/	200	8.21$\times {10}^{-3}$	0.81
MC	400	4.07$\times {10}^{-3}$	1.01
	800	2.17$\times {10}^{-3}$	0.90
	1600	1.19$\times {10}^{-3}$	0.87
	50	2.04$\times {10}^{-2}$	-
	100	1.22$\times {10}^{-2}$	0.74
HLLE/	200	8.24$\times {10}^{-3}$	0.56
WENO5	400	4.43$\times {10}^{-3}$	0.90
	800	2.61$\times {10}^{-3}$	0.76
	1600	1.52$\times {10}^{-3}$	0.78

Table 3. Table 3 Errors of the Brio Wu test with Resistivity.

	$L^1$ error 1	$L^1$ error 2	$L^1$ order
$\eta =0.01$	1.41$\times {10}^{-3}$	9.70$\times {10}^{-4}$	0.54
$\eta =0.004$	1.45$\times {10}^{-3}$	8.09$\times {10}^{-4}$	0.84
$\eta =0.007$	1.36$\times {10}^{-3}$	8.27$\times {10}^{-4}$	0.71

Table 4. Table 4 L 1 superscript 𝐿 1 L^{1} self-convergence error of the Brio Wu test with thermal conductivity

	$L^1$ error 1	$L^1$ error 2	$L^1$ order
$\kappa =0.001$	1.68$\times {10}^{-3}$	8.63$\times {10}^{-4}$	0.96
$\kappa =0.0002$	2.07$\times {10}^{-3}$	1.06$\times {10}^{-3}$	0.96
$\kappa =0.0006$	1.86$\times {10}^{-3}$	8.58$\times {10}^{-4}$	1.08

Equations

$${\partial }_t\rho +{\partial }_i\left(\rho v^i\right)=0,$$

$${\partial }_t\left(\rho v^i\right)+{\partial }_j\left(\rho v^i{v}^j-\frac{B^i{B}^j}{{\mu }_0}+p_T{\delta }^{ij}\right)=-\rho g{\delta }^{iz},$$

$${\partial }_{t}E+{\partial }_j\left[\left(E+p_T\right)v^j-\frac{B^j(\vec{B}\cdot \vec{v})}{{\mu }_0}\right]$$

$$={\partial }_j\left\{q^j-{\left[\frac{\eta }{{\mu }_0}\vec{J}\times \vec{B}\right]}^j\right\}-\rho v^{z}g,$$

$${\partial }_t{B}^i+{\partial }_j\left(v^i{B}^j-B^i{v}^j\right)=-{\left[\nabla \times \eta \vec{J}\right]}^i,$$

$$\nabla \cdot \vec{B}=0,$$

$$E=\frac{\rho v^2}{2}+\rho e+\frac{B^2}{2{\mu }_0},$$

$$p=(\Gamma -1)\rho e,$$

$$\vec{q}=\kappa \nabla T,$$

$$\vec{q}=\kappa T^{5/2}\left(\vec{B}\cdot \nabla T\right)\vec{B}/B^2,$$

$${\partial }_t\vec{U}+{\partial }_i\overrightarrow{F^i}=\vec{S},$$

$$\vec{U}$$

$$\overrightarrow{F^i}$$

$$\displaystyle\vec{S}=\left[\begin{array}[]{c}0\\ -\rho\vec{g}\\ -\rho\vec{v}\cdot\vec{g}-\nabla\cdot\left(\vec{q}+\frac{\eta}{\mu_{0}}\vec{J}\times\vec{B}\right)\\ -\nabla\times\eta\vec{J}\end{array}\right].$$

$$\displaystyle\begin{array}[]{llll}\partial_{i}\rightarrow\frac{1}{l_{0}}\partial_{i}\,,&\partial_{t}\rightarrow\frac{1}{t_{0}}\partial_{t}\,,&l\rightarrow\frac{l}{l_{0}}\,,&t\rightarrow\frac{t}{t_{0}}\,,\\ \rho\rightarrow\frac{\rho}{\rho_{0}}\,,&\eta\rightarrow\frac{\eta}{l_{0}v_{0}\mu_{0}}\,,&\vec{q}\rightarrow\frac{\vec{q}}{\rho_{0}v_{0}^{3}}\,,&\vec{B}\rightarrow\frac{\vec{B}}{B_{0}}\,,\\ \vec{v}\rightarrow\frac{\vec{v}}{v_{0}}\,,&p\rightarrow\frac{p}{\rho v_{0}^{2}}\,,&p_{T}\rightarrow\frac{p_{T}}{\rho v_{0}^{2}}\,,&E\rightarrow\frac{E}{\rho v_{0}^{2}}\,,\end{array}$$

$$\frac{\mathrm{d}{\vec{U}}_{(i,j,k)}}{\mathrm{d}t}=-\frac{{\vec{F}}_{(i+1/2,j,k)}^x-{\vec{F}}_{(i-1/2,j,k)}^x}{\Delta x}$$

$$-\frac{{\vec{F}}_{(i,j+1/2,k)}^y-{\vec{F}}_{(i,j+1/2,k)}^y}{\Delta y}-\frac{{\vec{F}}_{(i,j,k+1/2)}^z-{\vec{F}}_{(i,j,k-1/2)}^z}{\Delta z}$$

$$+{\vec{S}}_{(i,j,k)},$$

$$\Delta t=C_{\text{cfl}}\times {\text{min}}_{ijk}\left(\frac{\mathrm{\Delta }x}{|{\lambda }_{ijk}^{n,x}|},\frac{\mathrm{\Delta }y}{|{\lambda }_{ijk}^{n,y}|},\frac{\mathrm{\Delta }z}{|{\lambda }_{ijk}^{n,z}|}\right),$$

$$\displaystyle\left.\begin{array}[]{c}\vec{U}_{(0,j,k)}=\vec{U}_{(1,j,k)}\\ \vec{U}_{(N_{x},j,k)}=\vec{U}_{(N_{x}-1,j,k)}\end{array}\right\}\text{for \ outflow \ BC},$$

$$\displaystyle\left.\begin{array}[]{c}\vec{U}_{(0,j,k)}=\vec{U}_{(N_{x}-1,j,k)}\\ \vec{U}_{(N_{x},j,k)}=\vec{U}_{(1,j,k)}\end{array}\right\}\text{for \ periodic \ BC},$$

$$\displaystyle\left.\begin{array}[]{c}v^{x}_{(0,j,k)}=-v^{x}_{(1,j,k)}\\ v^{x}_{(N_{x},j,k)}=-v^{x}_{(N_{x}-1,j,k)}\end{array}\right\}\text{for \ reflecting \ BC}.$$

$$J_x$$

$$J_y$$

$$J_z$$

$$\nabla \cdot {\left(\eta \vec{J}\times \vec{B}\right)}_{(i,j,k)}=$$

$$\frac{[\eta (J_z{B}_y-J_y{B}_z){]}_{(i+1,j,k)}-[\eta (J_z{B}_y-J_y{B}_z){]}_{(i-1,j,k)}}{2\Delta x}$$

$$+\frac{[\eta (J_x{B}_z-J_z{B}_x){]}_{(i,j+1,k)}-[\eta (J_x{B}_z-J_z{B}_x){]}_{(i,j-1,k)}}{2\Delta y}$$

$$+\frac{[\eta (J_y{B}_x-J_x{B}_y){]}_{(i,j,k+1)}-[\eta (J_y{B}_x-J_x{B}_y){]}_{(i,j,k-1)}}{2\Delta z},$$

$${\left(\nabla \cdot \vec{q}\right)}_{(i,j,k)}$$

$$c_{f,s}=\frac{1}{2\rho }\left(\Gamma p+\vec{B}\cdot \vec{B}\pm \sqrt{{\left(\Gamma p+\vec{B}\cdot \vec{B}\right)}^2-4\Gamma p{B}_x^2}\right).$$

$${\vec{U}}_L\text{if}{S}_L>0,$$

$${\vec{U}}^{\text{HLLE}}=\frac{S_R{\vec{U}}_R-S_L{\vec{U}}_L+{\vec{F}}_L-{\vec{F}}_R}{S_R-S_L}\text{if}{S}_L\le 0\le S_R,$$

$${\vec{U}}_R\text{if}{S}_R<0,$$

$${\overrightarrow{F^x}}^{\ast }=\frac{S_R{\vec{F}}_L^x-S_L{\vec{F}}_R^x-S_L{S}_R({\vec{U}}_R-{\vec{U}}_L)}{S_R-S_L}.$$

$${\vec{U}}_L\text{if}0<S_L,$$

$${\vec{U}}_L^{\ast }\text{if}{S}_L\le 0\le S^{\ast },$$

$${\vec{U}}_R^{\ast }\text{if}{S}^{\ast }\le 0\le S_R,$$

Full text

1. Introduction

The theory of magnetohydrodynamics -the study of interactions between magnetic fields and conductive fluids in low frequencies- is of great importance to understand the dynamics of the plasma in the solar atmosphere (Priest & Hood, 1991). Since the plasma in this region is highly influenced by strong magnetic fields, some of the oscillatory phenomena observed in its magnetic structures can be modeled as magnetohydrodynamic waves. For instance, some oscillations detected by TRACE in a coronal loop could be interpreted as a stationary kink mode (Nakariakov et al., 1999). Moreover, in Tomczyk et al. (2007) was proposed that Alfvén waves detected by CoMP (Coronal Multi-Channel Polarimeter) were a possible mechanism to heat the corona. However, Van Doorsselaere et al. (2008) discussed that these waves were better explained as kink magnetoacoustic modes.

On the other hand, according to some data obtained from the Atmospheric Imaging Assembly, on board the Solar Dynamics Observatory, Morton et al. (2012) observed a transversal kink mode propagating in a UV/EUV solar jet. In a similar way, De Pontieu et al. (2007) pointed out that some Alfvén waves detected by the Solar Optical Telescope, on board the Japanese Hinode satellite, could have the enough energy to accelerate solar winds. Moreover, using the data taken with this telescope, magnetohydrodynamic waves propagating along magnetic flux tubes in the solar photosphere were detected (Fujimura & Tsuneta, 2009). In addition, using the Swedish Solar Telescope (SST) and the Solar Optical Universal Polarimeter (SOUP), torsional Alfvén waves were detected in a bright-point group near the solar disk center (Jess et al., 2009). Also, from images of the Rapid Oscillations in the Solar Atmosphere instrument at the Dunn Solar Telescope, it was reported how oscillations in the pressure, in small-scale photospheric magnetic bright-points, are able to generate kink modes in the outer solar atmosphere (Jess et al., 2012). Complete reviews of observations of magnetohydrodynamic waves in solar regions like the corona, sunspots, prominences, coronal mass ejections, solar flares and solar winds can be found in Nakariakov & Verwichte (2005); Khomenko & Collados (2015); Okamoto et al. (2007); Vršnak et al. (2013); Shibata & Magara (2011); Ofman (2010), respectively.

Through numerical MHD simulations, it is possible to process and analyze the data taken by the solar missions. This is one of the reasons for which a big number of MHD codes have been developed. Many of them are very robust, with adaptive mesh refinement methods (AMR), several Riemann solvers, slope limiters and integrators. For instance: Athena (Stone et al., 2008), a grid-based code for astrophysical MHD, developed primarily for studies of the interstellar medium, star formation, and accretion flows; Flash (Fryxell et al., 2000), a high performance application code, which has evolved into a modular, extensible software system from a collection of unconnected legacy codes for diverse applications in hydrodynamics, MHD, radiation transfer, diffusion and conduction; Pluto (Mignone et al., 2007), a code for the solution of hypersonic astrophysical flows; Zeus (Norman, 2000), a numerical code for the simulation of fluid dynamical flows in astrophysics, including a self-consistent treatment of the effects of magnetic fields and radiation transfer; Enzo (Bryan et al., 2014), a code based on block-structured AMR for modeling astrophysical fluid flows; Nirvana (Ziegler, 2008), a parallel computational MHD code with AMR for astrophysical problems; MAP (Jiang et al., 2012), a code with AMR and parallelization for astrophysics. Moreover, some codes have different numerical schemes like the one developed by Arber et al. (2001), which is based on control volume averaging with a staggered grid, instead of using approximate Riemann solvers.

Furthermore, there are codes focused only on solar atmosphere phenomena, for instance: VAC (Shelyag et al., 2008), a code for the simulation of the interaction of an arbitrary perturbation with a magnetohydrostatic equilibrium background, designed to reproduce physical processes in stable gravitationally-stratified plasma; Surya (Fuchs et al., 2011), a high-resolution, well-balanced finite-volume-based massively parallel code, to study the propagation of waves in a stratified non-isothermal magnetic atmosphere; MURaM (Vögler et al., 2005), a 3D magnetohydrodynamics code, for applications in the solar convection zone and photosphere; LFM (Lyon et al., 2004), a global magnetospheric simulation code; Bifrost (Gudiksen et al., 2011), a massively parallel numerical code designed for the research of the stellar atmospheres from the convection zone to the corona; CAFE (González-Avilés et al., 2015; Lora-Clavijo et al., 2015), a 3D MHD code for the analysis of solar phenomena within the photosphere-corona region.

Multiple MHD scenarios of solar physics interest have been evolved with the codes described before, for instance, Murawski & Musielak (2010) have studied the small amplitude-Alfvén waves in the solar atmosphere, finding their spatial and temporal signatures in the linear regime. Konkol et al. (2010) evolved a 2D potential arcade in a gravitationally stratified corona, in order to explore the influence of thermal conduction on the attenuation of the fundamental stationary slow magnetoacoustic modes. Moreover, Murawski et al. (2011) simulated the formation of macrospicules created by a localized pulse below the transition region. On the other hand, the microflares in the corona and the chromosphere were modeled by Jiang et al. (2012), changing the bottom boundary conditions. Shelyag et al. (2011) tested the generation of small-scale vortex motions in the solar photosphere. Jess et al. (2012) simulated some oscillatory phenomena in the photosphere as periodic intensity fluctuations with magnetoconvection processes. In Vigeesh et al. (2012) a MHD wave propagation in the solar magnetic flux tubes was carried out by evolving a strong, axially symmetric magnetic flux tube in a stratified solar model atmosphere, whose lower boundary is located at photospheric levels. Later, Wedemeyer-Böhm et al. (2012) reported the imprints of convectively driven vortex flows in the photosphere and the corona; Shelyag et al. (2013) also analyzed plasma motions in photospheric magnetic vortices; Jelínek & Murawski (2013) reproduced magnetoacoustic-gravity waves, at the gravitationally stratified corona, by using a realistic temperature profile and curved magnetic field lines. Furthermore, Kelvin-Helmholtz instabilities at the coronal mass ejection boundaries in the solar corona were found by Möstl et al. (2013).

Many other numerical studies in this direction have been developed. In Baumann et al. (2013) some reconnection events in the corona have been reproduced in 3D. Impulsively generated fast magneto-acoustic perturbations were presented by Pascoe et al. (2013) using 2D evolutions. Moreover, in Murawski et al. (2013) 3D simulations of magnetohydrodynamic-gravity waves and the production of vortices in a magnetized solar atmosphere were tested. Later on, Wedemeyer et al. (2013) analyzed the propagating shock waves in 3D radiation simulations of M-type dwarf stars. Impulsively generated Alfvén waves have been also reproduced in an isolated solar arcade, which is gravitationally stratified and magnetically confined (Chmielewski et al., 2014). Recently, solar prominences embedded in sheared magnetic arcades were evolved by Terradas et al. (2015). Besides, the coronal heating problem was also discussed in more complex scenarios, including radiative MHD terms (Hansteen et al., 2015). In addition, transverse MHD waves in a prominence flux tube were taken into account by Antolin et al. (2015). Also, in Yang et al. (2015), the excitation of fast-propagating magnetosonic waves along coronal magnetic funnels was studied. In more recent works, 3D MHD models of solar flares were carried out by Janvier et al. (2015). Murawski et al. (2016) have also analyzed the response of a solar small-scale and weak magnetic flux tube with photospheric twisting motions. Finally, the footpoint-driven transverse waves propagating in a coronal plasma with a cylindrical density structure were considered by De Moortel et al. (2016).

Motivated by the works mentioned above, with the code described in this paper, our goal is to perform numerical simulations of the propagation of magnetohydrodynamic waves in the solar atmosphere, in scenarios with electrical resistivity and heat flow, although the amplitude of those terms will not be dominant. Specifically, this code solves the resistive MHD equations with heat flow terms and with a constant gravity in 3D. It uses the HLLE (Harten et al., 1983) and HLLC (Toro et al., 1994) Riemann solvers. The HLLC method is implemented following the adaptation proposed by Li (2005). The slope limiters are second-order MINMOD (Roe, 1986), MC (van Leer, 1977) and fifth-order WENO schemes (Titarev & Toro, 2004). In addition, in regard of preserving the magnetic field divergence-free, $\nabla\cdot\vec{B}=0$, we follow the Flux Constrained Transport method (Balsara, 2004; Evans & Hawley, 1988; Tóth, 2000), which was slightly modified in order to deal with the resistivity terms, and thus keep the divergence of the magnetic field to the order of the machine round-off error.

This paper is organized as follows: in section 2, we describe the equations to be solved numerically, that is, the magnetohydrodynamic equations with resistivity, heat flow, and constant gravity. Then, these equations are written in a conservative form, which is ideal for the numerical used here. In section 3, we present the numerical methods to solve these equations. Since the purpose of the code is to study wave phenomena in the solar atmosphere with resistive and heat flow terms, in section 4, we show the numerical tests to validate the accuracy and the efficiency of the algorithms. Finally, in section 5, we discuss the results and present our main conclusions.

2. Resistive Magnetohydrodynamic Equations with Heat Flux

The MHD equations for a resistive medium with heat flux and a constant gravitational field in the negative direction of the z-axis are the continuity of mass, the Newton’s second law for a fluid element, the conservation of energy and the Faraday’s law, given respectively by

[TABLE]

with the Gauss’s law for the magnetic field

[TABLE]

where $\rho$ is the mass density of the fluid, $p$ the pressure, $\vec{v}$ the velocity vector field, $\vec{B}$ the magnetic field, $\vec{J}$ the current density vector, $\eta$ the electric resistivity, $\vec{q}$ the heat flux vector, $g$ the acceleration due to gravity, $\mu_{0}$ the permeability of free space and $E$ the total energy density

[TABLE]

which is expressed as the sum of the kinetic, internal ($\rho e$) and magnetic energy densities. In the equation of energy (3), it has been assumed that the fluid obeys an equation of state for an ideal gas, that is,

[TABLE]

where $\Gamma$ is the adiabatic index.

Furthermore, the heat flow vector $\vec{q}$ for the thermal conduction is modeled in two ways, the first one corresponds to the isotropic dependence with the temperature, given by

[TABLE]

where $\kappa$ is the thermal conductivity, and $T$ is the temperature. The second model for the heat flux vector is the classical model for a magnetized plasma (Spitzer, 2006)

[TABLE]

where the thermal conduction is transferred only along the magnetic field lines. This approximation is suitable for solar applications, since the conduction in the perpendicular direction is $2\times 10^{-31}n^{2}T^{-3}B^{-2}$ (in SI units and $n$ the total number of particles per unit volume) times smaller (Priest, 2014) due to the strong magnetic fields. It is important to notice that with the inclusion of the heat flux vector, the equations have now a diffusive term, which makes more difficult the numerical calculations, so in regards to the stability of the solution, we choose a time step proportional to the square of the spatial resolutions.

The system of equations (1)-(4) can be written in the compact form

[TABLE]

where $\vec{U}$ is a state vector of conservative variables defined as

[TABLE]

$\vec{F^{i}}$ are the fluxes along each axis

[TABLE]

and $\vec{S}$ is a source vector

[TABLE]

With the purpose to achieve an optimal performance in the numerical processes of the code, the equations are adimensionalized with the following conversions

[TABLE]

where $v_{0}=l_{0}/t_{0}$ and $B_{0}=v_{0}^{2}\mu_{0}\rho_{0}$. With (25), the system of equations (10) can be re-written in the same form except for the term of the magnetic permeability $\mu_{0}$, which disappears.

3. Numerical methods

This code solves the time-dependent system of partial differential equations (10) in an uniform grid. The MHD equations are solved in time by using the method of lines, which is equipped with different time integrators. Specifically, it has implemented second and third order total variation diminishing Runge-Kutta, fourth order regular Runge-Kutta and iterative Cranck-Nicholson time integrators. These equations are then discretized following the finite volume approximation

[TABLE]

where $\vec{F^{x}}_{(i\pm 1/2,j,k)}$, $\vec{F^{y}}_{(i,j\pm 1/2,k)}$ and $\vec{F^{z}}_{(i,j,k\pm 1/2)}$ are the numerical fluxes at the interfaces of a cell and will be calculated using the High Resolution Shock Capturing methods (HRSC). Specifically they are HLLE, HLLC and ROE (only for hydrodynamics) schemes which are coupled with the slope limiters of second-order MINMOD (Roe, 1986), MC (van Leer, 1977) and the fifth-order weighted Essentially Non Oscillatory WENO5 (Titarev & Toro, 2004).

In order to obtain stable evolutions and to calculate the convergence, the time step is chosen following the Courant-Friedrichs-Levy condition (Titarev & Toro, 2005)

[TABLE]

where $C_{\text{cfl}}$ stands for the Courant number and $|\lambda_{ijk}^{n,d}|$ is the speed of the fastest wave present at time level $n$ traveling in the $d$ direction. It is worth mentioning that the Courant number is chosen depending on the space dimension, that is, $0<C_{\text{cfl}}\leq 1/2$ for 2D and $0<C_{\text{cfl}}\leq 1/3$ for 3D. Furthermore, in the cases where the resistivity and the heat flux are considered, the time step is chosen using the same formula but in this case the spatial resolutions are squared, since in those cases the equations have now a diffusive component.

On the other hand, the default four boundary conditions (BC) implemented in the code are outflow, periodic, reflecting and fixed in time. For instance, for the boundaries at the face $x=x_{\text{min}}$ and $x=x_{\text{max}}$, they are calculated as

[TABLE]

and for the reflecting we use outflow conditions in all the variables except for the velocity in the x-direction, which reflects as

[TABLE]

For the fixed in time conditions, we reset the values at the boundary with the initial values. Finally, we parallelized the code with OpenMP (Dagum & Menon, 1998) to improve the computational times of the calculations.

3.1. Finite differences for the current, resistive and heat flux terms

For the numerical calculations, it is necessary to compute the terms related to the resistivity and heat flux with finite differences. The components of the current density vector $\vec{J}$ are obtained from Ampére’s law

[TABLE]

and the partial derivatives are computed by using simple second-order finite differences. It is worth mentioning that we have dropped the term $\mu_{0}$ since the equations implemented are dimensionless. The fifth component of the source vector (21) and the heat flux vector $\vec{q}$ from equation (8) are calculated following the schemes presented in Jiang et al. (2012), which are given by

[TABLE]

[TABLE]

For the second model of the heat flux vector (9), simple second-order centered finite differences are used. Finally, the three last components of the source vector (21), $\nabla\times\eta\vec{J}$, are not written in a finite differences scheme, since the equations (4) will be solved through the Flux Constrained Transport technique to preserve the magnetic field divergence-free and will be discussed in section 3.3.

3.2. Approximate Riemann Solver

Since the purpose of the code is to study scenarios with low resistivity and low heat transport, that is, small amplitudes of the resistivity and thermal conductivity, we have assumed that the eigenvalue structure of the system of equations (1)-(4) is not altered by the inclusion of the resistivity and heat flux terms. This means the eigenvalues are the same that those obtained in Powell (1994) for the ideal MHD.

The HRSC methods compute the fluxes in equation (26) by solving a local Riemann problem on both sides of every inter-cell (with the slope limiters MINMOD, MC, etc…). Let us summarize both methods, the HLLE and HLLC for the 1D case (x-direction) since the 3D version can be easily extended by doing a cyclic permutation of the variables.

The HLLE scheme (Harten et al., 1983) is based on the eigenvalues $\lambda_{i}$ associated with the Jacobian matrix of the system of equations (10). The first one is an entropy wave traveling with speed $\lambda_{1}=v_{x}$, the second and third ones are two Alfvén waves traveling with speeds $\lambda_{2,3}=v_{x}\pm c_{a}$, and the other four are magneto-acoustic waves, two slow $\lambda_{4,5}=v_{x}\pm c_{s}$, and two fast $\lambda_{6,7}=v_{x}\pm c_{f}$, where $c_{a}=\frac{B_{x}}{\sqrt{\rho}}$ is the Alfvén velocity, $c_{f}$ and $c_{s}$ are the fast and slow magnetosonic speeds,

[TABLE]

The solution is approximated by three constant states separated by two waves, one at the left, one intermediate and one at the right side, given by

[TABLE]

where $\vec{F}_{L}=\vec{F}^{x}(\vec{U}_{L})$ and $\vec{F}_{R}=\vec{F}^{x}(\vec{U}_{R})$ are the fluxes evaluated in the left and right states respectively. $S_{L}$ and $S_{R}$ are the velocities of the fastest waves moving to the left and to the right, that is, $S_{L}=\text{min}(0,\lambda_{i}^{R},\lambda_{i}^{L})$ and $S_{R}=\text{max}(0,\lambda_{i}^{R},\lambda_{i}^{L})$. Finally, the fluxes for the intermediate state are given by

[TABLE]

An important advantage of the HLLE method is that it uses only the eigenvalue information of the system of equations, so it is easier to implement compared to others. However, this scheme dissipates more than others, since it ignores the contact discontinuity of the characteristic waves.

On the other hand, the HLLC scheme (Toro et al., 1994) for hydrodynamics, is a modification of HLLE method, where the contact discontinuity is solved. This method estimates an intermediate contact wave with velocity $S^{*}$, so the solution is divided into the following substates:

[TABLE]

Now, by applying the Rankine-Hugoniot conditions across $S_{L}$ and $S_{R}$, the fluxes can be written in terms of the state variables $\vec{U}_{L}^{*}$ and $\vec{U}_{R}^{*}$ as

[TABLE]

where it has been assumed that the pressure and velocities remain unchanged across the contact discontinuity wave

[TABLE]

The extension of this method to the MHD case is straight forward, but in the work of Li (2005) it was pointed out that the consistency condition (the integral form of the conservation laws Toro (2009))

[TABLE]

is not satisfied under the MHD extension. However, Li (2005) has shown that the preceding condition is satisfied with the following restriction for the magnetic fields

[TABLE]

in such a way that the intermediate state $\vec{U}_{L}^{*}$ is

[TABLE]

where

[TABLE]

and a similar solution for the right state vector $\vec{U}_{R}^{*}$ can be found.

3.3. Divergence of magnetic field: Flux Constrained Transport

When the MHD equations are solved numerically, the accumulation of errors makes the divergence of the magnetic field grow in time. Therefore, several techniques have been implemented to prevent this issue (Tóth, 2000). In our code, the Flux Constrained Transport method (Balsara, 2004; Evans & Hawley, 1988) is applied, which consists in a special discretization of Faraday’s equation. Now, from Maxwell equations, Faraday’s law can be re-written in terms of a vector $\vec{\Omega}$

[TABLE]

that according to the Ohm’s law, can be expressed as

[TABLE]

Applying finite central differences in (58), the discretized evolution equations for the area-average components of the magnetic field on each face of the numerical cell, centered in $(i,j,k)$, are

[TABLE]

[TABLE]

[TABLE]

where the values of $\Omega$ are defined at the cell vertex and computed, at the inter-cells, as simple averages of the numerical fluxes $F^{ij}=v^{i}B^{j}-v^{j}B^{i}$ ,

[TABLE]

[TABLE]

[TABLE]

with the fluxes $F^{ij}$ computed with the HLLE or HLLC schemes along each spatial direction. Furthermore, in the last equations, we have included the average of the second term in the right-hand-side of (59), $\eta\vec{J}$, since it is not presented in the definition of the fluxes. It is worth mentioning that we have chosen to express those terms in that way, so the method could be used with any approximate Riemann solver and preserve the $\nabla\cdot\vec{B}$ to the order of the machine round-off-error, as it will be shown in the next section for all the tests considered in this work.

The last step is to calculate the magnetic field in the center of the cell by averaging the values obtained in equations (3.3) to (3.3), that is

[TABLE]

Finally, we compute $\nabla\cdot\vec{B}$ in each step of the evolution using the cell corner centered definition, since the centered one is not conserved by the Constrainted Transport Method (Tóth, 2000)

[TABLE]

where the averages of the derivatives are given by

[TABLE]

4. Tests

With regard to analyzing the accuracy of the code, in this section, we present some tests to validate the implementation of the algorithms. We start with the 1D test Brio-Wu (Brio & Wu, 1988) in the regime of ideal MHD, which is compared with the exact solution to calculate its convergence. Later we repeat the same test with resistivity and with isotropic heat flux. Secondly, we depict three ideal MHD tests in 2D, the Current Sheet test (Hawley & Stone, 1995), the MHD rotor (Balsara & Spicer, 1999) and the Cloud Shock Interaction test (Dai & Woodward, 1998). Also, to test the resistivity terms and the second model for the heat flux (9) we run the 2D magnetic reconnection test and the thermal conduction test from Jiang et al. (2012). Finally, we implement two ideal MHD evolutions for wave propagation in the solar atmosphere, the first one is the 2D Transverse Oscillations in Solar Coronal Loops (Del Zanna et al., 2005), and the second one is a 3D simulation of gravity waves in the solar atmosphere using a realistic temperature profile (Murawski et al., 2013). All of those evolutions were done using a uniform grid in Cartesian coordinates. The results were obtained with the HLLC Riemann solver with different slope limiters, as indicated in each case in table 1, along with the information of the numerical domain, the number of points in the grid, the Courant number and the adiabatic index $\Gamma$. For the integration in time, we use the third order TVD Runge-Kutta. Finally, to check the accuracy and to verify that there are not non-physical phenomena, we present the maximum value of the divergence of the magnetic field in each of the 2D and 3D simulations.

4.1. Test 1: Brio-Wu

The Brio & Wu (1988) test is the extension of the Sod’s shock-tube problem with magnetic fields, which can verify the ability of the code to capture the characteristic waves of a MHD problem, that is, the shocks, rarefactions and contact discontinuity waves. The initial state, given by $p=1$, $v_{x}=v_{y}=v_{z}=B_{z}=0$, $B_{x}=0.75$ and

[TABLE]

consists of a left and right constant states. The parameters for the numerical simulation are listed in table 1. The results are presented in figure 1, where we display the density, pressure, velocity components $v_{x}$ and $v_{y}$ and the magnetic field component $B_{y}$ at the time $t=0.1$, using the MINMOD slope limiter. In this figure, we also plot the non-regular exact solution from Takahashi & Yamada (2013) in order to compare the results. Moreover, with the exact solution for the density, we calculate the $L^{1}$ norm of the error

[TABLE]

where $N$ is the number of grid points. To compute the order of the error, $L^{1}$ needs to be calculated with $2N$ points, then with them, the order is given by

[TABLE]

These two values are obtained for each one of the Riemann solvers and slope limiters implemented in the code. The results are displayed in table 2 for 50, 100, 200, 400, 800 and 1600 grid points.

4.2. Test 2: Brio-Wu with Resistivity

A version of the Brio-Wu Test with a uniform resistivity has been implemented. Using the same configuration than before, we run the test using different resistivities. In figure 2, we plot the density for each value of $\eta$, in order to visualize the effect of this term. It can be seen that for small values of the resistivity, the effect is a dissipation of the discontinuities, and for bigger values, the contact discontinuity is displaced to the left, meanwhile, the shock discontinuity moves to the right, as it is shown in the inset plots of figure 2.

Since the analytical solution, in this case, does not exist, we calculate a self-convergence of the $L^{1}$ norm of the error. To do this, we first calculate the error between the solutions with $N$ and $2N$ points, secondly, we estimate the error comparing the solutions with $2N$ and $4N$ points, and with both values, we calculate the order of the error with (80). This process is repeated for the resistivities, $0.01$, $0.004$ and $0.007$, which are depicted in table 3.

4.3. Test 3: Brio-Wu with Isotropic Thermal Conduction

To conclude with the Brio-Wu tests, we have implemented a new version of this problem by adding the effect of the thermal conduction model (8), where the heat flux is modeled for an isotropic medium, proportional to the gradient of the temperature . Repeating the simulation with the same values than before, we obtain the results of figure 3, where we plot the density for six different values of the thermal conductivity $\kappa$, $0.0$, $0.0002$, $0.0006$, $0.0008$, $0.001$ and $0.002$. We can note that the effect of the heat flux is a dissipation and attenuation of the discontinuities, which is more notorious at the contact wave, as it is shown in the inset plot; however no changes are presented in the global topology of the solution, as it was the case of the test with resistivity. Also, the self-convergence norms of the error are calculated for three different values of the thermal conductivity $\kappa=0.001,0.0002,0.0006$ and are depicted in table 4.

4.4. Test 4: Current Sheet

The first 2D test of ideal MHD is the Current Sheet problem (Hawley & Stone, 1995). This test is implemented in order to prove the robustness of the algorithms of integration in the code. The initial state consists in two current sheets with density $\rho=1$, pressure $p=0.3$, velocity components $v_{x}=0.1\sin{2\pi y}$, $v_{y}=v_{z}=0$ and the magnetic field components $B_{x}=B_{z}=0$,

[TABLE]

The parameters for the numerical simulation are listed in table 1. The boundary conditions are periodic in the planes $|x|=0.5$, $|y|=0.5$, and outflow in $|z|=0.5$. In figure 4, we plot a color map of the $B_{y}$ component and the magnetic field lines at times $t=0.5$, $t=2.0$, $t=3.0$ and $t=5.0$, where the magnetic reconnection can be tracked in time with the formation of magnetic islands that grow in time, causing a loss in the magnetic energy, which is converted into heat. From this results with low resolution, we can prove that the numerical methods implemented are good enough to describe standard non linear dynamics as the compressions and rarefactions raised in this test. The maximum divergence of the magnetic field in each step of this evolution is shown in figure 12(a).

4.5. Test 5: MHD Rotor 2D

The second 2D test is the ideal MHD Rotor, proposed by Balsara & Spicer (1999). It describes a rapidly spinning cylinder in a light fluid, with the following initial set up $p=1$, $B_{x}=5/\sqrt{4\pi}$, $B_{y}=B_{z}=v_{z}=0$,

[TABLE]

where $f(r)$ is a smooth function connecting the two regions, $r_{0}$ and $r_{1}$ are the inner and outer radius.

For the evolution, we have used the information given in table 1 and outflow boundary conditions in all directions. The results of the simulation are presented in figure 5, they correspond to the color map of the density with the magnetic field lines (top plot), and the color map of the pressure (bottom plot), both in the time $t=0.15$. Here we find a very similar behavior than the one obtained by other codes, see for instance Mignone et al. (2007), finding that the initial torsional Alfvén waves, sent by the rotor into the fluid around, reduce their angular momentum, followed by the compression of the magnetic field, which distorts its initial circular shape into an oval one. In addition in figure 12(b), we plot the maximum divergence of the magnetic field as a function of time.

4.6. Test 6: Cloud Shock Interaction

The third 2D test, is very interesting from the astrophysical point of view, since it describes the disruption of a cloud with high density with a strong shock wave. It was proposed by Dai & Woodward (1998) and has been widely used to check the algorithms behavior with superfast flows. Here we have implemented the version of Tóth (2000), given by the initial state $B_{x}=0=v_{y}=v_{z}=0$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the cloud is centered in $x=0.8$, $y=0.5$ and has a radius $0.15$. The discontinuity consists in a fast shock wave and a rotational discontinuity in the $B_{z}$ component. For the results presented here, we show the numerical information in table 1 and outflow boundary conditions in all the faces. In the figure 6, we plot the color map of the density with the magnetic field lines for the times $t=0.03$ and $t=0.06$. The maximum violation in the divergence of the magnetic field is plotted in figure 12(c).

4.7. Test 7: Magnetic Reconnection

In order to evaluate the implementation of the resistive terms in the algorithms, we reproduced a 2D test taken from Jiang et al. (2012). It consists in an application for the magnetic reconnection, given by an initial state at rest, $v_{x}=v_{y}=v_{z}=0$, with a density and pressure with constant values $\rho=1$ and $p=0.1$, and the magnetic field $B_{x}=0$,

[TABLE]

[TABLE]

where $L_{r}=0.05$ is a characteristic length of the region with resistivity, which is modeled as

[TABLE]

with an amplitude $\eta_{0}=0.05$ in the domain $(-L_{r},L_{r})\times(-4L_{r},4L_{r})\times(0.0,1.0)$. For the simulation, it was used the information listed in table 1 and outflow boundary conditions everywhere. In figure 7, we plot a density map of the thermal pressure with the magnetic field lines at times $t=0.2$, $t=0.5$ and $t=1.5$, where the phenomena of magnetic reconnection due to localized resistivity can be tracked in time. In early times the phenomena is hardly noticeable but it is intensified later. To evaluate the reconnection rate, we plot in figure 8 the ratio between the inflow velocity $v_{\text{in}}$ and the Alfvén speed, that is $v_{\text{in}}/c_{a}$, which is the reconnection rate Petschek (1964), calculated along the line at $y=2.5,|x|<0.25,z=0$. Here it was found that the rate tends to have a maximum value around $t=1$ as the resistivity is increased. On the other hand, once the system reaches a steady state, we can see that the bigger the resistivity, the bigger the reconnection rate. As it is expected, without resistivity there is no reconnection rate (solid black line). However, as the magnitude of the resistivity increases, the rate is fairly the same. Finally, to verify the violation of the Gauss’s Law for the magnetic field, we plot the maximum divergence of the magnetic field in figure 12(d).

4.8. Test 8: Thermal Conduction

For the evaluation of the heat flux model (9), we have implemented the thermal conduction evolution proposed by Jiang et al. (2012), where a small hot region in the bottom of the domain is advected by the magnetic field lines. The initial state has a constant pressure $p=0.1$, is at rest $v_{x}=v_{y}=v_{z}=0$, and with the magnetic field $B_{x}=B_{y}=1$, $B_{z}=0$, which is oriented forming a $45^{\circ}$ angle. The density is

[TABLE]

and describes a small hot region located in the bottom. The information for the simulation is given in table 1. The boundary conditions are periodic at $|x|=0.5$, inflow in the bottom surface $y=-0.5$ and outflow in the rest. The results are presented in figure 9, where the temperature is depicted in times $t=0.01$, $t=0.12$ and $t=0.4$. It can be seen that the heat flux advects the temperature gradient along the diagonal magnetic field lines. In addition, the maximum violation in the divergence of the magnetic field is plotted in figure 12(e).

4.9. Test 9: 2D - Transverse Oscillations in Solar Coronal Loops

Here we present the evolution of MHD wave propagation in a solar phenomenon, in order to evaluate if the code is able to reproduce this type of simulations performed by other authors. Following Del Zanna et al. (2005), we test the propagation of Alfvénic pulses in the solar coronal arcades with the initial magnetic field $B_{y}=0$,

[TABLE]

where $B_{0}=40$ G is the magnitude of the magnetic field at the photospheric level in the footpoints $x=\pm L/2$. The angular frequency is $k=\pi/L$ and $L$ is a characteristic length, $L=50$ Mm $=50\times 10^{3}$ km. On the other hand, the equation of hydrostatic equilibrium is given by

[TABLE]

and through the equation of state, the pressure can be written as $p=2k_{\beta}\rho T/m_{p}$, where $m_{p}$ is the mass of the proton and $k_{\beta}$ is the Boltzmann constant. The pressure can be solved from (127)

[TABLE]

where the temperature profile from the photosphere to the corona is modeled by the step function

[TABLE]

which locates the transition region at a height of $z_{t}=2$ Mm and with a width of $z_{w}=0.2$ Mm. The temperatures in the corona and photosphere are $T_{\text{cor}}=1.2\times 10^{6}$ K and $T_{\text{phot}}=6000$ K. The value of the density in the base of the corona is $\rho=10^{-15}$ g cm*-3* and the magnitude of the gravity acceleration is $g=274$ m/s. The fluid is at rest initially $v_{x}=v_{z}=0$, but its transverse velocity is perturbed by a pulse located in $x=0$ and $z=L/2$ as follows

[TABLE]

where $r_{0}=1$ Mm and the Alfvén speed at the corona is $v_{0}=1$ Mm/s.

For the simulation, we use the information listed in table 1 and outflow boundary conditions were implemented in all the faces. In figure 10, we present the transverse velocity $v_{y}$ (km/s) and the magnetic field lines at times $t=8$ s, $t=24$ s, $t=64$ s and $t=80$ s. In this plots, we can follow the propagation of the twin Alfvénic pulses along the magnetic field lines and the fast and slow magneto-acoustic waves in the outer and inner arcs. Furthermore, we observe the reflection with the bottom surface as was described in Del Zanna et al. (2005). The maximum violation of $\nabla\cdot\vec{B}$ (in Tesla/km) is plotted in figure 12(f), where we can see how its initial value is not absolute zero, since it is calculated with the approximation of finite differences for the initial values of the magnetic fields (125) and (126).

4.10. Test 10: 3D - Gravity Waves

Until this point we have only shown evolutions in 1D and 2D, therefore we present here the 3D simulation of gravity waves in the solar atmosphere, which was performed using the Flash code (Fryxell et al., 2000) by Murawski et al. (2013). The initial state is given by a constant magnetic field in the $z-$ direction

[TABLE]

with a magnitude of $B_{0}=23$ G. The equation of equilibrium is given by (127) and through the equation of state, for this particular configuration, the pressure can be written as $p=k_{\beta}\rho T/(1.24m_{p})$. In a similar way as the Transverse Oscillations test, the pressure can be solved integrating the temperature

[TABLE]

where $p_{0}$ is the pressure at the level of reference in $z_{0}=10$ Mm at the solar corona. Assuming a temperature profile similar to VAL-IIC (Vernazza et al., 1981), the pressure and density can be obtained for the equilibrium state. The perturbation of the velocity is a vertical Gaussian pulse

[TABLE]

with amplitude $A_{v}$ = 3 km/s, width $\omega=100$ km and is located in $x=0$ and $z_{0}=500$ km.

For the simulation, we use the information given by table 1 with boundary conditions fixed in time at $z=-0.25$ Mm, $z=5.75$ Mm, and outflow in the other faces. In figure 11, we plot the vertical velocity $v_{z}$ at times $t=35.71$ s, $t=71.43$ s, $t=107.14$ s and $t=157.14$ s in the plane $y=0$. In this plots, we can track the propagation of the longitudinal magnetoacoustic-gravity waves caused by the initial perturbation in the vertical velocity $v_{z}$. The fast and slow modes are coupled and expanded in a quasi-isotropic way across the z-direction. Finally, the maximum violation in the divergence of the magnetic field (in Tesla/km) is plotted in figure 12(g).

5. Discussion and conclusions

We have presented in this work a resistive MHD code with heat flux terms. Here we have summarized how the numerical methods were implemented, emphasizing in the HLLC Riemann solver and the slight modification of the Flux Constrained Transport method, insights of the inclusion of the resistivity. We have shown that this code is able to reproduce the standard MHD tests in 1D and 2D since we have obtained the expected behavior of the solutions in all the cases. The resistivity and heat conduction were verified in an independent way and proved to be capable of displaying the expected phenomena, for instance, the magnetic reconnection and the dissipation of the energy. Furthermore, evolutions of MHD waves in solar physics problems were carried out to evaluate the efficiency of the numerical schemes in their description.

We conclude that one of the bests qualities of our results is the successful implementation of the Flux Constrained Transport method, since we managed to preserve $\nabla\cdot\vec{B}$ in the order of the machine round-off error, $1\times 10^{-12}$, as was displayed for the 2D and 3D tests. This is especially important in the cases where resistivity and heat flux were considered, since it shows that the numerical scheme was well adapted for them. We would like to point out that the numerical techniques need to be carefully adjusted to describe the physics laws, and this constraint should not be sacrificed in order to program a faster code. This is the reason why we preferred to build a unigrid code, where no unphysical behaviors could appear when AMR technics are implemented. Finally, it is worth mentioning that the code has been programmed in such a way that new approximate Riemann solvers, reconstructors, divergence-free methods and time integrators, can be easily added.

Acknowledgments

We want to thank professor K. TAKAHASHI for providing us with the exact solution of the Brio-Wu problem. A. N wants to thanks the financial support from COLCIENCIAS, Colombia, under the program “Becas Doctorados Nacionales 647” and Universidad Industrial de Santander. F.D.L-C gratefully acknowledges the financial support from VIE-UIS, grant number 1822. G.A.G. was supported in part by VIE-UIS, under Grants No. 1347 and No. 1838, and by COLCIENCIAS, Colombia, under Grant No. 8840.