[mesa-users] accretion

[Bill Paxton]

Hi Dean,

Consider this fairy tale.  We don't need to worry about a CFL limit for advection here because we aren't actually doing a time integration of the time derivative that includes the advection term.  For the near-surface cases where df/dt at constant relative mass is very tiny, then the time derivative will be dominated by the spatial term.  That means we are using values from different locations at the same time to form an "instantaneous" value of the time derivative to calculate an instantaneous eps_grav for use in the instantaneous energy balance equation (dL/dm = ...).  Hurray!   no CFL timestep worries after-all because we aren't doing any time integration of the Ds/Dt term; we're just using an instantaneous value of it to check instantaneous self-consistency for energy conservation.  This is completely different than the time integration of convection and the CFL limits it must respect.  So I should stop worrying about advection terms in the eps_grav calculation, right?

Well, maybe I should still worry a little.  When the dt/dt term isn't tiny, then it is approximated by an average value by a simple 1st order difference.  If we don't have an advection term, then we are consistently using average time derivatives since (f_final - f_init)/(t_final - f_init) is the average rate of change and is not necessarily the instantaneous rate of change at the final time.  When we split the Lagrangian into a time part and a space part, we necessarily introduce a mix of average and instantaneous.  But maybe that's okay.  Maybe it is even good.  But we should keep it in mind; the Lagrangian time derivatives in the split scheme will be instantaneous values near the surface (where df/dt at constant q is tiny and the spatial term dominates) and will change to average values in places where df/dt|q begins to dominate.

There is also still a potential concern about increased numerical round off errors introduced by the spatial derivatives.   In the Sugimoto&Nomoto form, we use a spatial derivative of entropy wrt relative mass which will be approximated by (S(k-1) - S(k))/dq_bar(k) where S(k-1)is the entropy at the center of cell k-1, S(k) is at the center of k, and dq_bar(k) is the relative mass difference between the cell centers.  We've seen that dq's can get very small so the difference between the S's at adjacent cells will be correspondingly small, so the difference may suffer from round off.  In Townsley&Bildsten you replace the dS/dq by a difference between adiabatic and actual temperature gradients, both measured at the same location and hence not dependent on the size of the cell or round off from subtraction of nearly equal values.  Excellent!  But do we need to enforce that this is only used in radiative regions?  Out near the surface we can get superadiabatic convection regions -- will they screw up the Townsley&Bildsten scheme?

We'll also need to take care to minimize round off problems in calculating the time difference part of the Lagrangian time derivative.  If we are doing differences in values at same relative mass then we should try to make the cell have the same relative mass coordinates at the end of the step as at the start.  Similarly if we are doing time difference at constant mass, then we should try to make the cell have the same mass coordinates at start and end.  If we can manage that, then we can arrange to have the newton solver directly give us the difference df for the cell without needing to do the subtraction df = f_final - f_init.   That works when f is one of the "basic" variables being solved directly in the newton iterations -- e.g., lnT and lnd (or lnPgas) are basic variables so we can do this trick, but entropy is not basic in that way -- it comes from the eos as a function of lnT and lnd.  But we can get Ds/Dt for eps_grav and still avoid the roundoff by rewriting the time derivative of entropy in eps_grav in terms of the time derivatives of lnT and lnd: Ds/Dt = cp*((1-grada*chiT)*DlnT/Dt - grada*chiRho*Dlnd/Dt).    Note that we can do this if the cell has either constant mass coordinates or constant relative mass coordinates.  The mesa adjust_mass routine arranges it so that cells near the surface have constant relative mass coordinates and ones in the core have constant mass coordinates, with a typically small transition region where the coords are changing in both absolute and relative mass.    Perhaps we'll be okay if we limit the round off problems to that transition region.  

Enough for now.....

-Bill