[mesa-users] optical_depth -- and a message for MESA BLACK BELT USERS

[jingluan at caltech.edu]

>> Hi,
>>
>> Let's see if there are any black belts out there who can help on this.
>>
>> But first you need to provide enough info to make it possible for them.
>> Make a saved model a few steps before the bad one.
> Thank you, Bill :-)
>
> I start with WD0.6-Hrich-purecooling.mod, which is attached. I evolved a 7
> Msun star from pre-main sequence. When the helium core is 0.5Msun, I
> started losing mass on a thermal timescale until the total mass is
> 0.6Msun. Then I terminated the mass-loss and evolves the 0.6Msun star
> until no nuclear burning is going on, i.e., it is purely cooling down
> white dwarf which is saved as this model I attached here.
>
>> Make a modified inlist that starts by loading that newly saved model.
>
> inlist_project is also attached. I change the net to 17_to_fe56.net which
> I also used to make the starting model from pre-main sequence star.
>
>> Confirm that the problem still happens when you run with the new inlist.
>> Send the starting model and the modified inlist.
>> If possible, provide a plot to show where things are bad.
>> explain what is wrong in enough detail so someone else can see the
>> problem.
>
> What I check is quite simple. I first find the shell number kphot for the
> photosphere. I expect s% tau(kphot)~0.66. But s% tau(kphot)~360.
>
> then I thought maybe s% tau(:) is not the optical depth. Then I add some
> lines in src/run_star_extras.f(also attached) to calculate the optical
> depth (tau_phot):
>
>          kphot=0
>          tau_phot=0.
>          do k=1, s% nz
>             if (s% R(k)/rsol<= s% photosphere_r) then
>                kphot=k
>                tau_phot=tau_phot+s% opacity(k)*s% rho(k)*(s% R(k)-s%
> R(k+1))
>                exit
>             end if
>          end do
>
> And it turns out that tau_phot is an even crazily larger number than s%
> tau(kphot).
>
> Then I thought maybe mesa had not yet converged but it did not report any
> err or warning information which made me worried that "Do I need to check
> manually every variable in s even if mesa seems to run happily?"
>
> Thank you very much!
>
> Sincerely,
> Jing
>
>> Include the log file for that model.
>> include version number info for mesa, the compiler, and the operating
>> system.

Sorry that I forgot this info:

mesa 2664; gfortran 4.6.0; Ubuntu 10.04.4, 64-bit version.

Sincerely,
Jing




