The code solves the mass, momentum and energy conservation equations for a compressible fluid numerically using a finite difference scheme. The model is one-dimensional, since the plasma is confined in a semicircular loop where plasma motion and heat flow are allowed only along the magnetic field lines. Symmetry conditions allow us to consider only one half of the loop:

where *n* is the hydrogen number density; *t* is the time, *s* is the field
line coordinate; *v* is the plasma velocity; is the mass of hydrogen
atom; *p* is the pressure; *g* is the component of gravity parallel to field
line; is the viscosity; is the ionization fraction
where is the electron density; *T* is the temperature; is
the thermal conductivity ( K; is the Boltzmann constant; is the hydrogen
ionization potential; are the radiative losses per unit
emission measure (Raymond et al. 1977); is the volumetric
power input to the solar atmosphere:

is the heating term, uniform in the corona and tabulated in the chromosphere, required to balance radiative and conductive losses for a static atmosphere; the second term is the impulsive heating triggering the flare.

The numerical scheme has been described previously by Peres et al. (1982) and Peres & Serio (1984). The mass and momentum equations are solved with an explicit numerical scheme, with upwind differencing schemes used for the advective terms (Roache 1976). The energy equation is solved with an implicit scheme, which prevents numerical instabilities (Richtmyer & Morton 1967; Fletcher 1988). Initial profiles of temperature and pressure are derived following the hydrostatic model of Serio et al. (1981), and matched with one of the chromospheric models of Vernazza et al. (1981).

As Dwyer (1983) pointed out, one of the key problems of numerical fluid dynamics is the simultaneous calculation of both the grid and the solution; in particular Dwyer suggests that the generation of the grid is perhaps more important than the approximation used for the spatial differencing.

The PH code has been designed on Eulerian coordinates. Lagrangian coordinates are appropriate to grant automatically a good resolution of the continuity equation but do not offer a significant advantage in the resolution of the fast moving steep transition region, one of the most critical aspects of the coronal loop dynamics and thermodynamics.

The use of a re-adaptive Eulerian grid generation method is probably the fittest choice to solve differential equations whose solutions present very steep gradients rapidly varying in time.

Another fundamental requirement of any finite-difference numerical grid scheme is the regularity of the spacing, i.e. the gradual change of the spacing along the grid (Thompson et al. 1985; Hawken et al. 1991; Tamamidis & Assanis 1991). Lack of regularity, namely when the cell size varies abruptly along the grid, typically leads to numerical errors even worse than those caused by inadequate spatial grid resolution, as many authors have pointed out (Fisher et al. 1985a,b; Thompson et al. 1985; Gan et al. 1991; Hawken et al. 1991; Tamamidis & Assanis 1991) and we ourselves have found.

In our approach we have essentially maintained the numerical scheme originally devised by Peres et al. (1982) for a time-independent spatial grid (with step changing according to a geometric progression of which the uniform grid is a special case). We have reorganized the code to account for the recurrent change of grid and succeeded in satisfying both the requirements of small relative change of the physical variables and of the cell sizes along the grid, with the twofold advantage of keeping both the number of points and the computing time low.

As a first fundamental condition the grid must satisfy the requirement that the fractional variation of temperature, density and pressure between two adjacent points remains below a fixed value ,

in the transition region and in the corona. Since this condition must be
satisfied throughout the evolution of the solutions, the grid has to be modified
according to Eq. (8), as discussed in the following, and the variables
must be linearly interpolated on the new grid points whenever required.

**Figure 1:** Profiles of temperature vs. field line coordinates along one
half of a symmetric loop of total length and base
pressure (solid line). The dots show the cell size
required to resolve spatially the temperature with a fractional jump of
exactly 10%. The dashed line yields the cell size derived according to our
algorithm. The inset shows an enlargement of the base of the transition
region

Figure 1 (click here) shows a typical temperature profile vs. field coordinate (solid line), computed according to the static loop model of Serio et al. (1981), based on hydrostatic equilibrium and energy balance, as well as the irregularly distributed cell sizes that would be required to yield exactly a 10% jump of temperature between any two adjacent points (dots), and the cell sizes determined with our algorithm, discussed below in more details (dashed line). The loop shown is a typical active region one with half-length , pressure at the base of the transition region 6 dynes and maximum plasma temperature at the top .

We have ascertained the importance of grid regularity with a simple test: if the initial hydrostatic atmosphere is used as initial condition for the numerical code, and , we expect no evolution apart from some numerical noise; and in fact so was the case for the logarithmic grid fixed in time with the previous PH code (Peres et al. 1982), where the hydrostatic profiles remain essentially unchanged for integration times much longer than the characteristic dynamic and cooling times (i.e. longer than thousands of seconds). The non-uniform grid with spacing as in Fig. 1 (click here) (dots) instead yields significant changes in the profiles of temperature and density, because the grid is too irregular.

We have therefore made our grid more regular; in particular we have generated a grid finer in the transition region, where the highest resolution is needed, and whose spacing increases smoothly on both sides. The grid is split into three sections. In the first, the spacing decreases along field coordinates according to a geometric progression from the loop base to the transition region. The second is a small uniform grid (generally ten points), of high resolution placed just around the location where the smallest cell is required, typically at the base of the transition region (see the inset of Fig. 1). Finally, in the third section, the spacing increases according to another geometric progression from the base of the transition region to the loop apex.

The whole grid is shown in Fig. 1 (click here) (dashed line). The two geometric
progressions on the sides of the region with uniform spacing (in the
first and in the third section) may have different ratios. The size of
the *k*-th cell on either size of the uniform grid
is determined as , .
The code selects the values *q* of the two ratios as the largest
enforcing
condition (8).
We dedicate particular care to the choice of the
cells adjacent to the uniform grid, in order to make the size of the cells
vary as gradually as possible; the first cell of each of
the two
geometric progressions satisfies the condition:

and is determined so that the two progressions span exactly the
remaining distance on each side of the loop. The initial temperature, pressure
and
velocity profiles are then linearly interpolated from the old grid onto the
new grid.

After each step of time integration the algorithm checks the grid and
verifies if the point
where the smallest cell size is required by the condition (8) is still inside
the section with
uniform spacing, otherwise it generates a new grid with the same criterion
described above. Since the value of the smallest size , its
location along the loop, and the geometric ratio *q* change frequently
during the flare evolution, the number of points used along the loop is
not constant but varies according to the requirements imposed by the conditions
(8) and (9).

The small uniform grid at the base of the transition region has two advantages:

- i)
- the computation of a new grid is not required at every time step;
we therefore reduce the computing times and limit the errors inevitably
introduced by interpolation.
- ii)
- it gradually matches two sections with different spacing,
reducing numerical errors significantly, as discussed above.

