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:
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:
where 
is the distance between two grid points and 
is
the sound speed
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.