CAFE-AMR: A computational MHD Solar Physics simulation tool that uses AMR

The set of MHD equations we are using to describe the dynamics of a fluid element are the following:

We assume the thermal pressure $p$ of the plasma is given by the equation of state of an ideal gas so that the total energy density is

In order to maintain the solenoidal constraint, that is $\nabla\cdot\mathbf{B}=0$, we implement the divergence cleaning scheme proposed in (Dedner et al., 2002), where the divergence of $\mathbf{B}$ is related to the scalar function $\psi$ that satisfies equation (5). With this correction, numerical errors that produce a non zero divergence of the magnetic field are transported and diffused through the numerical domain, and at the same time the addition of the equation for $\psi$ maintains the hyperbolic-parabolic structure of the equations. The coefficients $C_{h}$ and $C_{p}$ in Eq. (5) determine the damping of the divergence of the magnetic field.

The system of equations (1), (2), (3), (4), and (5) is solved as an Initial Value Problem (IVP) defined on the domain described in Cartesian coordinates $D\times[0,t_{f}]:=[x_{min},x_{max}]\times[y_{min},y_{max}]\times[z_{min},z_{max}]\times[0,t_{f}]$, provided initial and boundary conditions are given. We define a base cell centered discretization of $D$ as follows:

For the system of equations to be solved we define the vector of conservative variables ${\bf U}=(\rho,\rho{\bf v},\varepsilon,{\bf B},\psi)^{T}$. The finite volume implementation uses a Godunov type of shock capturing flux computation scheme, which prescribes the evolution formula for ${\bf U}$ from time $t^{n}$ to $t^{n+1}$ as:

Time evolution is implemented using the method of lines; we have two second order implementations to do so. A TVD second order Runge-Kutta method consistent with cell centering in the space-time volume with an appropriate reconstructor for the inter-cell boundary values; among the reconstructors available in CAFE-AMR, in the tests below we use the linear piecewise second-order Minmod, methods. The other time evolution scheme we have worked out leans on the corner transport Upwind (Colella, 1990) that implements the GLM-MHD equations, using the schemes proposed in (Mignone & Tzeferacos, 2010).

Following the guidelines in (Dedner et al., 2002), for the hyperbolic cleaning, the divergence present in the source term is implemented directly into the equations with the addition of the exact solution for the normal component of the magnetic field of the intercell face and $\psi$. We use the operator splitting on the parabolic part of the source term, which has an exact solution.

The selection of the sectors within the numerical domain to be refined need a refinement criterion. There is a certain liberty to define this criterion conveniently, according to the physical scenario that will be simulated. Typically, an arbitrary number of positive finite scalar functions $\chi(\mathbf{U}^{\ell})$ are defined in every cell of the coarse grid, then every cell $i,j,k$ is flagged to be refined if this scalar function has a value greater than a certain threshold $\chi_{r}$:

Among others, in this code we implement two principal refinement criteria:

• •

  Choosing as the evaluating functions the same numerical approximation of the gradients of every conservative variable as the one used to linearly interpolate their values in a new refined grid. That is:

  $$\chi_{m}=\overline{\nabla}(U_{m}),$$

  (11)

  where $U_{m}$ is the $m-th$ conservative variable.

• •

  The refinement criterion used in (Mignone et al., 2011), based on the second derivative error norm defined in (Löhner, 1987), where

  $$\chi(\mathbf{U})=\sqrt{\frac{\sum\limits_{d\in\{x,y,z\}}|\Delta_{d,+1/2}\sigma-\Delta_{d,-1/2}\sigma|^{2}}{\sum\limits_{d\in\{x,y,z\}}(|\Delta_{d,+1/2}\sigma|+|\Delta_{d,-1/2}\sigma|+\epsilon\sigma_{d,ref})^{2}}},$$

  (12)

  and $\sigma\equiv\sigma(\mathbf{U})$ is a function of the conservative variables, $\Delta_{d,\pm 1/2}\sigma$ are the undivided forward and backward finite differences along the direction $d$, for example

  $$\Delta_{x,\pm 1/2}\sigma=\pm(\sigma_{i\pm 1}-\sigma_{i}).$$

  The last term that appears in the denominator is given by:

  $$\sigma_{x,ref}=|\sigma_{i+1}|+2|\sigma_{i}|+|\sigma_{i-1}|,$$

  $$\sigma_{y,ref}=|\sigma_{j+1}|+2|\sigma_{j}|+|\sigma_{j-1}|,$$

  $$\sigma_{z,ref}=|\sigma_{k+1}|+2|\sigma_{k}|+|\sigma_{k-1}|,$$

  adding the $\sigma_{d,ref}$ term multiplied by the factor $\epsilon<1$ prevents the refinement of regions where small ripples are found.

