[mesa-users] mesa release 3844

[Bill P. <pa...@ki...>]

Hello and Happy New Year!


The first mesa release of 2012 has some cool new features:

	-- pressure instead of density as a primary variable in star
	-- Paczynski grey atmospheres for surface BCs and in pulsation models
	-- star/test_suite/astero example
	-- Hooke-Jeeves minimization option for astero model searches

Please read on if you'd like details!

Cheers,
Bill




pressure instead of density as a primary variable in star

Ideally it wouldn't matter whether you used density or pressure as a primary variable
in the stellar structure solutions -- for a given temperature and composition you always
can go back and forth from a density to a gas pressure, so it seems you can pick which you
use for reasons such as performance.  Since the microphysics packages tend to use
density, we use have used it in star as well.  In other words, the solver for star treats
density as a primary variable and pressure is indirectly determined as a function of
density, temperature, and composition.   If the equations were solved precisely,
then this extra level in the calculation of pressure would be harmless.  But the equations
are only solved approximately, so the indirection means that approximately correct values
for density and temperature can lead to anomalous pressure profiles -- i.e., cases
that show (tiny) violations of the expected hydrostatic balance.  The violations only
seem to show up near extreme jumps in density (such as at a sharp H/He boundary),
and the size of the violations seem to be negligible, so previously we haven't done anything 
about them.  But now the asteroseimologists have entered the scene!  They'd like to
have a perfect pressure profile as a starting model for their pulsation codes, and the
violations of hydrostatic balance, however small, are a cause for concern.

So now we have the option to replace density by gas pressure as a primary variable.
The good news it that it does seem to fix the problems with the pressure profiles.
Here's an example: this is a 25 Msun star crossing the gap at a time when Teff = 1e4.
There is no overshooting or semiconvection, so it has an extremely sharp H/He boundary;
in fact, we have single cell step function transitions in both abundances and in density;
as expected, logT is smooth across the transition.

the transition region is very thin in both mass and radius



pressure near H/He boundary is smooth when using lnPgas as primary variable (left);
however pressure is ragged when using lnd as primary variable instead of lnPgas (right).


note that the magnitude of the irregularities in P are tiny -- they are probably negligible for most applications.
but for asteroseismology, they might matter -- I don't know yet if the smoother pressure profile actually helps.
Please let me know if you find evidence that it does!

BTW: thanks to Aldo Serenelli for alerting me to the potential importance of using
pressure in place of density as a primary variable.   I fought his suggestion at first,
perhaps because I unconsciously anticipated the hard work that was going to be
required to make the necessary changes to mesa!

in star/public/run_star_defaults

   ! lnPgas variables
      
      ! lnPgas_flag is true if we are using lnPgas variables in place of lnd
      change_lnPgas_flag = .true.
      new_lnPgas_flag = .true.


Note that the default is to use gas pressure in place of density.
To go back to using density, edit your inlist to say
      change_lnPgas_flag = .true.
      new_lnPgas_flag = .false.





Paczynski grey atmospheres

Certain asteroseismology cases require a detailed stellar model out to very small optical depths (tau = 10^-4 or so).
There are even cases with small but significant convection zones out in those optically thin regions (e.g., white dwarfs).
Now, with prompting from Mike MacDonald and Hideyuki Saio, mesa can supply these models by using a scheme
described by Paczynski in this 1969 paper:


Key features are that this takes into account radiation dilution when tau < 2/3,
and calls the MLT package to get gradT allowing for convection.


To use it as the atmosphere boundary condition in star, set 
         which_atm_option = 'Paczynski_grey'   
               ! create an atmosphere for given base conditions.
               ! inspired by B. Paczynski, 1969, Acta Astr., vol. 19, 1.

The atmosphere integration is repeated until it gives a matching radius;
the tolerance for the match is given by this parameter:
            Paczynski_atm_R_surf_errtol = 3d-4

The steps for the integration are in units of log tau.  This sets the max step size:
         create_atm_max_step_size = 0.1 ! in units of log10_tau

That takes care of using Paczynski grey atmospheres for the surface boundary condition.
Next step is to include the details of the atmosphere in the input stellar model
for the pulsation code.  To do that, set one of these (from run_star_defaults):
   add_atmosphere_for_fgong = .false.
   add_atmosphere_for_osc = .false.
Alternatively, to get pulsation output created whenever a profile is written,
set these (from star_defaults)
         write_pulse_info_with_profile = .false. ! if true, write pulse info file when write profile
         pulse_info_format = 'FGONG'
         add_atmosphere_to_pulse_info = .false. 

Here are some plots showing a model output including a Paczynski atmosphere
This is the same model as shown above: a 25 Msun star crossing the gap at a time when Teff = 1e4.
The interior-atm boundary for this case is set at the photosphere (tau = 2/3, log10(tau) = -0.176).
Note that the profiles are smooth out to tau = 10^-4 and density of 10^-12







star/test_suite/astero example


We are making progress with getting mesa in shape to do asteroseismology.
The recent program at KITP had a big impact!  I've added an example case to
the star test_suite showing a typical setup for using mesa to search for a
model that matches given data from observations.  Take a look here:

	mesa/star/test_suite/astero

In addition to the usual inlists (inlist and inlist_astero), you'll find 2 more:
inlist_model_search_controls and inlist_oscillation_controls.  As you might
guess, the 1st sets the search parameters and the 2nd has information
about the data to be matched.  The adipls package is used for doing the
frequencies, so there are control files for it as well.    There are also minor
changes to make/makefile and src/run.f.  The change to run is to call
the main driver routine in mesa/astero/run_star_adipls.
	
My hope is that this will be a useful starter-kit for users who'd like to try
using mesa to search for models to match the great data coming from
KEPLER and CoRot.  Let me know if you'd like to give it a try.

BTW: the example is for HD49385 -- take a look at this lovely paper:




Hooke-Jeeves minimization option for astero model searches


Searching for a best-fit model for a set of asteroseismic data is a minimization problem
for a function f of multiple parameters (e.g., chi^2 using mass, mixing length, Y, and Z/X).
Modern algorithms for minimization can get improved performance by using the values
of previous evaluations of f to estimate the shape of the search domain to decide where
to look next.  The BOBYQA algorithm is a state-of-the-art example of such routines.
This works great for "smooth" search domains but can get lost with domains that
include "rough" spots  such as invalid parameter combinations inside the allowed
search hypercube.  In such cases, it can be better to give up the possible improved
performance to get a more robust routine.  The old (1961) Hooke-Jeeves algorithm
is just such an option.  It only uses the function values to make better-or-worse
decisions; it doesn't use the values to estimate slopes in the search domain.
So occasional "invalid" parameter combinations don't cause it problems -- 
they can simply be given a very large artificial value causing them to be rejected.

For applications with a smooth search space (such as solar model calibration?),
BOBYQA is still one of the best options available.  But for a search space with
potential "drop-outs" such as a preliminary parameter scan for asteroseismology,
Hooke-Jeeves seems to be a better choice.   Both are now available in the astero
package.


We are using the f90 implementation from John Burkardt's website (thanks!!!).
http://people.sc.fsu.edu/~jburkardt/f_src/toms178/toms178.html

      !    Robert Hooke, Terry Jeeves,
      !    Direct Search Solution of Numerical and Statistical Problems,
      !    Journal of the ACM,
      !    Volume 8, Number 2, April 1961, pages 212-229.