B. Appendix B

B.1. Sets of variables without central singularities

  With the lagrangian independent variable , and for the variables and in place of R and L, one seeks exponents , and in such a way that the gradients of pressure, radius and luminosity are O(1) in a neighborhood of the center, M=0. The gradients are written:

In the vicinity of the center the density, , is , therefore R=O(M^1/3) and ; similarly, L=O(M) and and, at the limit :

The singularities at M=0 are avoided if the exponents p, r and l are non negative integers; in Table 3 (click here) the exponents , and , deduced from the previous set of equations, are given for the values 0, 1 and 2, for p, r and l.

   

p r=0 r=1 r=2 l=0 l=1 l=2

0 2/3 2 4 6 2/3 4/3 2

1 1/3 1 2 3 1/3 2/3 1

2 2/9 2/3 4/3 2 2/9 4/9 2/3

From Table 3 (click here), with the exponents , and , i.e., with the variables M^2/3, R² and L^2/3, one ensures p=0, r=0 and l=0, therefore , and as expected.

Obviously other sets of exponents can be employed, the best of them is chosen according to the following numerical criterion: working with a scheme of, let say, third order, between two grid points the dependent variables are quadratic polynomials in ; any of them, for example, is written: . In a neighbourhood of , the gradient are ; with p=0 one has:

while, if p=1:

The last approximation only uses the term in , therefore information is lost. Hence the set of exponents , and , i.e., the variables M^2/3, R² and L^2/3, are optimal.

   

p r=0 r=1 r=2 l=0 l=1 l=2

0 2 2/3 4/3 2 2/3 4/3 2

1 1 1/3 2/3 1 1/3 2/3 1

2 2/3 2/9 4/9 2/3 2/9 4/9 2/3

From table 4 (click here), where the exponents for eulerian variables are given, one sees that the highest accuracy cannot be obtained with the set of variables of Eq. (26 (click here)).