Additionally, when a simulation uses more refinement levels, in order to properly nest them after refinement level $\ell$, cells are flagged according to the refinement criterion, cells in the level $\ell+1$ that contain cells of level $\ell+2$ are also flagged, and finally, level $\ell$ cells that have flagged neighbors of the same level are flagged as well. This procedure is made recursively from the most refined grid to the coarsest one.

A mentioned before, the refinement factor we use this work is $r=2$, so that a coarse cell will be covered with $2^{d}$ cells of the refined domain, where $d$ is the dimension of the domain. In Figure 1 it is shown how the cell volume of refinement level $\ell+1$ is $1/2^{d}$ of a cell in level $\ell$.

The refined cells will be centered in positions such that the region covered by them has the same boundaries as the coarse cell. For example, in two dimensions, four refined cells that cover a coarse one will have their centers at positions:

In what follows we use $D^{\ell}$ to denote the domain covered with the grid of refinement level $\ell$.

Likewise in (Porth et al., 2014) and (Matsumoto, 2007), data reconstruction onto the refined grid is done by a linear interpolation, which lowers the numerical diffusion of shock discontinuities, as well as prevents the generation of spurious shocks on cells that abut on coarser ones. Explicitly, the conservative variables $\mathbf{U}$ in the domain $D^{\ell+1}$ are:

The approximation of the gradient operator, $\overline{\nabla}$, acting on the arbitrary scalar function $\varphi$, is the one presented in (Matsumoto, 2007), given by:

The AMR strategy we implement generates multiple self similar grids, each one of them evolving independently in time by the scheme (7) and only sharing information through their boundaries at the beginning of the evolution from time $t^{n}$ to $t^{n+1}$.

The Godunov type of finite volume time evolution used for every grid has the Courant-Friedrichs-Lewy (CFL) stability condition that relates the time step $\Delta t$ with the grid resolution by

Following the indications in (Berger & Oliger, 1984) and (Berger & Colella, 1989), time step size $\Delta t^{\ell+1}$ of refinement level $\ell+1$ is related with $\Delta t^{\ell}$ of refinement level $\ell$ by

A complete evolution cycle of the system from $t^{n}$ to $t^{n+1}$ is considered as done once all of the grids synchronize at the time $t^{n+1}$. This requires that the domain $D^{\ell}$ must be evolved $2^{\ell}$ times to catch up with the coarsest domain. Figure 2 shows how grids of multiple refinement level advance several times to do so.

In (Keppens et al., 2003), it was found that choosing a non-homogeneous refinement ratio independently of the overall scheme can optimize certain aspects of the simulation, leading to faster computing times. However, our observations indicate that the gain from this approach, especially with a refinement ratio greater than 2, is outweighed by the benefits of using the recursive code structure of the fixed ratio. The recursive approach facilitated our control of proper refinement nesting and hierarchical time stepping.

Various AMR codes prove that in many physical simulations, homogeneous time step is preferred, as seen in the Nirvana code (Ziegler, 2008), which employs a single time step, not only for the better handling of simulations with non zero diffusion source terms, but also to realize faster parallelization schemes. In our tests, hierarchical time stepping yielded good results.

Except for the boundaries of the physical domain, the evolution of refined domains $D^{\ell}$ starts with the implementation of boundary conditions. This is done by filling ghost cells that extend the domain $D^{\ell}$, filled with appropriate flux values, see e.g. (LeVeque, 2002). Filling of ghost cells depends on the kind of boundary they are extending. As shown in Figure 3, we identify three types of boundaries of domain $D^{\ell}$:

