CPM: Issues with mass_conv_core reporting

[rjfarmer]

From Anne:

It seems like that with the new CPM there are issues with reporting history items like mass_conv_core and conv_mx1_top are not reported properly and do not map to the actual convective/semiconvective boundary.

Investigating

[evbauer]

I don't think this issue should be closed yet? The comment in controls.defaults is good enough for going ahead with a release, but I think there's still an underlying issue that we want to improve. Can I reopen?

[warrickball]

Can someone create a minimal working example of the problem so it's easier for other to get started on this? I failed to find anything in the mailing lists. I just ran this basic inlist through r12778 because it's convenient. I can switch to a specific dev version if someone tells me which.

&star_job

      create_pre_main_sequence_model = .true.
      pre_ms_T_c = 9d5

/ ! end of star_job namelist                                                                                                                          

&controls
      history_interval = 1
      use_dedt_form_of_energy_eqn = .true.
      initial_mass = 2.0

! mixing controls                                                                                                                                     

      use_Ledoux_criterion = .true.
      do_conv_premix = .true.

! When to stop                                                                                                                                        

      xa_central_lower_limit_species(1) = 'h1' ! isotope name as defined in chem_def                                                                  
      xa_central_lower_limit(1) = 1e-3

/ ! end of controls namelist                                                                                                                          

Here's a comparison of different convective boundary options. I presume the relative jaggedness of mass_conv_core is the problem.

[warrickball]

The default controls include these:

num_cells_for_smooth_gradL_composition_term = 3
threshold_for_smooth_gradL_composition_term = 0

num_cells_for_smooth_brunt_B = 2
threshold_for_smooth_brunt_B = 0d0

I think this allows MESA to smooth gradL and brunt_B across the semiconvective boundaries, which can leave a small gradL excess on the convective side of the core boundary. When MESA computes the mixing properties, it sees the convective boundary at the bottom of this smoothing region, which is (with defaults) ~2 mesh points below where gradr == grada.

There are two workarounds. The first, used in the MESA V inlists, is to increase the thresholds to some non-zero values that prevent smoothing across the boundaries. Specifically, the MESA V inlists used

num_cells_for_smooth_gradL_composition_term = 50
threshold_for_smooth_gradL_composition_term = 0.02

num_cells_for_smooth_brunt_B = 50
threshold_for_smooth_brunt_B = 1E-4

Another option is to disable the smoothing entirely (i.e. set both num_cells... controls to 0). Here's what that does to the example above: where mass_gradr_grada_core is interpolates in gradr-grada to find the core boundary using this code snippet in run_star_extras.f90:

         vals(1) = 0d0
         if (s% gradr(s% nz) < s% grada(s% nz)) then
            vals(1) = 0d0
         else
            vals(1) = s% xmstar/Msun
            do k = s% nz-1, 1, -1
               if (s% gradr(k) < s% grada(k)) then
                  x1 = s% q(k+1)
                  x2 = s% q(k)
                  y1 = s% gradr(k+1) - s% grada(k+1)
                  y2 = s% gradr(k) - s% grada(k)
                  a = (y2-y1)/(x2-x1)
                  vals(1) = (x1-y1/a)*s% xmstar/Msun
                  exit
               end if
            end do
         end if

mass_conv_core is still a bit jagged but that looks at least in part like a mesh resolution issue. Using the thresholds from MESA V has the same effect.

I also found the control recalc_mix_info_after_evolve, which I thought would help but now it looks like it doesn't make a difference.

[warrickball]

I was testing on r12778 but @annethoul showed me that the situation is actually better in r14427. In that case, recalc_mix_info_after_evolve does make a difference and produces a pretty clean convective core mass for a 5Msun model. The 2Msun example isn't as smooth, especially, during the convective core's early growth but this can be mitigated by increasing the mesh resolution. I'm waiting to hear from @annethoul about how these additional options (recalc... and the thresholds) perform in other worrisome areas.

The only things that changed in star/private/conv_premix.f90 between r12778 and r14427 are the addition of controls for timing information (r13297) and line 197 (added r12877):

     s% need_to_setvars = .true.

