History value "Lnuc" is confusing (and maybe subtracting photodisintegrations twice)

[simonguichandut]

We recently found a situation where log_Lnuc in the history file went to -99 (meaning the non-log value is negative (?)), despite ongoing burning. A different quantity, power_nuc_burn, stays positive.

Looking through the code, it is a bit confusing what these quantities should mean. Here's the column descriptions:

mesa/star/defaults/history_columns.list

Lines 239 to 243 in ea1ca9b

	!log_power_nuc_burn ! total thermal power from all burning, excluding photodisintegrations
	log_LH ! log10 power_h_burn
	log_LHe ! log10 power_he_burn
	log_LZ ! log10 total burning power including LC, but excluding LH and LHe and photodisintegrations
	log_Lnuc ! log(LH + LHe + LZ)

My reading of these is that log_power_nuc_burn and log_Lnuc should be the same. Both are total burning power excluding photodisintegrations. But the actual definition of log_Lnuc is:

mesa/star/private/history.f90

Lines 2228 to 2230 in ea1ca9b

	case(h_log_Lnuc)
	power_photo = dot_product(s% dm(1:nz), s% eps_nuc_categories(iphoto, 1:nz)) / Lsun
	val = safe_log10(s% power_nuc_burn - power_photo)

If power_nuc_burn already excludes photodisintegrations, then this subtraction removes it an extra time (which could explain the -99 we were finding).

Then we go to the definition of power_nuc_burn:

mesa/star/private/report.f90

Lines 95 to 102 in ea1ca9b

	s% power_nuc_burn = 0d0
	do k=1,nz
	if (s% op_split_burn .and. s% T_start(k) >= s% op_split_burn_min_T) then
	eps_nuc = s% burn_avg_epsnuc(k)
	else
	eps_nuc = s% eps_nuc(k)
	end if
	s% power_nuc_burn = s% power_nuc_burn + eps_nuc*s% dm(k)

It seems to simply be $\int \epsilon_{\rm nuc}dm$. Does $\epsilon_{\rm nuc}$ include or exclude photodisintegrations?

Finally, adding to the confusion, there's this "ghost" history value (not listed in history_columns.list):

mesa/star/private/history.f90

Lines 2860 to 2861 in ea1ca9b

	case(h_lg_Lnuc)
	val = safe_log10(s% power_nuc_burn)

Notice this lg_Lnuc instead of log. This one does not subtract the photodisintegrations.

[Debraheem]

Hello, thanks for raising this! To clarify:

log_Lnuc = log(LH + LHe + LZ), and does not include photodisintegrations, hence the photodisintegrations are subtracted. There does not appear to be any double counting as s% power_nuc_burn (calculated from eps_nuc) includes photodisintegration, and therefore to remove it, it must be subtracted as is done in case(h_log_Lnuc).

You are correct that when Lnuc becomes negative, log_Lnuc becomes undefined .i.e. - 99.

I think you've also found a documentation bug for lg_Lnuc, thanks for catching that! I can add it back into the default history_column.list file with an approprate clarification on its definition. I can also add a clarification for log_Lnuc.

Perhaps this motivates providing a non-logarithmic power_nuc_burn value to the history_columns.list. Is that what you're looking for?

[Debraheem]

you can also use the following to compute it yourself

  !power_h_burn ! total thermal power from PP and CNO, excluding neutrinos (in Lsun units)
  !power_he_burn ! total thermal power from triple-alpha, excluding neutrinos (in Lsun units)
  !power_photo
  !power_z_burn

[simonguichandut]

If power_nuc_burn does include photodisintegrations (I agree, I think it does), then isn't this incorrect?

mesa/star/defaults/history_columns.list

Line 239 in ea1ca9b

!log_power_nuc_burn ! total thermal power from all burning, excluding photodisintegrations

lg_Lnuc (with lg) represents the same thing as power_nuc_burn, so adding it would only make things confusing, I think.

If you want to add a non-logarithmic value it should be Lnuc (=power_nuc_burn - photodisintegrations), since that is the one that can go negative.

[Debraheem]

"! total thermal power from all burning, excluding photodisintegrations "
This comment is incorrect. It includes photodisintegration, so there is no error, but the documentation needs corrected.

"lg_Lnuc (with lg) represents the same thing as power_nuc_burn, so adding it would only make things confusing, I think."
Yes, agreed.

"If you want to add a non-logarithmic value it should be Lnuc (=power_nuc_burn - photodisintegrations), since that is the one that can go negative."
Okay.

[Debraheem]

It appears that lg_Lnuc is the value displayed in the terminal, which could easily be confused with log_Lnuc.

         case(h_lg_Lnuc)
            val = safe_log10(s% power_nuc_burn)

[simonguichandut]

Thanks for having a look! Your PR clears up the confusion.