• •

  Shared boundaries with another domain of the same refinement level $\ell$.

• •

  Shared boundary with domains of lower refinement level, which may be only of level $\ell-1$ because of how new sub grids are generated.

• •

  Shared boundary with the physical boundary.

According to these types of boundaries ghost cells are filled in the following way:

• •

  If ghost cells have the same position as cells contained on another same refinement level grid, simply copy data values contained on those cells.

• •

  If ghost cells of a domain $D^{\ell}$ are contained on a coarse cell from a domain $D^{\ell-1}$ linearly interpolate the data values.

• •

  If ghost cells correspond to a physical boundary, implement the appropriate boundary conditions, such as outflow or rigid, to name a few.

In the case of boundaries between consecutive refinement levels, the values of the variables at ghost cells of $D^{\ell}$ are constructed using linear interpolation from the values at cells in $D^{\ell-1}$ in the same way as equation (14), that is

where it is necessary to define the coarse grid values $\mathbf{U}^{*\ell}$, as

Once the variables in the coarse grid and its refinement are synchronized, the values of the conservative variables in the refined grid must be injected onto the coarse one. This is done by overwriting the data values contained in the coarse cell, with the average of the data values in the refined cells that cover it.

In three dimensions this is done using the formula

Finite volume methods are based on the conservation of quantities at a given cell. For example, the conservation of mass equation (1) on its finite volume form, equation (7), at cell $(i,j,k)$ can be written as:

Equation (26) establishes that each of the numerical fluxes introduces an amount of mass gain or loss across the faces of a cell. We focus our attention on the face boundary at $x=x_{i+1/2}$ where the mass gain is

By expressing (28) in terms of (27), an equivalence of the numerical fluxes of the cells in refinement level $\ell-1$ and the cells in refinement level $\ell$ across the shared boundary is exposed:

In order to enforce (29) we define a flux correction given by

This code works with a self similar block structure, as seen in Figure 5, where every block has the following data:

• •

  A grid discretized with a regular spacing using $N_{x}\times N_{y}\times N_{z}$ cells plus ghost cells used at the boundaries.

• •

  The physical location of the center of the block.

• •

  Allocatable arrays. These may contain the values of the conservative variables, or the flux corrections given in (31) to the cells abutting on the boundaries of the block.

• •

  A refinement level $\ell$ flag that indicates the physical volume of the cells on its grid as well as the time step $\Delta t^{\ell}$.

• •

  Pointers to other blocks that share physical boundaries which have the same refinement levels. This blocks are referred to as neighbors.

• •

  Pointers to other blocks referred to as children. These specific blocks cover a quadrant, in two dimensions, or an octant, in three dimensions, contained inside the original block.

• •

  Pointer to another block which has a portion of its domain covered by this block. This block is referred to as father.

• •

  Several logical flags. They have various uses, such as defining boundary conditions and deciding whether memory will be used, to name a few.

All the blocks are ordered hierarchically in a QuadTree, in two dimensional problems, and an OcTree in three dimensional problems. In Figure 6 a QuadTree is represented. In this data hierarchy, the refinement of the domain is generated in the following way. First the blocks are analyzed by quadrants:

1. 1.

   Every time a cell of a given quadrant satisfies the refinement criteria (10), the whole quadrant is flagged to be refined.

2. 2.

   If a Block’s quadrant associated child has children of its own, the quadrant is flagged to be refined.

3. 3.

   If a quadrant is flagged to be refined, all of its neighbors are also refined, even if they correspond to a different block.

This whole procedure is done from the highest refinement level to the base level so that the blocks are properly nested, according the scheme in subsection 3.1. Once the refinement process has ended and all data have been interpolated, the whole neighbor pointer structure is updated. Octrees are refined in the same way.

Finally, if a data block is flagged to abut on a coarser one, a flux correction array is allocated. This array adds on the corrections in order to calculate the terms of Eq. (31). Flags that indicate whether coarse data will be interpolated to fill ghost cells in children blocks are also activated at the end.

As the data structure suggests, time evolution is done recursively beginning on all the blocks of a given base level to an user specified last level of refinement. The time evolution control flow that holds synchronization of all blocks is the following:

