moddump_nucleationtogrowth

[AlagueroA]

Description:

Added a moddump utility that starts a dust-growth simulation from a dust nucleation dump. Minor modifications have been made in growth.f90 and dust_formation.f90 for that purpose.

Components modified:

• Setup (src/setup)
• [ x] Main code (src/main)
• [ x] Moddump utilities (src/utils/moddump)
• Analysis utilities (src/utils/analysis)
• Test suite (src/tests)
• Documentation (docs/)
• Build/CI (build/ or github actions)

Type of change:

• Bug fix
• Physics improvements
• [ x] Better initial conditions
• Performance improvements
• Documentation update
• Better testing
• Code cleanup / refactor
• Other (please describe)

Provide the dust-nucleation conditions as the initial conditions for a dust-growth simulation.

Testing:

Did you run the bots? no

Did you update relevant documentation in the docs directory? no

Did you add comments such that the purpose of the code is understandable? yes

Is there a unit test that could be added for this feature/bug? no

If so, please describe what a unit test might check:

Create a dust nucleation file, use moddump_nucleationtogrowth.f90 on that file, and check that the code runs for one timestep using the created dump.

Related issues: #

[gemini-code-assist]

Code Review

This pull request introduces support for converting dust nucleation data to dust growth data, including updating growth rate calculations with a dynamical timescale and adding a new utility moddump_nucleationtogrowth.f90. The review feedback highlights a critical bug in the new utility where a scalar dust_to_gas ratio is overwritten in a loop and incorrectly used for scaling, suggesting an average calculation instead. Additionally, feedback recommends adding defensive checks for the allocation of the nucleation array, simplifying the dynamical timescale calculation, and protecting against potential division-by-zero and out-of-bounds errors.

Important

The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.

[gemini-code-assist]

In the loop over gas particles, dust_to_gas is a scalar variable that gets overwritten for each particle. Consequently, after the loop, massoftype(igas) is scaled using only the dust_to_gas ratio of the very last particle (np_gas), which is incorrect and can lead to wrong total gas mass scaling (especially if the last particle has zero dust). We should compute the average dust_to_gas ratio over all gas particles and use that to scale massoftype(igas) when ndusttypes == 1.

  dust_to_gas_avg = 0.
  do i=1,np_gas
     if (ndusttypes > 1) then
        dustfrac(1:ndusttypes,i) = dust_to_gas*dustbinfrac(1:ndusttypes)
     else
        if (nucleation(5,i)<tiny(nucleation(5,i))) then
            dust_to_gas = tiny(dust_to_gas)
        else
            dust_to_gas = nucleation(5,i) * mc / mass_per_H
        endif
        dust_to_gas_avg = dust_to_gas_avg + dust_to_gas
        call set_dustfrac(dust_to_gas,dustfrac(:,i))
     endif
  enddo

  if (ndusttypes > 1) then
     massoftype(igas) = massoftype(igas)*(1. + dust_to_gas)
  else
     dust_to_gas_avg = dust_to_gas_avg / np_gas
     massoftype(igas) = massoftype(igas)*(1. + dust_to_gas_avg)
  endif

[gemini-code-assist]

Add a defensive check to ensure that the nucleation array is allocated before accessing it, preventing potential segmentation faults if this utility is run on a dump file without nucleation data.

 excitation_HI = 1

 if (.not. allocated(nucleation)) then
    print *, "Error: Nucleation array is not allocated. This dump does not contain nucleation data."
    stop
 endif

[gemini-code-assist]

In get_dynamical_timescale, the calculation 1/(xyzmh_ptmass(ivwind,1)/xyzmh_ptmass(iReff,1)) can be simplified to xyzmh_ptmass(iReff,1)/xyzmh_ptmass(ivwind,1) to avoid an extra division. Additionally, we should protect against out-of-bounds array access if nptmass < 1 (by importing and checking nptmass), and prevent division by zero if xyzmh_ptmass(ivwind,1) is zero. Finally, the intent of tdyn should be out rather than inout since it is only written to.

 subroutine get_dynamical_timescale(i,tdyn)
 use part,           only:xyzmh_ptmass,ivwind,iReff,Omega_k,nptmass
 use dim,            only:do_nucleation
 integer,   intent(in) :: i
 real,      intent(out) :: tdyn
 
 if (do_nucleation .and. nptmass >= 1) then
    if (xyzmh_ptmass(ivwind,1) > 0.0) then
       tdyn = xyzmh_ptmass(iReff,1)/xyzmh_ptmass(ivwind,1)
    else
       tdyn = 1/Omega_k(i)
    endif
 else
    tdyn = 1/Omega_k(i)
 endif

 end subroutine get_dynamical_timescale

[gemini-code-assist]

Add a declaration for dust_to_gas_avg to support calculating the average dust-to-gas ratio across all gas particles.

 real    :: dust_to_gas,smincgs,smaxcgs,sindex,dustbinfrac(maxdusttypes),udens,dust_to_gas_avg