B.2. Setting a grid point at a given location.

  At each LMR the derivative of the temperature gradient is a discontinuous function. The numerical technique used for the integration of the equations of internal structure is based on an approximation of the unknown functions by piecewise polynomials (see Appendix B.5 (click here)). Between two mesh points, each polynomial piece is a differentiable function, therefore their derivatives are allowed to be discontinuous only at grid points where the polynomial pieces are connected; as a consequence, each LMR must be localized precisely on a knot; this allows, with an ad hoc choice of the basis, to have different values on left and right for the derivatives of the unknown functions. For sake of clarity in the discussion and the notation, one assumes that there is only one LMR where a discontinuity arises (labeled hereafter with the subscript "d''); the extension to several discontinuities is straightforward. The location of each discontinuity on a knot is obtained by a slight change of the algorithm used for the automatic mesh refinement; it is, basically, an inverse linear interpolation: The grid points are localized in such a way that, from a mesh to the next, the change of the distribution function is constant; inside each mesh, indexed by i, the relative weights of the variables are defined by the distribution factors and known within a factor , , which was, up to now, implicitly assumed to be the unity for each mesh, but it could be adjusted (iteratively) in order to set each LMR on the closest grid point; in other words, locally, for both sides of each LMR, instead of using the constant value C(t) the repartition constant is adjusted to the amounts just needed for fitting a grid point right on the LMR. Doing that, instead of Q(q), for , one uses a "weighted distribution function", defined by pieces, as:

 
Analogously one also defines:

The pieces of the function being connected by continuity: , . Therefore, instead of Eq. (20 (click here)), one uses:

 
Note that and are step functions. Due the facts (see Appendix B.5 (click here)), that at grid points first, the differential equations are not written and, second, discontinuous derivatives are allowed for, the discontinuities of are implicitly managed. Let, at time t, a solution of Eq. (21 (click here)) be known and Q_i⁽⁰⁾(q) be the weighted distribution function defined by Eq. (A1 (click here)) with . Here each quantity with a "(j)'' superscript is related to the j-th iterated model, At the location , where the discontinuity occurs, the distribution function, calculated with Eq. (22 (click here)), amounts to ; let l be the index of the nearest mesh point from the discontinuity () (in this discussion one simplifies by assuming that the index l does not change as it does, sometimes, in the calculations); since Q^(j)(q) is a strictly monotonous function, there is only a unique value such that , therefore, for , , the weights can be defined by: and, for :

Hence, at the iteration , the LMR will be closer, than previously, to the grid point with the index l. After convergence (), , . The jump of the distribution function, between q_l-1 and will satisfy:

just that is needed. For j=0 the weights are initialized to unity and updated between each Newton-Raphson iteration of the quasi-static equilibrium problem with Q^(j)(q) and q^(j)_D derived from the locations of the LMR at the (j-1)-th iteration, for j>0. Figure 7 (click here) shows a typical displacement of a grid point toward a LMR (see also Fig. 5 (click here)). In this simplified case only three iterations have been necessary to set a grid point on the LMR.

  
Figure 7: Automatic adjustment of a grid point on the limit between a radiative and a convective zone. Two iterations have been needed to push a grid point towards the LMR; during the process the radiative gradient slightly decreases and the LMR moves right. The adiabatic gradient , full line, remains unchanged. The initial locations of the grid points (empty squares), are moved to the right at the first iteration (empty circles) and slightly right again at the second on (full circles)

This iterative process is only of first order. It will converge only if the functions considered are sufficiently smooth and, moreover, if the LMR is sufficiently well defined, as in the case illustrated in Fig. 7 (click here); it is the case for the external LMR of the Sun and for the convective core of stars of intermediate mass. On the other hand, for semi-convective zones or convective cores with diffuse LMR, as for the convective core of the young Sun, the closest index l, from the LMR, is not well fixed and the algorithm may not converge; fortunately, in such cases, the jump of the slope of the gradient is almost zero, therefore a mislabelling of a small jump in the slope of the gradient has almost no consequence.

B.3. Precise external boundary conditions.

  As described Section 2.4.3, the precise reconstitution of an atmosphere consists in a differential problem, with boundary conditions at three different levels: , and ; recall that, at , defined by Eq. (10 (click here)), there is a open inner limit. In the numerical solution of Eq. (11 (click here)), one assigns, among the grid points allowed for the atmosphere, a fixed index to the moving limit ; that is done by the piecewise logarithmic linear mapping: of in defined by:
 
Here the slopes and are defined by:

hence one has:

, at the bottom of the atmosphere where ,

, at the inner limit where ,

, at the upper limit of the atmosphere where ;

the index is a open parameter; with , will lie in the vicinity of the middle of the interval . From this mapping results a constant mesh size, it is an advantage of the algorithm. For the numerical solution it is convenient to add one more unknown: the temperature T, which fulfills the scalar equation , and the boundary condition . With the following variables, employed for the reconstitution of the atmosphere: , , , , and , Eq. (11 (click here)) are written:

 
where (respt. ) if (respt. ). The boundary conditions are:

here:

The numerical solution of the non-linear differential problem Eq. (A4 (click here)), with its boundary conditions at three different levels, is solved using the B-spline collocation method (see Appendix B.5 (click here)). The iterative process is initialized either, with an approximate model of the solar atmosphere for PMS or ZAMS initial models or, along the evolution, with the atmosphere of the previous model . The solution gives , , and their derivatives with respect to and . Typically, one uses from to grid points for the interval . With B-splines of order 3, there are from 59 to 99 points where Eqs. (A4 (click here)) are fulfilled; due to the superconvergence (see Sect. B.5 (click here)), the order of the scheme is 4.

   
Table 5: Coefficients in Eq. (A5 (click here)) for the IRK formulas Lobatto IIIC (see text) of order p=1,2,4 with s=1,2,3 steps. As usual, the are in the first column, the b_i in the last row, the remainder entries are the coefficient a_ij of the IRK matrix

B.4. Stable integration with IRK formulas.

  With stiff problems, care must be taken for saving the stability; the procedures are briefly recalled here in the case of only one chemical element; let y₀ (respt. y) be the abundance of that element at time t (respt. ). With the coefficients reproduced in Table 5 (click here), the IRK scheme is written:
 
where T_j, and y_j are, respectively, the temperature, density and abundance at intermediate times , ; s is the number of steps of the IRK formulae. This standard formulation for IRK schemes does not give satisfactory results for the evolution of the chemical composition: negative values for the abundances are observed; that is due to the amplification of the errors by large Lipschitz constants. The round-off errors are reduced if the first Eq. (A5 (click here)) is rewritten for the smaller quantities z_i=y_i-y₀ (Hairer & Wanner, loc. cit., Sect. IV.8) as:
 
The matrix of any IRK Lobatto IIIC formulas is invertible therefore, with , one has d_j=0 if and d_s=1; then the second Eq. (A5 (click here)) is simply written:
 

B.4.1. Instantaneous mixing with IRK formulas.

In MZ the abundances have constant values as described Sect. 3.3.2 (click here); for the chemical species y, let be this mean. At time , the mean, , is approximated using the mid-point rule:

here the sum involves all the meshes inside the MZ; y_0,k+1/2 is the abundance at the mid-point, and:
 
Therefore Eq. (A6 (click here)) is written:
 
where (respt. ) is the temperature (respt. the density) at the intermediate point at time ; then is written:
 
The extension to several chemical elements of Eq. (A6 (click here)) to Eq. (A10 (click here)) is straightforward. Note that, with the use of the mid-point rule, Eq. (A8 (click here)) for the numerical approximation the integral of Eq. (28 (click here)) over the MZ, is only of second order; a two-points Gaussian formula (i.e. third order accuracy) will be more appropriate. The error is not important in the solar case since the thermonuclear reactions are not really efficient in MZ.

In CESAM, Eq. (A6 (click here)) to Eq. (A10 (click here)) are solved iteratively using Newton-Raphson scheme with an analytical jacobian updated at each iteration; the initial values of the z_i are fixed to zero; the quadratic iterative process is halted as soon as the corrections on the z_j are less than 10^-8 in relative values; most of the time, the convergence is reached in less than five iterations.

B.4.2. Conservation of nucleus and charge numbers.

  For a number of chemicals , let at time t:

be any linear relationship between the abundances with fixed coefficient , e.g. conservation of baryon and charge numbers are relationships of that form, then:

Owing to the linearity of Eq. (A6 (click here)) and Eq. (A9 (click here)), for elements, with obvious notations:

and:

therefore, with IRK formulae, the number of baryons and the number of charges are exactly conserved up to the accuracy of the closure of the Newton-Raphson scheme.

B.5. Integration of the diffusion equation

  One seeks the X_i as piecewise polynomials of order with discontinuous first derivative at each LMR, i.e., . For the calculations it is convenient to project them on their normalized B-spline basis, (see App. B (click here)), , is the dimension of the basis; therefore is written:
 
One looks for the weak solutions of Eq. (30 (click here)) using the Petrov-Galerkin's formalism (Schatzman et al. ; Marchouk & Agochkov, ; Quarteroni & Valli, ): Let be the B-spline basis without discontinuity, i.e., , constructed on the grid of - schematic representations of nodal vectors of and can be found App. B.5.2 (click here); for , multiply both sides of Eq. (30 (click here)) with , and integrate over the whole star:
 
with the standard notation for the inner product of two functions f and g. Due to the continuity of the fluxes F_i, the standard integration by parts gives:
 
Owing to the boundary conditions, at the center , at the outer limit of the envelope and to the continuity of X_i and F_i, the set of Eq. (A12 (click here)) is written:
 
here all the quantities are taken at the time except X_i^(t) taken at time t. This fully implicit scheme is of order in space and of first order in time.

One emphasizes the following three facts: (i) in Eq. (A12 (click here)) the discontinuities of the first derivatives of the X_i are implicitly taken into account by the integral formulation, (ii) the use of the B-spline basis is convenient because the algorithms remain simple even if the locations of the LMR change with time, (iii) the Petrov-Galerkin formulation only involves first order derivatives, (iv) the fully implicit scheme has no restrictive time step condition (Richmyer & Morton 1967).