`SMOOTHSPEC` Call in `ADD_AGB_DUST` Dominates Run Time of `SSP_GEN`

[elijahmathews]

Context

When FSPS computes the SSPs, it calls GETSPEC to generate a stellar spectrum for each star on the grid (for each mass, metallicity, and age). Within GETSPEC, additional logic gets applied to AGB stars (presumably stars with phase == 4 or phase == 5, and provided that the AGB dust model is enabled with a nonzero dust value):

  !add circumstellar dust around AGB stars
  IF ((phase.EQ.4.OR.phase.EQ.5) &
       .AND.add_agb_dust_model.EQ.1.AND.pset%agb_dust.GT.tiny_number) THEN
     CALL ADD_AGB_DUST(pset%agb_dust,spec,mact,&
          logt,LOG10(lbol),logg,zlegend(pset%zmet),ffco,lmdot)
  ENDIF

What does ADD_AGB_DUST do? It essentially serves as a way to apply a dust screen (from DUSTY) on top of AGB stars' spectra. In the end, a dust transfer function dusty is computed, and this array is multiplied onto the stellar spectrum tspec to attenuate it.

However, there's a step immediately preceding this that complicates that a bit - see this set of lines at the very end of the ADD_AGB_DUST subroutine:

  !implement the dusty spectra (which are in units of 
  !flux_out/flux_in) into the AGB spectra
!  tspec = tspec * dusty
  CALL SMOOTHSPEC(spec_lambda,tspec,10000.d0,30000.d0,100000000.d0)
  tspec = tspec * dusty

Essentially, SMOOTHSPEC gets called before the dust is applied and smooths the stellar spectrum with a $\sigma_\textrm{vel} = 10^4\ \textrm{km}\ \textrm{s}^{-1}$ velocity smoothing kernel (for $3\ \mu\textrm{m} \leq \lambda \leq 10\ \textrm{mm}$), which almost certainly has the effect of smoothing the stellar spectrum from the MIR onwards into a nearly flat continuum.

Impact

Physical

Obviously, this process has the effect of completely wiping out any spectral features better than $R \approx 13$ in the spectral library (MILES/BaSeL/C3K) for AGB stars (phase == 4 or phase == 5). I can't quite see the reason for this (assume absorption line features have no chance making it through the dust screen, or we just don't trust the models?), but I suspect there may be an intentional reason for this (that SMOOTHSPEC call has been there for over 11 years at this point, untouched).

However, if this isn't intentional (shouldn't we be trusting the libraries?), then this line should probably be removed. In my pretty brief testing, it seems like this changes CSP fluxes on the order of ~4% or so, which is quite close to the standard 5% systematic uncertanties that are typically assumed when using FSPS (it may even contribute more than this - this was just for the parameter vector I happened to be testing with).

Computational

This is the largest impact, even if it only affects the SSP caching. During SSP caching, GETSPEC gets called for every single star in the grid, which means that SMOOTHSPEC (a pretty slow process on its own) also gets called for every single AGB star in the grid.

Profiling seems to indicate smoothspec takes up roughly ~90% of the compute time during SSP caching (when using MIST+MILES). In other words, the addition of SMOOTHSPEC in ADD_AGB_DUST alone slows down SSP generation by about a factor of 10.

Profiling Results

		     Call graph (explanation follows)


granularity: each sample hit covers 4 byte(s) for 0.06% of 17.59 seconds

index % time    self  children    called     name
                0.00   17.54       1/1           main [2]
[1]     99.7    0.00   17.54       1         MAIN__ [1]
                0.29   15.99       2/2           ssp_gen_ [3]
                0.18    0.90       1/1           sps_setup_ [9]
                0.00    0.17       1/1           compsp_ [11]
                0.01    0.00     323/323         check_array_1d.5 [25]
                0.00    0.00     535/535         check_val.3 [36]
                0.00    0.00       2/2           __fsps_context_MOD_fsps_context_set_pset [92]
                0.00    0.00       1/1           __fsps_context_MOD_fsps_context_create [98]
                0.00    0.00       1/2           __fsps_context_types_MOD_fsps_context_state_destroy [93]
                0.00    0.00       1/2           __sps_utils_MOD_sps_takedown [95]
