Suggestion for reading in the gasIIdust() in molinit.c

[cleevesli]

Hey all,

I use LIME to create line simulations from pre-tabulated chemical models, so it involves a lot of interpolation. The code has always been a bit slow even though I try to use fast interpolation techniques, and so finally last week I spent a little time trying to figure out why. The reason was actually pretty simple, in the kappa() function, the ordering of the for-loop essentially has the code retrieve the gas-to-dust mass ratio for every level in the molecule, which -- if gasIIdust is being interpolated -- ends up being one of the slowest steps.

As a suggestion, switching the ordering of the for-loop and then splitting out the 230 GHz "fix" as a second loop was a huge speed increase for me. It might be a second slower for those using a constant value of of gasIIdust, but for me it took a ~4 hour run down to < 20 minutes for a model with 7500 points.

Below is my suggestion. I can't think of any situation this wouldn't work for, gasIIdust shouldn't depend on molecular properties or frequency.

Best,
ilse

Original code in kappa in molinit.c:

  for(iline=0;iline<m[s].nline;iline++){
    for(id=0;id<par->ncell;id++){
      gasIIdust(g[id].x[0],g[id].x[1],g[id].x[2],&gtd);
      g[id].mol[s].knu[iline]=kappatab[iline]*2.4*AMU/gtd*g[id].dens[0];

      /* Check if input model supplies a dust temperature. Otherwise use the kinetic temperature. */
      if(g[id].t[1]==-1) {
        g[id].mol[s].dust[iline]=planckfunc(iline,g[id].t[0],m,s);
      } else {
        g[id].mol[s].dust[iline]=planckfunc(iline,g[id].t[1],m,s);
      }
    }
    /* Fix the normalization at 230GHz. */
    m[s].norm=planckfunc(0,par->tcmb,m,0);
    m[s].norminv=1./m[s].norm;
    if(par->tcmb>0.) m[s].cmb[iline]=planckfunc(iline,par->tcmb,m,s)/m[s].norm;
    else m[s].cmb[iline]=0.;
    m[s].local_cmb[iline]=planckfunc(iline,2.728,m,s)/m[s].norm;
  }

New version:

  for(id=0;id<par->ncell;id++){
    gasIIdust(g[id].x[0],g[id].x[1],g[id].x[2],&gtd);

    for(iline=0;iline<m[s].nline;iline++){
      g[id].mol[s].knu[iline]=kappatab[iline]*2.4*AMU/gtd*g[id].dens[0];
      /* Check if input model supplies a dust temperature. Otherwise use the kinetic temperature. */
      if(g[id].t[1]==-1) {
        g[id].mol[s].dust[iline]=planckfunc(iline,g[id].t[0],m,s);
      } else {
        g[id].mol[s].dust[iline]=planckfunc(iline,g[id].t[1],m,s);
      }
    }
  }

  // fix the normalization at 230GHz
  for(iline=0;iline<m[s].nline;iline++){
    m[s].norm=planckfunc(0,par->tcmb,m,0);
    m[s].norminv=1./m[s].norm;
    if(par->tcmb>0.) m[s].cmb[iline]=planckfunc(iline,par->tcmb,m,s)/m[s].norm;
    else m[s].cmb[iline]=0.;
    m[s].local_cmb[iline]=planckfunc(iline,2.728,m,s)/m[s].norm;
  }

[smaret]

Hi Ilse,

@allegroLeiden has made some extensive changes in this part of the code in preparation for the 1.6 release. Could you please check that this bottleneck is still present with the HEAD version of the code?

[tlunttil]

Ilse's version looks fine to me. In general the number of calls to user-supplied functions should be minimised in LIME because for all we know, they can be extremely slow. This has been done in the main iteration and ray tracing parts, but I guess no-one has looked closely at the initialisation phase because it (apart from qhull) usually takes a very small fraction of the total running time. In particular, I think all test cases have gasIIdust=constant.

I don't think any of the recent optimisations have touched this part of the code.

[smaret]

As a side note, computing chemistry on-the-fly at each grid point (i.e. within the model.c) would probably be faster (and more precise) than pre-compiling the chemistry on a fixed grid and interpolating the results at each grid point. This can be done quite easily with Astrochem (sorry for the shameless self advertisement :-)

[cleevesli]

:-) Maybe at some point in the future! At the moment my grid points aren't independent because of self-shielding... On the LIME end, I'm using a shepherd's method interpolation which has been pretty fast and reliable. I was surprised myself to see the "kappa" function being one of the largest time sinks, and happy it was such an easy fix. I only wish I'd found it sooner. :)

Since I have other bits and pieces running around I'd rather not muck up the source code, could someone else incorporate my fix in?

[tlunttil]

If we'll make a PR with this fix, we might as well check that there aren't other similar potential simple-to-fix time sinks left and do the changes all at once. I'll have a look.

[smaret]

I don't think any of the recent optimisations have touched this part of the code.

@tlunttil PR #128 has I think.

[tlunttil]

@tlunttil PR #128 has I think.

You're right, that code was moved from molinit.c to grid.c and the loop order was changed. Call to gasIIdust remains (needlessly) inside the inner loop, though, so the issue remains. But it would be a single-line fix.

[tlunttil]

gasIIdust was the only parameter with the "too many calls to user defined functions" problem, as far as I could see (after #145 is merged, which among other things changes how B is stored). Velocity still gets called a lot of times during ray tracing and density during grid construction, but these are unavoidable without doing things entirely differently.

Actually, maybe it would be easiest if @allegroLeiden makes the one line change to fix the problem when he merges #145 to master? I think he'll have to look at that part of code when doing the merge anyway, because it's different between master and #145.

[imckstewart]

Will do.

[allegroLeiden]

Actually, maybe it would be easiest if @allegroLeiden makes the one line change to fix the problem when he merges #145 to master?

Done.

[smaret]

This has been implemented with a7c30b8.