1. 1.

   Apply physical boundary conditions at refinement level $\ell$.

2. 2.

   Copy data values from same level neighbors that share a boundary into ghost cells.

3. 3.

   Fill ghost cells values to level $\ell+1$ blocks abutting level $\ell$ blocks by interpolation in a synchronous way following equation (23).

4. 4.

   Advance level $\ell$ blocks a time-step $\Delta t^{\ell}$.

5. 5.

   Advance recursively $\ell+1$ blocks a time-step $\Delta t^{\ell}/2$.

6. 6.

   Fill ghost cell values to level $\ell+1$ blocks abutting level $\ell$ blocks by interpolation in an asynchronous way following equation (23).

7. 7.

   Advance recursively $\ell+1$ blocks a time-step $\Delta t^{\ell}/2$.

8. 8.

   Inject values of children blocks into the corresponding region they cover on level $\ell$ blocks using equation (25).

9. 9.

   Apply numerical flux corrections on level $\ell$ blocks following (31).

Flux corrections are implemented by first collecting in allocatable arrays, both, the cells abutting on $\ell-1$ grids and the corresponding cells on coarse blocks, afterwards these data are used to construct the terms $\delta F^{\ell-1}$ of equation (30).

Kelvin-Helmholtz instabilities are standard tests for hydrodynamical codes, and in the context of Solar Physics these can occur at the boundaries between different regions of the solar wind with varying velocities or densities as well as in small scale solar jets (Skirvin et al., 2023). As shown in (Mishin & Tomozov, 2016). These instabilities can lead to the formation of turbulent structures and the mixing of plasma, which can affect the propagation and properties of the solar wind. This is a MHD test, in which the discontinuities in counter-flow regions develop spiral-like shapes on, that are a challenge for an AMR method to follow.

The initial conditions for this test are set to:

This test is defined in the periodic domain $[-0.5,0.5]\times[-0.5,0.5]$, with a base resolution of $160\times 160$ cells, time step resolution given by a CFL factor $0.75$. For this test we use the Roe flux formula and the CTU time stepping scheme. We present three simulations with different refinement criteria as well as different number of refinement levels.

In a first simulation we use criterion (12) with a threshold value of $\chi_{r}=0.1$ using the density $\rho$ as the evaluating function in (10). We use five refinement levels, which gives an equivalent resolution of $5120\times 5120$ cells. In Figure 7 we show a snapshot along with a zoomed region at the interface region, the addition of resolution across the interface manages to develop spiral like shapes in small spatial scales. In Figure 8 we illustrate how the AMR method follows the discontinuity even with a slim refined domain while also maintaining the overall structure of the less dynamical regions.

In Figure 9 we show the results of a second simulation where we use three refinement levels, giving an equivalent resolution of $1024\times 1024$ cells. In this case, we compare the use of refinement criterion (11) with $\chi_{r}=10$ and criterion (12) with $\chi_{r}=0.05$, acting on the density $\rho$. In the refinement structure, it is shown that, since the criterion (11) uses all of the conservative variables as evaluating functions to set the refinement regions, the refined grids are set to cover all of the mixing regions whilst in criterion (12) the refinements are located in smaller regions. This results in a difference on small scale structure of the simulation. It is important to keep in mind this in order to define appropriate refinement strategies.

In presence of acceleration, when a fluid in low and high density phases, with the low density supporting the weight of the high density phase, initially in hydrostatic equilibrium, small perturbations give rise to the hydrodynamic Rayleigh-Taylor instability.

The analysis of this type of instability is used to learn how an initial linear perturbation develops a complex non-linear behavior. Rayleigh-Taylor instabilities in solar physics serve to gain insights into the physical processes responsible for solar eruptions, energy release, and the dynamics of the Sun’s atmosphere, (see e.g. Jenkins & Keppens (2022)). At first, in the linear regime, the growth rate of these instabilities depends on the wave number and frequency of the perturbation, for example the amplitude grows faster for shorter wavelengths. After the perturbation evolves into the nonlinear regime, bubbles from the low density region rise and drops from the high density shell fall, crossing the initial interface surface and the dynamics leads to the mixture of the two initial phases. Its complexity is a reason why this is a challenging test of fluid dynamics codes.

