one-zone burner at constant density

[sunnywong314]

I was trying to use the one-zone burner at constant density from $MESA_DIR/net/test , with an initial composition of pure helium, logT = 9 and logRho = 5, and I noticed that the temperature didn't evolve much (changed only by 1d-8) when I used a net like mesa_125.net . The simplest net that shows this problem is one with only he4 and c12. However, the temperature does evolve (logT changes by 0.2ish) if I use approx21.net, or if I use other nets with use_3a_fl87 = .true. (e.g., by adding the line g% use_3a_fl87 = .true. near line 236 in $MESA_DIR/net/test/mod_one_zone_burn.f90).

I dug deeper with the he4+c12-only net. I added some print statements and found that eps_nuc returned by eval_net (called by burner_derivs in $MESA_DIR/net/private/net_burn_const_density.f90) is zero , while dxdt(:) gives the right rates. I should also note that eval_net is also called by burner_jakob but returns the right eps_nuc.

I used the same conditions (pure He, logT=9, logRho=5 and he4+c12-only net) with sample_net (which also lives in $MESA_DIR/net/test) and got the right eps_nuc, so I think this behavior seems to only affect the one-zone burner. I'm using 24.08.1, so reverse 3a and detailed balance with >2 reactants are also not the problem.

I'm attaching my inlist and net here.
he4_and_c12.net.txt
inlist_one_zone_burn_pureHe.txt

Any pointers/suggestions are greatly appreciated!

Version

• MESA r24.08.1
• MESA SDK VERSION: aarch64-macos-23.10.1
• OS: MacOS Ventura 13.1

[pmocz]

I need @fxt44 or @Debraheem to help me understand the code logic in net/private/net_derivs.f90.

I see that the function get_derivs() calls other subroutines like update_special_rates(), get1_derivs(), get_general_1_to_1_derivs(), get_general_2_to_1_derivs(), get_general_2_to_2_derivs(), get_basic_2_to_2_derivs(), get_basic_2_to_1_derivs(). There is a flag just_dydt. The first two subroutines (update_special_rates(), get1_derivs()) always update the variable eps_nuc_MeV. But the other subroutines (e.g. get_general_1_to_1_derivs,...) exit out with if (just_dydt) return before eps_nuc_MeV is updated.

My questions are, is this the intended behavior, and why does this happen?

[Debraheem]

I'm not at my computer yet, but i think "just_dydt" is the flag for for split burn. If it's true you skip calculating the derivatives like "d(dyidt)dyj" and all the jacobian terms for composition derivatives are 0. The derivs functions are only called when fully coupled.

[Debraheem]

It looks like star/private/struct_burn_mix.f90 which is where the bruning is done in star. Here, when split_burn is turned on in do_burn -> burn1_zone-> call net_1_zone_burn_const_density. This is where split_burn is activated in the evolution. Then it calls to net_get (from net/public/net_lib) but with "just_dxdt = .false.", meaning we don't want any mixed derivative information. So the constant density solver is used for one zone burning and the eps_nuc_MeV term from net_derivs is not used in this circumstance to calculate eps_nuc. Either eps_nuc_MeV is set somewhere else, or there is a separate variable doing the same job.

[Debraheem]

I believe net_1_zone_burn_const_density is called directly out of star into net and hands the eps_nuc_MeV information back to star in the form of the vairable "avg_eps_nuc" wherein struct_burn_mix.f90, there is line 916:
s% burn_avg_epsnuc(k) = avg_epsnuc,

this leads back to star/private/hydro_energy.f90, where we have on line 234+:

     else if (s% op_split_burn .and. s% T_start(k) >= s% op_split_burn_min_T) then
               eps_nuc_ad = 0d0
               eps_nuc_ad%val = s% burn_avg_epsnuc(k)
            else
               eps_nuc_ad = 0d0
               eps_nuc_ad%val = s% eps_nuc(k)
               eps_nuc_ad%d1Array(i_lnd_00) = s% d_epsnuc_dlnd(k)
               eps_nuc_ad%d1Array(i_lnT_00) = s% d_epsnuc_dlnT(k)
            end if

[Debraheem]

I added some print statements and found that eps_nuc returned by eval_net (called by burner_derivs in $MESA_DIR/net/private/net_burn_const_density.f90) is zero , while dxdt(:) gives the right rates. I should also note that eval_net is also called by burner_jakob but returns the right eps_nuc

This seems expected, as eps_nuc from eval_net is not used when the constant_density solver is called.

[pmocz]

Thanks, @Debraheem for the explanations!

I see that if we move the guards if (just_dydt) return down a few lines in the functions {get_general_1_to_1_derivs(), get_general_2_to_1_derivs(), get_general_2_to_2_derivs(), get_basic_2_to_2_derivs(), get_basic_2_to_1_derivs()} in the file net/private/net_derivs.f90 so that all switch cases of get_derivs() update the variable eps_nuc_MeV, then the temperature in Sunny's one-zone test change by a few percent, more in line with expectations. But I do not know if this is the right thing to do, and haven't tested if it changes star sims