Renaming dust_density array to dust_density_o_n_grains for clarity

[cpinte]

Type of PR:
Bug fix / new physics / modification to existing code / other

Description:
Give a short description of the change.

Testing:
Describe how you have tested the change

Did you run the botscheck that the code and comments follow the code of conduct? yes/no

Did you update relevant documentation in the docs directory? yes/no

[gemini-code-assist]

Code Review

This pull request renames the dust_density variable to dust_density_o_n_grains across the codebase to clarify its definition as density divided by the number of grains. It also updates the Makefile to support HDF5 on Apple Silicon (Darwin arm64) and removes several debug print statements. The review feedback identifies multiple instances in density.f90 and SPH2mcfost.f90 where the new division by n_grains(k) could cause a division-by-zero error if a grain size bin is empty, suggesting the addition of safety checks to ensure numerical stability.

[gemini-code-assist]

Potential division by zero if n_grains(k) is zero. While n_grains is typically positive in MCFOST, it's safer to check for zero values to avoid runtime crashes or NaNs, especially since this division is a new addition in this PR.

                    if (n_grains(k) > 0.0_dp) then
                       dust_density_o_n_grains(k,cell_map(i,j,phik)) = sph_dens(i,jj,phik,1) / n_grains(k)
                    else
                       dust_density_o_n_grains(k,cell_map(i,j,phik)) = 0.0_dp
                    endif

[gemini-code-assist]

Potential division by zero if n_grains(k) is zero. A safety check is recommended here as well to handle potentially empty grain size bins.

                    if (n_grains(k) > 0.0_dp) then
                       dust_density_o_n_grains(k,cell_map(i,j,phik)) = sph_dens(i,jj,phik,n_a) / n_grains(k)
                    else
                       dust_density_o_n_grains(k,cell_map(i,j,phik)) = 0.0_dp
                    endif

[gemini-code-assist]

Potential division by zero if n_grains(k) is zero. A safety check should be added to handle cases where a grain size bin might be empty in the distribution.

                    if (n_grains(k) > 0.0_dp) then
                       dust_density_o_n_grains(k,cell_map(i,j,phik)) = &
                            ( sph_dens(i,jj,phik,l) + f * ( sph_dens(i,jj,phik,l+1) -  sph_dens(i,jj,phik,l) ) ) / n_grains(k)
                    else
                       dust_density_o_n_grains(k,cell_map(i,j,phik)) = 0.0_dp
                    endif

[gemini-code-assist]

Potential division by zero if n_grains(k) is zero. This new division should be protected with a check to ensure numerical stability.

                   if (n_grains(k) > 0.0_dp) then
                      dust_density_o_n_grains(k,icell) = rho_dust(1) / n_grains(k)
                   else
                      dust_density_o_n_grains(k,icell) = 0.0_dp
                   endif

[gemini-code-assist]

Potential division by zero if n_grains(k) is zero.

                   if (n_grains(k) > 0.0_dp) then
                      dust_density_o_n_grains(k,icell) = rho_dust(ndusttypes+1) / n_grains(k)
                   else
                      dust_density_o_n_grains(k,icell) = 0.0_dp
                   endif

[gemini-code-assist]

Potential division by zero if n_grains(k) is zero.

                   if (n_grains(k) > 0.0_dp) then
                      dust_density_o_n_grains(k,icell) = (rho_dust(l) + f * (rho_dust(l+1)  - rho_dust(l))) / n_grains(k)
                   else
                      dust_density_o_n_grains(k,icell) = 0.0_dp
                   endif

[cpinte]

/gemini review

[gemini-code-assist]

Code Review

This pull request refactors the codebase by renaming the dust_density variable to dust_density_o_n_grains to clarify its physical definition. It also introduces logic to automatically detect variable dust sizes in FITS and Phantom input files, and updates the Makefile for Apple Silicon HDF5 support. Review feedback identified a bug in src/density.f90 where uninitialized loop indices are used in an else branch, potentially leading to incorrect memory access. Additionally, the logic for reading dust types from Phantom HDF5 headers was found to be fragile and requires more robust handling of optional fields.

[gemini-code-assist]

The variables i, j, and phik are loop indices defined only within the if (n_grains(k) > 0.0_dp) block. In this else branch, they are either uninitialized or contain stale values from a previous iteration. Furthermore, using cell_map here only zeroes out a single cell, whereas the intention is likely to zero out the density for grain size k across all cells to avoid using uninitialized or invalid data.

        else
           dust_density_o_n_grains(k,:) = 0.0_dp
        endif

[gemini-code-assist]

The logic for determining ndusttypes from ndustsmall and ndustlarge is fragile. The got variable is overwritten by the second call, and if ndustsmall is present but ndustlarge is not, ndusttypes will be incorrectly set to 0. Additionally, ndustsmall and ndustlarge should be initialized to 0 to avoid using junk values if they are not found in the HDF5 file.

         ndustsmall = 0
         ndustlarge = 0
         call read_from_hdf5(ndustsmall,'ndustsmall',hdf5_group_id,got,ierr)
         call read_from_hdf5(ndustlarge,'ndustlarge',hdf5_group_id,got,ierr)
         ndusttypes = ndustsmall + ndustlarge