In the case of the MHD, the fluid is considered a magnetized ideal gas that responds to a magnetic field, that changes the evolution of the instability (see e.g. (Stone & Gardiner, 2007)). In this test the magnetic field is set parallel to the interface between the two fluid states and inhibits the growth rate of the instability. In order to show that CAFE-AMR can handle this problem, we present the two dimensional Rayleigh-Taylor problem.

We solve the equations on the domain $x\in[-0.5,0.5]$, $y\in[-1.5,0.5]$, and define the interface as the line $y=0$. Initial conditions assume the density in lower and higher parts of the domain are $\rho_{L}=1$ and $\rho_{H}=4\rho_{L}$ respectively, the gravity constant is set to $g=-1.0$ and the system is in hydrostatic equilibrium, consistently with the pressure

We use periodic boundary conditions along the $x-$direction and fixed conditions at the top and bottom of the domain. The base resolution uses $64\times 128$ cells and we use five refinement levels, which results in an equivalent resolution of $2048\times 4096$ cells with the finest resolution. The refinement criterion used is that in Eq. (12) using hydrodynamic pressure, density and the $\hat{y}$ component of the velocity as evaluating functions together with a threshold parameter $\chi_{r}=0.1$. Numerical fluxes are constructed using the Roe flux formula and the minmod limiter. Time resolution is given by a CFL factor of $0.75$.

The perturbation that triggers the instability is added to the velocity field at the interface line, assuming the profile

Figure 10 presents the evolution of the density and refinement structure profiles of the instability at $t=1.5,2.0,3.0$ for the unmagnetized case. As the denser shell begins to fall along the vertical direction, and due to the difference in velocities of the low and high density regions, the flow develops Kelvin-Helmholtz instabilities which in turn develop further into a highly structured interface at latter times. The refinement strategy selected on this case adequately tracks the overall structure of the interface.

Figure 11 presents the density and refinement structure profiles in the slightly magnetized case $\alpha=0.06$, the Kelvin-Helmholtz instabilities are suppressed because the topology of the field lines interrupts the torsion of the fluid interface, the magnetic field works against gravity and helps the denser shell to remain in the upper part which also stops the generation of KH-instabilities. In Figure 12 we show that the magnetic pressure is acting overall only at the interface surface. In Figure 13 we show that because of the augmented precision, we can capture at the lowest region of the interface there is a small scale variation in both density and magnetic pressure.

Solar flares are modelled as events generated magnetic reconnection (see e.g. (Priest, 2014)). In (Yokoyama & Shibata, 2001) two dimensional MHD simulations were employed to show the effect that localized magnetic resistivity has at triggering this process. The AMR strategy used to approach such simulations is advantageous as locality of phenomena is important. With this in mind, we present the localized resistivity Current Sheet (LRCS) test. This test models the generation of flares in a stratified density profile were a current sheet magnetic field is developed with the presence of a localized resistivity with Gaussian profile. The simulation is carried out on a numerical domain of $(x,y)\in[-10,10][0,20]$; fixed boundary conditions at the bottom boundary and outflow boundary conditions on the other sides of the domain; we use an initial resolution of $80\times 80$ cells with four refinement levels giving an equivalent resolution of $1280\times 1280$; the refinement criterion used is (11) with a threshold of $\chi_{r}=1$; numerical fluxes are constructed using the HLLE formulas with the second order CTU time stepping. A $CFL=0.01$ factor was chosen so that time steps capture difussion time scales adequately.

Initial data is defined as follows. Density has the stratified profile

In Figure 14 it is shown that near the maximum of the resistivity there is an increase in magnetic pressure as well as a decrease in hydrodynamic pressure. The results is that matter enters to this region from the sides and is expelled along the vertical directions as illustrated with the $y$ component of the velocity in Figure 15. In this case we use precisely the velocity as the function to be tracked by the refinement criterion, which as seen in the plots, follows the discontinuities of $v_{y}$.

Within the subject of Space Weather, the most commonly studied system is the Solar Wind, including the stationary wind.