>>
>> -Bill
>>
>>
>>
>>
>> On Mar 23, 2012, at 10:18 AM, jingluan at caltech.edu wrote:
>>
>>> Dear Bill,
>>>
>>> Thank you so much for guiding me. Although I have not yet figured out
>>> what
>>> specific problem causes the strange tau(:), I will try to follow yours
>>> to
>>> see.
>>>
>>> I still have one question. If this weird tau were caused by
>>> non-convergence, I would NOT be able to find it out if I did not
>>> accidentally check the tau values. MESA does not report err information
>>> like 'convergence problem', normally I would accept this model as a
>>> converged one and trust the outputs by mesa, e.g. the figures plotted
>>> by
>>> pgplot. But actually there is something wrong in the model (e.g. the
>>> tau
>>> in this case) which is NOT revealed. I feel that in practice people do
>>> not
>>> check every variable's value in s to see whether it makes sense if MESA
>>> does not report error.
>>>
>>> Thank you so much!
>>>
>>> Sincerely,
>>> Jing
>>>
>>>
>>>>
>>>>
>>>>
>>>>
>>>>>>>  MESA USERS
>>>>
>>>> Please consider that if you are publishing papers using mesa,
>>>> it might be a good idea to be able to poke around inside the code
>>>> enough
>>>> to get an idea of what is going on.  If you can do that, you are a
>>>> mesa
>>>> black belt!
>>>> Are you a black belt yet?  Do you want to become one?
>>>> If so, take a moment to see if you could have answered this email
>>>> yourself.
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> On Mar 23, 2012, at 9:00 AM, jingluan at caltech.edu wrote:
>>>>>
>>>>> At the kth shell where it satisfies s% R(k)/rsol==s% photosphere_r, I
>>>>> found that s% tau(k)=361 in a model by mesa2664;
>>>>>
>>>>> Then I calculate a variable called tau_phot by summing s%
>>>>> opacity(j)*s%
>>>>> rho(j)*(s% R(j)-s% R(j+1)) over j=1,..,k; and I got tau_phot=5184;
>>>>>
>>>>> It is pretty weird. What could be the problem? The model has not yet
>>>>> been
>>>>> converged well? But mesa does not report any ierr information though.
>>>>
>>>>
>>>> Hi Jing,
>>>>
>>>> Good question!   Rather than just giving an answer, let's go through
>>>> the process of discovering how things work (which is what I just
>>>> did!).
>>>>
>>>> Question: why do I find weird value for tau before the model has
>>>> converged.
>>>>
>>>>
>>>> 1st step: grep to find where tau is set.
>>>> 	cd star/private
>>>> 	grep tau *.f -n > lst
>>>>
>>>> 	open the file lst and don't panic when you see a lot of lines
>>>> 	we don't care about places that are reading tau;
>>>> 	we want to find where it is being set.
>>>> 	this looks like the best bet:
>>>>
>>>> star_utils.f:566:            tau(k) = tau(k-1) + dtau
>>>>
>>>> 	that is in a routine called "get_tau" which is a sign that we're on
>>>> the
>>>> right track.
>>>> 	open star_utils and take a look at get_tau -- bingo!
>>>>
>>>>
>>>> 2nd step: find where and when get_tau is being called.
>>>>
>>>> grep 'call get_tau' * -n
>>>> report.f:203:         call get_tau(s, s% tau)
>>>>
>>>> 	open report and find the call on get_tau.  it is in a routine called
>>>> do_report
>>>>
>>>>
>>>> 3rd step: find where and when do_report is being called.
>>>>
>>>> grep 'call do_report' * -n
>>>> evolve.f:285:         call do_report(s, ierr)
>>>> read_model.f:67:         call do_report(s, ierr)
>>>>
>>>> 	we find 2 places: one in evolve and the other in read_model.
>>>> 	since we're doing evolution rather than reading, we want the one in
>>>> evolve.
>>>>
>>>> 	open evolve.f and find the call on do_report.
>>>> 	turns out that it is called from do_evolve_step,
>>>> 		after the model has converged (in do_struct_burn_mix) .
>>>>
>>>>
>>>> understanding the order of events in do_evolve_step is the hard part
>>>> of
>>>> this.
>>>> there's a lot of extra things in there for error checks and
>>>> timekeeping,
>>>> but
>>>> fundamentally here's the order in which things are done:
>>>>
>>>> 	prepare_for_new_step -- does remesh
>>>> 	prepare_for_new_try -- save stuff in case need to do retry or backup
>>>> later
>>>> 	do_winds -- calculate mdot for this step
>>>> 	do_adjust_mass -- add or remove material according to mdot
>>>> 	do_set_vars -- set variables for new mesh and mass
>>>> 	do_element_diffusion -- optional
>>>> 	save_pre_hydro_values -- save starting values for solver
>>>> 	do_struct_burn_mix -- solve for new model
>>>> 	smooth_newly_non_conv -- optionally smooth abundances
>>>> 	do_brunt_N2 -- optionally calculate brunt N^2 at each cell
>>>> 	do_report -- calculate information for logs and profiles
>>>>
>>>> One more subtle point that "black belt" mesa user's should know:
>>>> the diffusion coefficients for convective mixing, overshooting,
>>>> semiconvection,
>>>> thermohaline mixing are calculated in the call on do_set_vars from
>>>> do_evolve_step and remain fixed during the rest of the calls from
>>>> do_evolve_step;
>>>> in particular, they are not updated between iterations of the solver,
>>>> so
>>>> they
>>>> are not being treated in a fully implicit manner.  this is done for
>>>> numerical reasons only.
>>>>
>>>> Okay --- now I expect that there will be lots of mesa users out there
>>>> in the future who'll be able to dig out the answers to questions like
>>>> this.
>>>> Consider it part of your qualification for getting a mesa black belt.
>>>> ;-)
>>>>
>>>> Cheers,
>>>> Bill
>>>>
>>>>
>>>>
>>>>
>>>>
>>>
>>>
>>
>>
>