As for the influence of truncation errors in the new implementation, we focus on the energy Eq. (3) and, in particular, on the dominant thermal conduction term. The scheme has been originally devised to work on a geometric progression grid. In the section of the grid with uniform spacing (the one with the steepest temperature gradient) the error is of the order of:

In the non-uniform grid , apart from a factor of the order of unity. In the
region of the small uniform grid, is very small, but, since
the temperature derivatives are very big, the error is *de facto*
dominated by the factor in parentheses. Far outside that region,
is large but the temperature structure is much smoother,
almost flat; since the critical part of our model is anyway the
transition region, we take care that condition (8) is always
satisfied.

The adoption of a new grid has significant consequences also on the value of
the integration time step . In particular, is selected
as the minimum value among the following times evaluated for each *i*-th
point in the grid:

- the
*Courant-Friedrichs-Lewy*(CFL) time

where is the distance between two grid points and is the sound speed

- the
*radiative*time

- the
*conductive*time

From (10) (see, for example, Richtmyer & Morton 1967; Fletcher 1988) it is immediately evident that a reduction of leads to a reduction of The conductive time is indeed but we have seen that during most of the flare simulations the time step is in fact governed by the CFL condition.

Occasionally the time step decreases below s, however the code continuously monitors and increases it whenever possible. Indeed if the time step remained below s, we would typically need hours of CPU time to simulate 600 s of a flare on a DEC Alpha 3000, since an integration time step takes s (with 250 grid points along the loop). With this strategy computations take approximately one half this time. Our strategy is even more advantageous for runs covering the decay phase, where the time step can increase markedly.

We can also limit the time step imposing a lower size of the grid spacing , typically set at 100 cm, with the minor disadvantage of a resolution at the base of the transition region sometimes slightly less than adequate for Eq. (8).

Following the suggestions by Thompson et al. (1985), we also make sure
that the time step and, independently, the minimum space step may change by, at most, a factor *q*.
We have seen that such a condition helps in further reducing
the numerical fluctuations sometimes present.

web@ed-phys.fr