In order to simulate this type of processes in the Solar System, we first solve the MHD equations with a null magnetic field in a cubic domain $[-x_{max},x_{max}]\times[-y_{max},y_{max}]\times[-z_{max},z_{max}]$, with $x_{max}=y_{max}=z_{max}$, centered at the position of the Sun, and consider $x_{max}$ bigger than the major semi-axis of Earth’s orbit. For the position of Earth we use Heliocentric Earth Ecliptic (HEE) coordinates, in which Earth’s position moves along the positive part of the $x-$axis from perihelion to aphelion.

Solar Wind simulations are carried out in a particular domain of interest, the volume contained from a sphere of $r_{in}\sim 20R_{\odot}$ (see e.g. (Riley & Gosling, 1998; González-Esparza et al., 2003; Feng et al., 2021)), until the faces of the cubic numerical domain. The reason to use a sphere of radius $r_{in}\sim 20R_{\odot}$ as an inner boundary is that at that distance from the Sun’s surface, the characteristic speeds of the wind are all pointing outwards, and consequently the flow is in the wind regime, as opposed to the accretion regime which is possible at shorter distances. This allows the injection of the physical properties of a wind or a storm through this inner boundary without the danger of material entering back in an accretion process.

The permanent injection of plasma at the inner boundary can be simulated by maintaining the physical variables of the Solar Wind $\rho_{in},{\bf v}_{in},p_{in},{\bf B}_{in}$, with suitable values at $r=r_{in}$, constant in time. Notice however that in Cartesian coordinates this inner boundary is the boundary of a lego-sphere, not an actual sphere, and therefore the implementation of boundary conditions at each face of boundary cubes is delicate (Kleimann et al., 2009). What we do to overcome this problem is to fill-in the whole lego-sphere with values of the primitive variables characterizing the wind at all times with the values at $r=r_{in}$, then the plasma will propagate outwards and eventually will fill the interplanetary space until it settles down to the heliosphere density and velocity profile (Guzmán & Mendoza-Mendoza, 2022).

At the external boundary, namely the faces of the cubic domain, we would like the fluid to continue its way out through the solar system and then use outflow boundary conditions at the outer cubic boundary.

As a particular example we simulate the formation of one of the stationary Solar Winds cases in (González-Esparza et al., 2003), which is spherically symmetric, with no magnetic field, in a domain with inner boundary $r_{in}=17.18R_{\odot}$, density $\rho_{in}$=$2100~{}\mathrm{cm}^{-3}$, velocity ${\bf v}=v^{r}\hat{r}$, $v^{r}_{in}=2.5\times 10^{5}$m/s, temperature $T_{in}=5\times 10^{5}$K and zero magnetic field and adiabatic index $1.4$. We use the numerical domain to be $x_{max}=250R_{\odot}$ for the wind to expand and cover at least a volume that contains Earth’s orbit.

Translation from code to physical units. Following (Guzmán & Mendoza-Mendoza, 2022), code and physical units are defined in terms of fixed scales of the state variables, for density $\rho_{phys}=\rho_{0}\rho$, pressure $p_{phys}=p_{0}p$, length $l_{phys}=l_{0}l$, time $t_{phys}=t_{0}t$, velocity $v_{phys}=v_{0}v$, temperature $T_{phys}=T_{0}T$. Notice that there is no restriction on density scale $\rho_{0}$, that we choose to be one.

We define code units by fixing the length scale to one solar radius $l_{0}=R_{\odot}=6.9634\times 10^{8}~{}{\rm m}$, time to one hour $t_{0}=3600~{}{\rm s}$, which defines the velocity scale $v_{0}=\frac{l_{0}}{t_{0}}=\frac{6.9634\times 10^{8}}{3600}~{}{\rm m/s}=1.93428\times 10^{5}{\rm m/s}$. Density scale is fixed to $\rho_{0}=\frac{m_{H}}{cm^{3}}=1.003573\times 10^{-21}{\rm kg/m^{3}}$ which defines the magnetic field scale $B_{0}=v_{0}\sqrt{\mu_{0}\rho_{0}}$. Temperature scale is $T_{0}=m_{H}v_{0}^{2}/k_{B}$ with $m_{H}=\mu m_{P}=0.6\cdot 1.6726219\times 10^{27}$kg the Hydrogen molecular mass and $k_{B}$ Boltzmann’s constant, and pressure and energy scales are $p_{0}=\rho_{0}v_{0}^{2}$ and $E_{0}=\rho_{0}v_{0}^{2}$.