which appears before the return from subroutine do_conv_premix. So between that and the extra controls, we might have a solution to this issue.

[rjfarmer]

With the re-addition of the convective_bdy_weight controls is this fixed now or at least better?

[evbauer]

I think this is a somewhat different issue conceptually. I would want to check with @annethoul to be sure, but my recollection is that there was some inconsistency between what MESA reports as the core boundary and the location of material that gets mixed by semiconvection when CPM is on.

[annethoul]

Yes, Evan is correct. These issues are totally separate. convective_bdy_weight is something that actually interferes with the models, to evolve the convective boundary more smoothly from timestep to timestep. The problem described in this thread is that the models are good, and if you use profiles to figure out the convective (and "partially mixed" ) boundaries, you get the correct positions, while those boundaries are incorrectly reported in the history files. For MESA 5 figures, Rich figured out the boundaries using the profiles. But apparently it was not straightforward to do. Still waiting for a volunteer to write a piece of code that finds those boundaries from the profiles data (from I assume, kinks in the abundance profiles?). It is very easy to position those boundaries "by eye" looking at profiles, apparently much harder to code.

[evbauer]

@annethoul I think I recall you had some plots that pretty specifically isolate where the location should be based on the profiles vs where it is reported in the history files. Do you have those available to post here? Or better yet, could you generate a minimal example on the latest release?

[annethoul]

@evbauer I do not have those plots anymore. I will try to create such plots for you. It might take a few days...

[annethoul]

Here is one example (inlist based on the mass_conv_core test suite example, with better resolution). Looking at the history plot of mass_conv_core versus model number, one sees "peaks" or "dips" in the convective core mass. Those are due to tiny splittings of the core. CPM very quickly recovers from those small splittings, and by eye one can easily follow the evolution of the convective core mass, but in the history file, the value reported for the mass of the convective core is the first convective zone from the center. Note that those tiny splittings also mess up the profile of he4, and it is only when the core size increases and "swallows" the messed up region that we recover a decent profile for he4. The exact evolution of he4 is probably not totally correct due to those small splittings.

Let me know if this is unclear!

[annethoul]

history

[annethoul]

profile 1256

[annethoul]

profile 1257

[annethoul]

profile 1258

[annethoul]

profile 1350

[annethoul]

@wmwolf Note: all figures generated using Bill's all new and wonderful MESA Explorer tool!

[annethoul]

Also, kind of separate issue but linked to this one, once a "semiconvective" region (due to CPM) develops, MESA has no way to locate its boundary. Again, easy to locate "by eye" looking at e.g. he4 profiles, but hard to code.

[fxt44]

mass_conv_core test suite example

i don't see this test suite case in 15140 or 23.05.1. are you pulling it from another test suite case? or maybe i'm just being silly.

Those are due to tiny splittings of the core.

one of your hand-drawn figs explains it pretty well:

Note that those tiny splittings also mess up the profile of he4, and it is only when the core size increases and "swallows" the messed up region that we recover a decent profile for he4.

with crazy small timesteps and high mass resolution, it is possible to avoid the splitting until x(he4) drops below ~1e-3, which is still an issue if one cares about the central o16/c12. its more fun further out in the star :)

[annethoul]

@fxt44 my fault, it is called conv_core_cpm ... told you I am getting old ;-( I am amazed you kept those hand written figures... I do miss the KITP blackboards...

[evbauer]

Also, kind of separate issue but linked to this one, once a "semiconvective" region (due to CPM) develops, MESA has no way to locate its boundary. Again, easy to locate "by eye" looking at e.g. he4 profiles, but hard to code.

Oh interesting. This is what I thought the main point of this issue was. This part is strictly a reporting issue, in the sense that we know where the boundary ought to be intuitively, but still have some work to do to get MESA to report the right number.

The other issue is more of an actual physics issue that's having an impact on the evolution of the models, and probably should motivate some sort of improved resolution controls to prevent the splitting from happening.