[mesa-users] Understanding opacity blending: is smoothing function being applied?

[Warrick Ball]

Bill,

Great, thanks for the thorough reply.  Your reminder about results 
depending on engineering and design choices probably applies equally well 
to any evolution code.  But it's great that in MESA, not only are they 
available to be seen, but the code and controls are sufficiently 
transparent that an outsider can make sense of them!

I'll try hardcoding it back to the smoothed blend, just to see if it helps 
with my particular scenario.  If it turns out to be helpful, maybe it'll 
be worth adding a control, but I don't think it's necessary unless I 
discover that it helps.  (At which point I'll probably also try 
reintroducing the partial derivatives, for self-consistency.)

Cheers,
Warrick

PS:  While MESA's engineering decisions might have no basis in physics, 
they usually have a sound basis in engineering! :p


On Tue, 29 Oct 2013, Bill Paxton wrote:

> 
> On Oct 29, 2013, at 9:42 AM, Warrick Ball wrote:
>
>       I've been digging around in the opacity routines (mesa/star/kap/private) and, unless I am reading the
>       code incorrectly, it looks to me like opacities are not blended according to the formula presented in
>       Paper I (equation 1). ...... [see below for rest of email]
> 
> 
> 
> Hi Warrick,
> Great to have more people looking at the code!  I appreciated it very much, and I expect you will find bugs.  ;p  
> 
> The opacity blending is a case of an engineering decision that has no basis in physics --- it is forced on us by
> the need to do combine sources that cover partially overlapping regions of (T,Rho).
> 
> 
> 
> AN IMPORTANT REMINDER FOR ALL MESA USERS: results are only partly determined by the choice of physics parameters
> -- they are also sensitive to engineering decisions that have no basis in physics.  I think one of the values of
> mesa is that these engineering choices are out in the open for discussion and testing.  
> 
> 
> 
> In mesa/kap we have a patchwork of opacities from different sources that must be stitched together to give a more
> complete coverage of (T,Rho) space.   We need to spread out the transitions between sources over a range of values
> in order to avoid jumps in opacity that give bogus partials and create impossible convergence problems for
> mesa/star.  So, for example, a blend in logT will extend over a range logT_lo to logT_hi with source1 used for
> logT <= logT_lo, source2 used for logT >= logT_hi, and a combination of the sources used for intermediate values.
>  
> 
> The simplest way to combine the sources is just to do a transition that is linear in logT:
> kap1 = opacity according to source1 (for lower T)
> kap2 = opacity according to source2 (for higher T)
> 
> if (logT <= logT_lo)  then just use kap1
> if (logT >= logT_hi)  then just use kap2
> else blend as follows:
> 
> alfa = (logT - logT_lo) / (logT_hi - logT_lo) ! alfa is fraction of the high T source (kap2)
> kap = kap1*(1 - alfa) + kap2*alfa
> 
> In the old days of Paper 1, we used a variation on this:
> 
> f = (logT - logT_lo) / (logT_hi - logT_lo) ! f is fractional distance
> alfa = 0.5*(1.0 - cos(f*pi)) ! alfa is fraction of the high T source (kap2)
> kap = kap1*(1.0 - alfa) + kap2*alfa
> 
> The linear blend has the benefit of spreading the transition from source1 to source2 evenly over the entire blend
> range; it has the disadvantage of discontinuous slopes at the edges of the blend range.
> 
> The cosine form removes the boundary slope discontinuities but as a result it focuses the major part of the
> transition into a smaller central region of the blend range.
> 
> Which one is better?  Who knows???   At the time of Paper 1 I thought the cosine form was the way to go.  Since
> then I must have run into cases that worked better with the linear form leading me to make a switch to the simpler
> linear form.  Or maybe I just decided that the cosine form didn't actually make any improvement, so I might as
> well just use the simpler linear form.  Perhaps I should have introduced yet another control parameter for this
> choice!  I can certainly do that now if you'd like to experiment --- just let me know.
> 
> Here's one other disgusting little detail you might not have noticed.  We are returning partials of the opacity
> with respect to lnT and the blending factor depends on lnT, so we should include the terms for the partial of the
> blending factor in the partials of the opacity, right?  But we don't include it -- I tried it long ago, and it
> seemed to make things worse.  Perhaps I just screwed up, but I removed that extra term and the convergence in star
> improved.  I found the same to hold in the temperature dependent blends in the EOS.  
> 
> Cheers,
> Bill
> 
> 
> 
> 
>
>       Bill,
>
>       Me again...  I've been digging around in the opacity routines (mesa/star/kap/private) and, unless I am
>       reading the code incorrectly, it looks to me like opacities are not blended according to the formula
>       presented in Paper I (equation 1).
>
>       Let me start be restating that formula.  On page 7 of Paper I, it says that, in the blending region,
>       the opacities are blended by first defining
>
>       F = (log T - log T_L) / (log T_U - log T_L) ... (1)
>
>       and
>
>       S = (1 - cos(F*pi))/2 ... (2)
>
>       where T is the temperature, T_L is the cooler boundary of the blending region and T_U is the hotter
>       boundary.  Then the code returns
>
>       log(kappa) = S*log(kappa_U(R,T)) + (1-S)*log(kappa_L(R,T)), ... (3)
>
>       where kappa_U is the opacity source above the blend (e.g. OPAL), kappa_L the source below (e.g.
>       Ferguson) and R is the local value of the density parameter.
>
>       I had a look for this inside the code.  First, I did some grepping in mesa/star/kap:
>
>             grep cos private/*.f
>
>       private/kap_eval_co.f:95:            !alfa = 0.5 * (1.0 - cos(alfa * pi))
>       private/kap_eval_co.f:1255:          fac = 0.5*(1 - cos(3.14159*fac))
>       private/kap_eval.f:138:              !fac = 0.5 * (1.0 - cos(fac * pi))
>       private/kap_eval_fixed.f:94:         !alfa = 0.5 * (1.0 - cos(alfa * pi))
>       private/kap_eval_fixed.f:805:        fac = 0.5*(1 - cos(3.14159*fac))
>
>       (I've done some re-alignment for easy reading.)  That's the formula alright, but why is it almost
>       always commented?  I looked, for example, in kap_eval_fixed.f and found that the commented line
>       belongs to Get_kap_fixed_metal_Results, which is called by Get_kap_Results in kap_eval.f, which is in
>       turn called by kap_get_blend_1_2 in public/kap_lib.f.  And that's ultimately called (through
>       kap_get_Type1 or 2) in mesa/star.  What I couldn't find is any other occurence of the smoothing
>       formula (2).  I may have missed a comment in Paper II that this had changed.  So I don't think the
>       correction is being applied at the low temperature end.
>
>       For the record, there are the two *uncommented* occurences but both appear in the blend with Compton
>       scattering at the high-temperature end.
>
>       My question, I guess, is "why is this the case?"  I suppose it probably doesn't matter that much
>       (unless you're doing something crazy like calculating an outrageously fine sequence of models... :p),
>       but then why was the formula devised and implemented in the first place?  I haven't tried uncommenting
>       those and recompiling because I have a whole lot of runs in our cluster queue that would be disrupted
>       but I'll update this when I get a chance.
>
>       And I don't want to send you a "fix" when you had a good reason to change it in the first place!
>
>       Help and wisdom, as ever, greatly appreciated.  By the by, this is MESA 5527.
>
>       Cheers,
>       Warrick
> 
>
>       ------------
>       Warrick Ball
>       Postdoc, Institut für Astrophysik Göttingen
>       wball at astro.physik.uni-goettingen.de
>       +49 (0) 551 39 5069
> 
> 
> 
>


------------
Warrick Ball
Postdoc, Institut für Astrophysik Göttingen
wball at astro.physik.uni-goettingen.de
+49 (0) 551 39 5069