Initial conditions are available given for the density, velocity and temperature of the wind as follows:

The simulation uses a base resolution of $80\times 80\times 80$ cells with 3 refinement levels, producing an equivalent resolution of $640\times 640\times 640$ cells. We use the refinement criterion (11) with a threshold value $\chi_{r}=10$. Outflow boundary conditions are set at the borders of the numerical domain. Numerical fluxes are computed using the HLLE formula along with the mimod reconstructor. A CFL factor of $0.125$ was used throughout the whole simulation.

At initial time the wind variables have constant values within the sphere of radius $r_{in}$, which simulates a stationary pumping of plasma into the heliosphere. The gas fills the numerical domain with the front wave that propagates outwards until it reaches the external boundary and leaves the domain. In Figure 16 we show how the refinement criteria follows this front wave. Once the wave leaves the boundary the flow approaches a stationary state.

The result of this simulation appears in Figure 17, where we show snapshots of density, $x$ component of ${\bf v}$ and temperature of the gas along the $x-axis$, until the flow stabilizes in the whole domain by $t\sim 1000$ hours. In order to have a reference of the spatial domain, in these coordinates (HEE) the average Earth location is $x_{Earth}=214.83452R_{\odot}\hat{x}$.

After the wind has become stationary, we inject a toy-model of Coronal Mass Ejection (CME). We characterize the CME launched from the inner boundary at $r=20R_{\odot}$ by its density, velocity and widening angle $\theta=\pi/6$. The injection is implemented similarly as the stationary Solar Wind at the inner boundary provided density, velocity and temperature of the CME. As a test we inject a toy CME with the following properties $\rho_{\mathrm{CME}}=3\times\rho_{\mathrm{{wind}}}$ $v_{r,\mathrm{CME}}=2.5\times v_{r,\mathrm{wind}}$ and $T_{\rm CME}=2\times T_{\mathrm{wind}}$, conditions taken from (González-Esparza et al., 2003). We launch this CME centered at the solar equator, during a time window of 3 hours. In Figure 18 we show snapshots of density during the propagation of the CME. With this initial setup the initial burst propagates, its density profile smears out and widens. In Figure 19 we show that the structure of the highest refinement level tracks the perturbation.

Beyond the solar corona, the solar wind is fairly well and commonly modelled using ideal MHD. In this context, the macroscopic plasma quantities, such as density, velocity and temperature, can all be related to the topology of the field lines of the magnetic field. Current techniques applied to simulate the interplanetary heliosphere rely on this fact, for example the WSA (Arge et al., 2004) approximation and its implementation in the ENLIL model (Odstrcil, 2003).

In order to explain our implementation of the magnetic field we briefly review the bases of empirical models. The physical domain is partitioned into some basic spherical shells approximately located at given spheres: the photosphere at 1 $\mathrm{R}_{\odot}$, the Solar Corona (SC) at $\mathrm{R}_{SC}\approx 2.5\mathrm{R}_{\odot}$, the inner boundary at $R_{in}\approx 20\mathrm{R}_{\odot}$, and the rest of the numerical domain that includes a volume containing at least the Earth’s orbit, of $\approx 1\mathrm{AU}$. The magnetic field at the photosphere is reconstructed from magnetogram synoptic maps, used as input to construct analytic expressions of the magnetic field at the Corona via the potential field source surface method (Newkirk & Altschuler, 1969). Next, the super-radial expansion factor ($f_{s}$), that refers to the variation of field strength along a magnetic flux tube is calculated, for example as done in (Shiota et al., 2014):

Once $f_{s}$ is calculated, we use the empirical Wang and Sheely formula (Wang & Sheeley, 1990; Arge & Pizzo, 2000) that gives the wind radial velocity at the inner boundary as

Based on the Parker Solar Wind model (Priest, 2014), we use a constant temperature, as well as a magnetic field that fulfills conservation across radial flux tubes: