The endeavor of a numerical code of stellar evolution is, more or less, directed towards a given goal; either the purpose is the computation of models evolved for long time intervals as for the study of AGB stages (Han et al. 1994), in which cases the search of efficient and low-cost algorithms is the priority, or a high level of accuracy has to be reached, as for solar models which need high numerical accuracy in addition to precise physics (Christensen-Dalsgaard 1982; Reiter et al. 1994).
A first way for improving the numerical precision of a stellar model is the refinement of the mesh, a second is the use of a numerical algorithm of high order. Until recently the codes used for stellar evolution were, more or less, derived from the Henyey's code (Henyey et al. 1959, 1964) which uses finite differences and so, have a fixed order of accuracy, e.g. first or second; with this class of schemes, the only way to have various orders of accuracy is to write several algorithms.
For the integration of two-point boundary value problems, high accuracy
(typically 10-10 or even better) can be obtained with the multishooting
technique (Stoer & Bulirsch 1979, Sect. 7.3.5) employed in few
codes of stellar evolution (Wilson 1981; Reiter et al. 1994). They
are CPU consuming as reported by Embarek (1989) but they are well designed
for parallel computers (Reiter et al. loc. cit.). Solar models with such
impressive
numerical accuracy have been calculated with a very few number (
) of
fitting points. These methods, however, suffer of a conceptual gap between the
reachable accuracy in space and in time since the mixing of the
chemical species in Convective Zones (CZ) with such high-accuracy
schemes is, in practice, only possible if neither thermonuclear reactions nor
diffusion take place in CZ; therefore they are well suited for the calculation
of
Standard Solar Models (SSM).
Applications of solar models are
the computation of the eigenfrequencies; that needs a precise knowledge of the
solution around the nodes
whose locations change according to the modes, their
orders and degrees; therefore the solution needs to be restored at locations
well
distributed around the turning points, i.e., between the grid points, that needs
lengthy calculations with multishooting methods if the full accuracy
is preserved.
The goal of CESAM
is to perform a standard solar evolution from ZAMS to
present age with an accuracy of the order of 10-4, using
about 500 mass shells and 50 time steps (with a maximum time step of 100 million
years); according to the physics employed it
also calculates evolution from PMS to the onset of 4He
burning for various stellar masses. The method of integration does not use the
familiar finite differences; it is based
on the so-called splines collocation method derived by deBoor (1978,
Chapt. XVII). The basic idea is to seek piecewise polynomials as
solutions of differential problems; among all the bases of the linear space of
piecewise polynomials, the B-splines basis (deBoor, loc. cit.; Schumaker
1981) has been designed especially for calculation purposes. Its main
advantages are: (i) the order of the scheme has only one parameter dependence,
(ii) algorithms specially derived for spline calculations allow efficient and
stable computations, (iii) a convenient method has been designed for
connecting the pieces of polynomials with a given order of continuity,
(iv) only few changes in the algorithms allow
to get either the weak solutions (Galerkin's methods) or the
strong solutions (collocation method),
(v) ability of restoring exactly the numerical solution anywhere, (vi)
ability to handle discontinuous solution, (vii) no need
of a peculiar algorithm at the center.
However, due to the non-trivial and
unfamiliar algebra of B-splines, the algorithms are much more intricate than
with the standard finite differences.
With the use of a functional basis (i.e. presently, the B-spline basis) it is
a triviality to structure
the code in a "numerical space" where the differential problem is formally
solved and in a "physical space" where the equations are formed whatever the
numerical method used for their solution is; hence all the physics: Equation Of
State (EOS), opacities, nuclear reaction rates, external boundary conditions,
calculation of the convective flux, etc... are calculated in external
routines which can be, on need, user supplied or taken from the source itself
without any change in the "numerical space". Hence CESAM allows calculations of
stellar models with various physical assumptions, physical data, external
boundary conditions, numerical methods and numerical accuracy; it also includes
diffusion and mass loss.
For consistency the same order of accuracy in time and in space is desirable; with respect to time, the fourth order is a hardly reachable upper limit, that is due to (i) the stiffness of the problem: the evolutionary time scales of the various chemical species differ by more than eighteen orders of magnitude, (ii) the convective mixing can hardly be done with a high order scheme. For standard models, i.e., without diffusion, the numerical scheme currently employed is the first order Euler's backward formulae (Kippenhahn et al. 1968), with simplifications (Arnett & Truran 1969) or improvements (Wagoner 1969), (Christensen-Dalsgaard 1982) makes use of the second order trapezoidal rule; CESAM uses a stiffly stable Implicit Runge Kutta (IRK) scheme with order up to 4. With diffusion, Runge Kutta scheme (Alecian 1986), finite difference schemes (Cox et al. 1989), Crank Nicholson (Proffit & Michaud 1991; Charbonel et al. 1992) are employed; in CESAM the diffusion and thermonuclear equations, with mixing of chemicals in CZ, are solved as a whole using the implicit backward first order Euler's scheme in time and, in space, the Petrov-Galerkin's method with order up to four.
Some recent papers are based on calculations made with CESAM (Morel et al. 1990; Morel et al. 1993, 1994; Berthomieu et al. 1993; Goupil et al. 1993; Audard & Provost 1994; Lebreton et al. 1995; Chmielewski et al. 1995; Morel & Schatzman, 1996). Preliminary and short descriptions of the numerical methods have already been given (Morel 1989; Morel 1993a, 1993b) or referenced (Gabriel 1990); the algorithms have been adapted for the computation of evolution of giant planets (Guillot & Morel 1995). This paper is devoted to a detailed description of the numeric of CESAM.
In Sect. 4 the basic physics of stellar evolution is briefly recalled for references and notations; special emphasis is given to the mixing of CZ, on initial conditions and on external boundary conditions for which a detailed restoration of an atmosphere is allowed for. Sect. 3 (click here) is devoted to the numerical techniques; a first part is concerned with the solution of the two points boundary problem with emphasis on the choice of the integration variables and on the automatic location of grid points; in a second part, the schemes used for solving the stiff differential initial value problem of the evolution of the chemical composition with, and without, diffusion are described. In Sect. 4 (click here) some indications concerning the management of the code are given. Results of investigations on internal accuracy are presented in Sect. 5 (click here) and the discussion in Sect. 6 (click here). Three appendices are devoted to numerical technical discussions.