[mesa-users] Fwd: Solar calibration

[Jean-Claude P. <jc...@gm...>]

Hi Hocine,

we had a discussion on the mailing list regarding this problem last January. Bill sent a very nice email summarizing his findings (see below). It contains all the inlists/models, and should be a good starting point.

Cheers,

JC

On May 5, 2014, at 5:44 PM, h....@cr... wrote:

> Hi,
> I can know how to work the solar calibration, by minimizing the chi2
> spectroscopic.
> In particular I tried to calibrate the sun for different abundances tables
> (GS98, ASG04, L03, AGSS09, L09) and I compared the sound speed of each
> model with the seismic data and I get the same results for all the models,
> which should not be the case. I try to understand why, can you help me.
> best regards
> hocine



Begin forwarded message:

> From: Bill Paxton <pa...@ki...>
> Subject: Re: [mesa-users] solar calibration with Asplund abundances
> Date: January 18, 2014 6:17:35 AM GMT+11:00
> To: Jean-Claude Passy <jc...@gm...>, Aldo Serenelli <al...@ic...>
> Cc: mesa-users group group <mes...@li...>
> 
> [This email is so long it needs an abstract!!!]
> 
> Abstract: compare doing calibration by newton root find to the alternative method used in mesa of doing calibration by chi^2 minimization.  give mesa results using opacities for gs98 abundances showing a good chi^2 match.  with a09 opacities we cannot find a comparably good chi^2.  finally redo the gs98 case adding oscillation frequencies to the chi^2 -- tiny sigmas for frequencies, a few parts in 10^5, result in a large chi^2 in spite of matches to frequencies better than a part in 10^3.
> 
> 
> 
> 
> On Jan 14, 2014, at 8:37 AM, Jean-Claude Passy wrote:
>> 
>> I was wondering if anyone had tried a solar calibration but with the Asplund abundances. I gave it a try and found the following accuracies:
>> 
>> 1% in Teff,
>> 5% in L,
>> 2% in Rcz.
> 
> 
> On Jan 14, 2014, at 11:50 AM, Aldo Serenelli wrote:
> 
>> One point is: your accuracy is certainly too low. The three quantities solar models are forced to satisfy are:
>> 
>> L_sun, R_sun, Z/X
>> 
>> and this can be done to parts per 10^5 easily by all codes around. So, your calibration right now is not good. 
> 
> 
> 
> 
> Hi JC and Aldo,
> 
> Mesa doesn't currently include a solar calibration package like the ones in other codes such as Aldo mentioned where you can expect 1 part in 10^5 accuracy for the matched values.  My understanding is that those codes do solar calibration as a newton root find with 3 unknowns (such as initial X, initial Z, and MLT alpha) and 3 equations relating model results to observations (e.g., model L = observed L_sun, model R = observed R_sun,  model surface Z/X = observed solar surface Z/X).  If you have reasonable input physics, a good starting guess, and do enough iterations, then as Aldo notes you should be able to get a part in 10^5 residuals for the 3 equations using the newton method.  After getting convergence for the newton solve, you can take that calibrated model and use other observables to evaluate the values the model gives for things like surface Y, radius at base of convection zone, or sound speed profile.   The key point is that the calibration is based on 3 observables, not for the full available set.   (btw: I had something like this in mesa in the old days -- it is in version 2891, e.g. -- back before we started to work on supporting asteroseismology.  I never liked it because of the need for numerical partials and good initial guesses made it fragile, and I happily threw it out once the astero package came into being.)
> 
> An alternative approach is the one we have in mesa for doing asteroseismology:  collect as many observables as you want into a chi^2 function and adjust a set of as many parameters as you want to find a minimal chi^2.  Unlike the previous newton root solve, this doesn't require same number of parameters as equations, so you can easily use more parameters if you wish (e.g., you might want to add a parameter for exponential overshooting), and you can use more observation data in the chi^2 -- for example, it can include surface Y, radius of base of convection zone, or sound speed profile so your optimization search will take more things into account when looking for a best fit.   
> 
> Note that rather than aiming for part in 10^5 accuracy for 3 numbers, we are hoping to get a chi^2 approaching 1 for a larger set of data.  For most of the data that we are trying to match 1 sigma is at least a part in 10^3 so we aren't trying to get 1 in 10^5 accuracy.   In normal cases the accuracy for calibrated L, calibrated R, and calibrated Z/X will not be as good with this method as it will be with the newton.  The optimization process will be happy to sacrifice 1 in 10^5 accuracy if sigma is 1 in 10^3 and by reducing the match for one observable, the model can get a smaller overall chi^2.
> 
> Both of these approaches have uses.  For work on the sun, the 1st method could well be superior -- Aldo might want to enlighten us on that issue.  I've provided the 2nd in mesa since our efforts have been focused more on asteroseismology than the calibration of solar models.  Our so-called solar calibration is done using our asteroseismology chi^2 minimization tools, the Nelder-Mead Simplex algorithm in particular.   Following are some examples to show how it does.
> 
> 
> [Note: keep in mind that I may have used incorrect values for sigmas for some of these (e.g., what's sigma for logL?).   Before taking any of this too seriously, you should double check the inlists and rerun these for yourself!    If you notice something bogus, please let us know.   The following is meant as a demonstration of capabilities, not science!]
> 
> 
> 1. solar calibration with GS98 abundances
> 
> Here's an example using OP opacities based on GS98 abundances.   As expected we are not getting 1 in 10^5 for Teff, L, and surface Z/X, but we are getting about 1 sigma matches to almost all of the 6 terms I included in the chi^2; also note that I've included exponential overshooting as a 4th parameter in addition to initial Fe/H, initial Y, and MLT alpha.   With a bit more work, we could extend this to take into account the small uncertainties in solar age and mass, i.e., mass and age could be both parameters and also terms in the chi^2.  And of course we could go nuts and try to match the oscillation frequencies too by adding them to the chi^2 (more about that below).  Note the good results for surface He, radius of base of convective zone, and sound speed profile rms error.
> 
> solar min chi^2 best model using GS98 abundances for opacities
> 
> 4 params (with values for best model)
> 
> init Fe/H = 9.0020828465566499D-02
> init Y = 2.7897021506175745D-01
> mlt alpha = 1.8999105557276836D+00
> exponential overshooting f = 3.5708368876502033D-03
> 
> 6 terms in chi^2 (with target, sigma, and best model result for each)
> 
> best model chi^2 = 0.8 (not reduced for number of parameters)
> 
> Teff_target = 5777d0
> Teff_sigma = 65d0
> Teff = 5.7919731102619780D+03
> 
> logL_target = 0.00d0
> logL_sigma = 0.05d0
> logL = 4.7418072907512709D-03
> 
> surface_Z_div_X_target = 2.292d-2 ! GS98 value
> surface_Z_div_X_sigma = 1d-3
> surface_Z_div_X = 2.4662720344812331D-02
> 
> surface_He_target = 0.2485d0 ! Basu & Antia 2004
> surface_He_sigma = 0.0035
> surface_He = 2.5133918425690588D-01
> 
> Rcz_target = 0.713d0 ! (radius of base of convective zone) Basu & Antia 1997
> Rcz_sigma = 1d-3
> Rcz = 7.1356906450350521D-01
> 
> sound speed profile (rms error for r=0.094 to 0.94 Rsun)
> target = 0
> sigma 1d-3
> csound_rms = 8.6368382612381526D-04
> 
> other info for best model
> 
> initial h1 = 7.0124662477067556D-01
> initial he4 = 2.7894231804025127D-01
> initial Z = 1.9783160167566916D-02
> logR = -2.3333789864714455D-05
> delta_nu = 142  (target = 135, sigma = 1)
> 
> 5841
> 
> 
> 
> 
> 
> 
> 2. solar calibration with A09 abundances
> 
> Rerunning with opacities based on Asplund abundances (a09) fails to find as good a match.   I used the OP data provided by Radek Smolek and removed the small number of NaN's by hand.  Those were used in conjunction with the low T opacities from Ferguson for the a09 abundances.  the same parameters and targets were used, with surface_Z_div_X_target changed to 1.81d-2 and the lower bound for the FeH parameter dropped to -0.1.
> 
> solar min chi^2 best model using A09 abundances for opacities
> 
> 4 params (with values for best model)
> 
> init Fe/H = -3.0272265702073010D-02
> init Y = 2.7355443746545149D-01
> mlt alpha = 1.8934468858740710D+00
> exponential overshooting f = 3.8623001090813583D-02
> 
> 6 terms in chi^2 (with target, sigma, and best model result for each)
> 
> best model chi^2 = 20
> 
> Teff_target = 5777d0
> Teff_sigma = 65d0
> Teff = 5.8522865437621213D+03
> 
> logL_target = 0.00d0
> logL_sigma = 0.05d0
> logL = 2.1504566243727114D-02
> 
> surface_Z_div_X_target = 2.292d-2 ! GS98 value
> surface_Z_div_X_sigma = 1d-3
> surface_Z_div_X = 1.9108543180648075D-02
> 
> surface_He_target = 0.2485d0 ! Basu & Antia 2004
> surface_He_sigma = 0.0035
> surface_He = 2.4893202779985088D-01
> 
> Rcz_target = 0.713d0 ! (radius of base of convective zone) Basu & Antia 1997
> Rcz_sigma = 1d-3
> Rcz = 7.2282050191100233D-01
> 
> sound speed profile (rms error)
> target = 0
> sigma 1d-3
> csound_rms = 4.5104639148126815D-03
> 
> other info for best model
> 
> initial h1 = 7.1123501257440602D-01
> initial he3 = 2.7355443746545149D-05
> initial he4 = 2.7352708202170495D-01
> initial Z = 1.5210549960142587D-02
> logR = -6.4004483584790260D-04
> delta_nu = 1.4238961741912979D+02 (target = 135, sigma = 1)
> 
> 
> 
> 
> 
> 
> 
> 3. GS98 results when oscillation frequencies are included in the chi^2
> 
> now for giggles, let's see what happens if we rerun gs98 adding oscillation frequencies to the chi^2.vwe'll keep the same 4 parameters; the same 6 chi^2 terms from above will now be called the "spectral" chi^2;and we'll add 74 frequencies as a "seismic" chi^2. the total chi^2 will be the average of the spectral and seismic. as you might expect, the optimization process gives up quality in the spectral chi^2 in order to improve the seismic result. but in spite of that it isn't able to find a good overall match. in fact, even though the spectral chi^2 has jumped from 0.8 above to 8.8 now, we still have a chi^2 seismic of 1863, so our total chi^2 is a whopping 936 -- ugh.   one small bright spot is that the delta_nu value now is within 1 sigma of the observed. but the frequencies are so well know that they have very small sigmas and that gives us the big chi^2 seismic.
> 
> for example, here's the results for the lowest order radial frequency.
> obs = 1548.328 microHz
> model = 1547.604
> sigma = 0.037
> we are off by 20 sigma, so even though within 1 microHz (4 in 10^4), it gives a big contribution to the chi^2.
> 
> it might be possible to get a better model by giving the spectral chi^2 a higher weight than the seimsic one so that the tiny sigmas for frequencies don't make the seismic values dominate the total chi^2 as they do now.
> 
> Here's a summary of what we get when we try to match frequencies as well as the usual things.
> 
> solar min chi^2 best model using GS98 abundances for opacities, including frequencies in chi^2.
> 
> 4 params (with values for best model)
> 
> init Fe/H = 9.6964541300850082D-02
> init Y = 2.7055352082138318D-01
> mlt alpha = 1.5559254176170356D+00
> exponential overshooting f = 4.3651916415638134D-03
> 
> 6 terms in "spectral" chi^2 (with target, sigma, and best model result for each)
> 19 l=0, 19 l=1, 20 l=2, and 16 l=3 frequencies (74 total) in "seismic" chi^2
> total chi^2 average of the spectral and seismic
> 
> best model total chi^2 = 936
> chi^2 spectral = 8.8
> chi^2 seismic = 1863
> 
> Teff_target = 5777d0
> Teff_sigma = 65d0
> Teff = 5.4684562018480437D+03
> 
> logL_target = 0.00d0
> logL_sigma = 0.05d0
> logL = -9.5437344710497446D-02
> 
> surface_Z_div_X_target = 2.292d-2 ! GS98 value
> surface_Z_div_X_sigma = 1d-3
> surface_Z_div_X = 2.5276077534231767D-02
> 
> surface_He_target = 0.2485d0 ! Basu & Antia 2004
> surface_He_sigma = 0.0035
> surface_He = 2.4543028071293943D-01
> 
> Rcz_target = 0.713d0 ! (radius of base of convective zone) Basu & Antia 1997
> Rcz_sigma = 1d-3
> Rcz = 7.1533417507751329D-01
> 
> sound speed profile (rms error)
> target = 0
> sigma 1d-3
> csound_rms = 3.8666863830618473D-03
> 
> other info for best model
> 
> initial h1 = 7.0911880920827919D-01
> initial he3 = 2.7055352082138320D-05
> initial he4 = 2.7052646546930104D-01
> initial Z = 2.0327669970337723D-02
> logR = -1.8931209827114938D-04
> delta_nu = 1.3589598957704791D+02  (target = 135, sigma = 1)
> 
> 
> 
> 
> 
> 
> 
> One final reminder: the simplex optimizer does a bounded search, so when you set things up in the inlist_astero_controls, you need to set reasonable bounds as well as starting values.  If you set bounds that are too tight, you are sure to fail.   For example, if you switch from gs98 to a09 but keep the search bounds that assume gs98 you'll be unhappy with the result.
> 
> And as I mentioned above, if you would rather play the newton root find game instead doing minimization with simplex, feel free to go for it!
> 
> Cheers,
> Bill
> 
> 
> 
> 
> 
> On Jan 14, 2014, at 11:50 AM, Aldo Serenelli wrote:
> 
>> Jean-Claude, 
>> 
>> regardless of the wrongness of MLT and all other possible sources of errors, your solar calibration should reproduce other solar calibrations, provided the input physics is the same (or very similar). 
>> Your quoted helium abundance is certainly too high compared to all other solar calibrations in the market that incorporate gravitational settling in the model.
>> 
>> One point is: your accuracy is certainly too low. The three quantities solar models are forced to satisfy are:
>> 
>> L_sun, R_sun, Z/X
>> 
>> and this can be done to parts per 10^5 easily by all codes around. So, your calibration right now is not good. Since L_sun and R_sun are satisfied by construction, so it will be Teff. Keep in mind the Sun evolves very slowly, so a 5% difference in L is actually a very long period of time.
>> 
>> You should therefore check what is causing the difference. Here are some possibilities:
>> 
>> 1) improve accuracy
>> 2) is grav. settling included?
>> 3) opacities consistent with assumed (AGSS09) abundance?
>> 
>> You won't get Rcz right, of course, but your results should be overall much more consistent with published literature.
>> 
>> Aldo
>> 
>> On 14 Jan, 2014, at 5:37 PM, Jean-Claude Passy wrote:
>> 
>>> Hi all,
>>> 
>>> I was wondering if anyone had tried a solar calibration but with the Asplund abundances. I gave it a try and found the following accuracies:
>>> 
>>> 1% in Teff,
>>> 5% in L,
>>> 2% in Rcz.
>>> 
>>> Comparing with Serenelli et al. (2009):
>>> http://esoads.eso.org/abs/2009ApJ...705L.123S
>>> 
>>> these results are consistent in terms of Rcz, and better regarding the He surface abundances (here Ysurf = 0.2477). I used MESA version 5595 (inlists and output attached).
>>> 
>>> Has anybody tried to do the same thing and found better accuracies?
>>> 
>>> Cheers,
>>> 
>>> JC
>>> 
>>> <inlist_astero_search_controls3>
>>> <inlist_solar3>
>>> <sample_0001.data3>