-----------------------------------------------
                                                 <spontaneous>
[2]     99.7    0.00   17.54                 main [2]
                0.00   17.54       1/1           MAIN__ [1]
-----------------------------------------------
                0.29   15.99       2/2           MAIN__ [1]
[3]     92.6    0.29   15.99       2         ssp_gen_ [3]
                1.18   14.79  211116/211116      getspec_ [4]
                0.00    0.02     214/214         __fsps_imf_MOD_compute_imf_weights [20]
                0.00    0.00     214/214         add_remnants_ [28]
                0.00    0.00     214/214         mod_gb_ [38]
-----------------------------------------------
                1.18   14.79  211116/211116      ssp_gen_ [3]
[4]     90.8    1.18   14.79  211116         getspec_ [4]
                0.30   14.47   79300/79300       __fsps_dust_MOD_apply_agb_dust_screen [5]
                0.02    0.00  346228/96509046     __fsps_interpolation_MOD_find_interval [7]
-----------------------------------------------
                0.30   14.47   79300/79300       getspec_ [4]
[5]     84.0    0.30   14.47   79300         __fsps_dust_MOD_apply_agb_dust_screen [5]
                6.69    7.78   79300/79300       smoothspec_ [6]
                0.01    0.00  158600/96509046     __fsps_interpolation_MOD_find_interval [7]
                0.00    0.00   79300/79300       __fsps_dust_MOD_compute_circumstellar_optical_depth [89]
-----------------------------------------------
                6.69    7.78   79300/79300       __fsps_dust_MOD_apply_agb_dust_screen [5]
[6]     82.2    6.69    7.78   79300         smoothspec_ [6]
                4.09    0.00 80965300/96509046     __fsps_interpolation_MOD_find_interval [7]
                3.69    0.00 161772000/161790800     __fsps_integration_MOD_integrate_trapezoid_array [8]

		         ...cut for brevity...

Note: The profiling done above was done with my heavily-modified fork of FSPS, but in this fork __fsps_interpolation_MOD_find_interval is LOCATE and __fsps_integration_MOD_integrate_trapezoid_array is TSUM. None of my modifications should be affecting this particular slowdown.

In the call graph above, we see that SSP_GEN (ref [3]) and its children consume about $0.29\ \textrm{s}+15.99\ \textrm{s}=16.28\ \textrm{s}$ of total run time. One of its children is SMOOTHSPEC (ref [6]), which takes a total of $6.69\ \textrm{s}+7.78\ \textrm{s}=14.47\ \textrm{s}$ of run time. Thus, SMOOTHSPEC accounts for about ~90% of the total run time of SSP_GEN.

In summary:

[3] ssp_gen_ (100.0%)
 └── [4] getspec_ (98.1%)
      └── [5] apply_agb_dust_screen (90.7%) # a.k.a. ADD_AGB_DUST
           └── [6] smoothspec_ (88.9%)
                ├── [7] find_interval (25.1%) # a.k.a. LOCATE
                └── [8] integrate_trapezoid_array (22.7%) # a.k.a. TSUM

Proposed Changes

1. Optimization (If Smoothing Is Intentional)

If the physical requirement is that AGB spectra must be featureless continua in the MIR (to match DUSTY resolution or simulate optical depth or whatever), we should replace the expensive SMOOTHSPEC call with a faster approximation.

• Proposal: Replace the convolution with a generated Blackbody spectrum or a spline fit for $\lambda > 3\ \mu\textrm{m}$.
• Benefit: This would eliminate the convolution overhead, reclaiming the ~90% runtime cost in SSP_GEN.

2. Removal (If Smoothing Is Unnecessary)

If modern DUSTY grids or current physical assumptions no longer require the underlying star to be smoothed to $10,000\ \rm{km\ s^{-1}}$, we should simply remove the SMOOTHSPEC call.

• Benefit: Immediate 10x speedup in SSP generation and more accurate recovery of AGB molecular features in the MIR.