FPEs in `wd_cool_0.6M` from Skye

[warrickball]

I've got two FPEs in wd_cool_0.6M that are bubbling up from Skye. The top of the initial backtrace points to the first issue:

#1  0x110206a in __skye_ideal_MOD_compute_ideal_ele
        at /home/wball/mesa/mesahub/eos/make/helm_electron_positron.dek:3

The offending subroutine is

   subroutine compute_ideal_ele(temp_in, den_in, din_in, logtemp_in, logden_in, zbar, ytot1, ye, ht, s, e, p, adr_etaele, adr_xnefer, ierr)
      use helm_polynomials
      use eos_def
      real(dp), intent(in) :: temp_in, den_in, din_in, logtemp_in, logden_in, zbar, ytot1, ye
      type (Helm_Table), pointer, intent(in) :: ht

      ! Intermediates                                                                                                                                                                                                                                         
      real(dp) :: x, y, z, ww, zz, xni
      real(dp) :: temp, den, din, logtemp, logden
      integer :: k

      ...

      integer, intent(out) :: ierr

      temp = temp_in
      logtemp = logtemp_in
      den = den_in
      logden = logden_in
      din = din_in

      include 'helm_electron_positron.dek'
      ...

and line 3 of helm_electron_positron.dek is

        xne    = xni * zbar

zbar is an argument to the subroutine but xni hasn't been set. But nothing ever seems to use xne, so I commented out that line...

...which got me to the next FPE:

#1  0x1126ebb in __skye_coulomb_solid_MOD_solid_mixing_rule_correction                                                         
        at ../private/skye_coulomb_solid.f90:218               

where the offending line of code is in

      do i=1,n
         found = .false.

         do j=1,num_unique_charges
            if (unique_charges(j) == AZion(i)) then
               found = .true.
               found_index = j
            end if
         end do

         if (.not. found) then
            num_unique_charges = num_unique_charges + 1
            unique_charges(num_unique_charges) = AZion(i)
            charge_abundances(num_unique_charges) = AY(i)
         else
            charge_abundances(j) = charge_abundances(j) + AY(i) ! <— problem
         end if
      end do

charge_abundances is never initialised, so I think the FPE comes from adding something to the uninitialised data. Just looking at this loop, there are several possibilities of what the logic should be. I can compile MESA and run wd_cool_0.6M by adding charge_abundances = 0 before the loop but that might not be the intended behaviour. Another possibility is that the loop over j should exit when it finds the matching index, or that charge_abundances(j) should be charge_abundances(found_index).

[adamjermyn]

Nice sleuthing!

For the first FPE I'm going to copy-paste helm_electron_positron.dek into compute_ideal_ele and remove the offending line. HELM does actually use xne so we can't just comment it out uniformly, and I think it's better to have these actually diverge rather than having Skye continue pretending to be HELM.

The second FPE was an error in the problem line: it should be resolved by setting j -> found_index.

The fixes were simple enough that I've just pushed them to main, so I'm closing this issue. Feel free to reopen if you think it isn't addressed.