Repository navigation
Rearrangement of struct grid #81
Description
I'd like to instigate a discussion on restructuring of the 'grid' struct. There are two reasons for doing so: (i) to allow precalculation of a 'vfac spectrum' for each link or edge connecting two grid points, which could speed up processing dramatically; (ii) to improve the present somewhat disorganized memory layout.
Most of the data which is stored in 'grid' falls naturally into several arrays of varying dimension and size as follows:
| Type | Name | Dimensions |
|---|---|---|
| int | id | par->ncell |
| double | x | par->ncell * DIM |
| double | vel | par->ncell * DIM |
| int | numNeigh | par->ncell |
| int | sink | par->ncell |
| int | nphot | par->ncell |
| int | conv | par->ncell |
| double | dopb | par->ncell |
| double | dens | par->ncell * par->collPart |
| double | t | par->ncell * 2 |
| double | nmol | par->ncell * par->nSpecies |
| double | abun | par->ncell * par->nSpecies |
| double | mol.dopb | par->ncell * par->nSpecies |
| double | mol.binv | par->ncell * par->nSpecies |
The elements I have omitted are of two types: those whose size scales with the number of nearest neighbours of the point, and those whose size depends on the molecular species. A list of the first set is
| Type | Name |
|---|---|
| point | dir |
| struct grid | neigh |
| double | w // doesn't seem to be used. |
| double | a0,a1,a2,a3,a4 |
| double | ds |
The species-dependent set is contained in the element 'mol' of type struct populations. Their sizes are tabulated as follows:
| Type | Name | Dimensions |
|---|---|---|
| double | mol.pops | par->ncell * par->nSpecies * m[i].nlev |
| double | mol.knu | par->ncell * par->nSpecies * m[i].nline |
| double | mol.dust | par->ncell * par->nSpecies * m[i].nline |
| double | mol.partner.up | par->ncell * par->nSpecies * m[i].npart * m[i].ntrans[j] |
| double | mol.partner.down | par->ncell * par->nSpecies * m[i].npart * m[i].ntrans[j] |
I discuss here better memory arrangements for the three sets of attributes in three respective sections.
A) Nearest-neighbour quantities.
The present arrangement is unsatisfactory because (i) information concerning the links between pairs of points is duplicated, (ii) allocation of **neigh is very messy and requires many small chunks spread across the memory. I suggest that these quantities should be contained in (or replaced by) two additional separate objects as follows.
Object 1:
Let's call this object 'links' for the sake of having a handle. Others might prefer 'edges', or some other name. It would have elements, or othewise include, the following:
| Type | Name | Dimensions |
|---|---|---|
| double | dir | numberOfLinks * DIM |
| double | dirNormalized | numberOfLinks * DIM // These two replace point *dir |
| double | aCoeffs | numberOfLinks * numberOfCoeffs // replaces *a0, *a1 etc. (numberOfCoeffs is, in effect, hardwired to 5 in the present code.) |
| double | ds | numberOfLinks |
| int | gridPointI | numberOfLinks * 2 // a new element which links back to grid. |
We'd possibly want additional elements to contain vfac information for example.
Object 2:
To avoid having scattered, variable-size chunks of memory, we delete the original element 'struct grid **neigh', and replace it by a 1D array of integer indices to 'links', which has a size equal to the sum of numNeigh for all grid points. I'll call this array 'nnLinkI' for the sake of the example. If a grid point has N nearest-neighbours, access to the links to these neighbours is via N consecutive entries in nnLinkI. We would also need to add an element
| Type | Name | Dimensions |
|---|---|---|
| int | firstNN | par->ncell |
to the grid struct. As an example of how this works, a loop which accessed the velocity-polynomial coefficient 0 for all the links of a grid point would look something like
int ni, gi;
double coeff;
struct grid g;
// blah
for(ni=0;ni<g.numNeigh[gi];ni++)
coeff = links.aCoeffs[nnLinkI[g.firstNN[gi]+ni]][0];
The nnLinkI array stores each link index twice, thus there is some wasted memory; but this is less than with the previous arrangement, and it is made up for by the much simpler, monolithic memory layout.
B) Species-specific quantities.
The natural way this information breaks up is to have a list of length par->nSpecies, each element of which was a struct of the following 2D arrays:
| Type | Name | Dimensions |
|---|---|---|
| double | pops | par->ncell * m[i].nlev |
| double | knu | par->ncell * m[i].nline |
| double | dust | par->ncell * m[i].nline |
Let's call this struct 'gridMolData' for the sake of the argument.
The remaining two elements of struct 'populations', partner.up and partner.down, are more problematic because one of their dimensions, ntrans, varies in size with the collision partner. Hence in the present code, partner.up and partner.down are, like neigh, malloced in numerous small non-contiguous chunks. Because there are not expected to be many collision partners (1, less usually 2, seem to be usual values), probably the best way to get around this problem is to calculate a maximum number of transitions per collision partner, then add further 3D array elements to gridMolData as follows:
| Type | Name | Dimensions |
|---|---|---|
| double | up | par->ncell * m[i].npart * m[i].maxNumCollTrans |
| double | down | par->ncell * m[i].npart * m[i].maxNumCollTrans |
There is potentially some wasted memory if .npart>1 and the several collision partners have very different numbers of transitions, but again, the risk seems worth the gain.
(Some parenthetical notes:
- Of course all the possible transitions between the species' energy levels apply in theory to each collision partner: however, only the in practice fairly limited subset of transitions for which rate data is known figure in these data structures.
- m[i].npart, which is set by the number of collision partners given in the moldatfile for that species, can differ between species, and they all can be different from par->collPart, which is set by the number of densities returned by the user-supplied density function. This potential mismatch probably warrants some careful consideration.)
C) Overall arrangement.
I had thought in terms of retaining the present arrangement in which an array (actually a pointer) of size par->ncell of grid objects is used to store the data. I had not thought further than getting rid of the duplicated link information and fixing the scatter-gun mallocs of element .neigh. Tuomas Lunttila however suggested instead to have a single struct of arrays (or pointers) for the grid, and I think I like this better. It is a subject for discussion however. The whole thing is of course complicated by the fact that many of the array dimensions are not known at the time of declaration of the data structure.
Activity
ParfenovS commentedon Mar 3, 2016
I think, it will can better to not store collision coefficients for every grid point. It requires a lot of memory especially in the case of molecules with high number of levels and high number of grid points. One can replace these arrays in the following way.
coll_fac stores the factor that can be used for the collision coefficients interpolation on temperature.
coll_fac = (g[id].t[0]-part[ipart].temp[tnint])/(part[ipart].temp[tnint+1]-part[ipart].temp[tnint]); (see line 276 in the molinit.c file of the main code branch).
coll_t1 and coll_t2 store indexes of collision temperatures (taken from the moldatfile) that can used for interpolation of collision coefficients in the getmatrix procedure (in stateq.c file).
I have a working version of the code without mol.partner.up and mol.partner.down arrays. This version requires about 50% less memory than the current code version in the case of model from the example directory. The overhead due to "on-the-fly" interpolation of collision coefficients is almost negligible because all necessary variables are precomputed and stored.
tlunttil commentedon Mar 16, 2016
I completely agree with @ParfenovS. The cost of on-the-fly temperature interpolation is negligible and storing the temperature-interpolated collision coefficients for each grid point separately just wastes memory. Actually I'd never realized that the code currently does that.
With some rearranging of stateq calculation we could probably do even without coll_fac, coll_t1, and coll_t2 arrays and calculate everything on-the-fly without significant performance penalty. The three arrays are quite small however, so it's not a major saving to remove them. By the way, @ParfenovS, why do you need to store two temperature indices?
Not sure if it's really necessary use (big) precalculated tables for mol.knu and mol.dust either.
ParfenovS commentedon Mar 16, 2016
For me, with two indices it was easier to take into account cases when the temperature in a grid point is outside the range given in the moldatafile. The goal was to not introduce if-statements into stateq due to on-the-fly interpolation.
tlunttil commentedon Mar 16, 2016
Ok, I see. But couldn't you take care of that in molinit? If the temperature in a grid point is lower than the minimum in moldatafile, set coll_t1=0, coll_fac=0.0, if the temperature is above the range in moldatafile, set coll_t1=ntemp-2, coll_fac=1.0. Then you could always interpolate between coll_t1 and coll_t1+1.
ParfenovS commentedon Mar 16, 2016
Yes, I'm taking care of most of the things in molinit. After that I'm doing inter/extrapolation in the getmatrix procedure (stateq.c file) with a single expression that looks like (in simplified form):
rate = C[coll_t1] + (C[coll_t2]-C[coll_t1])*coll_fac ,
where C - collision coefficient from moldatafile.
If one always interpolates between coll_t1 and coll_t1+1 and if the temperature is above the range in moldatafile then the expression above will return C[ntemp-2]. But I need C[ntemp-1]. So, I decided that using two indices is a simplest way.
tlunttil commentedon Mar 16, 2016
Sorry, I may have misunderstood what you're doing, but wouldn't coll_fac = 1.0 in your expression give rate=C[coll_t2] or, since we'd be using coll_t2=coll_t1+1, rate=C[ntemp-1] when the temperature is above the range? Or are you extrapolating beyond the temperature range given in moldatafile instead of clamping to the closest values?
ParfenovS commentedon Mar 16, 2016
No, you are right. The expression will return rate=C[ntemp-1] with coll_fac = 1.0.
I mostly use closest value. And, actually, this should be the default method so one can safely use only one index or maybe to not use any indices. But I also was interested (and, probably, will be interested) in extrapolatation beyond temperature range. One should consider my realization of the interpolation only as an example.
tlunttil commentedon Jun 7, 2016
@ParfenovS: Is your version without collision coefficient arrays ready to be PRd? I was thinking of writing a PR with that change myself, but I noticed that you wrote that you already have a working version. In my opinion this change is something that could be included in the next release; it's 'safe' in the sense that it shouldn't have any effect on results, and when using molecules that have a lot of collisional transitions the memory savings can be huge.
imckstewart commentedon Jun 7, 2016
@tlunttil do you think that on-the-fly interpolation could be done without a significant slow-down? That is a lot of extra operations.
There are several points in the code where a tradeoff between memory and speed can be made. I am not sure how best to deal with these. I suppose we could add the possibility to optimize the code for one or the other as the user prefers, but this is a lot of added complication.
tlunttil commentedon Jun 7, 2016
@imckstewart: I haven't tried, but I would be surprised if there were any significant effect to performance especially if we use the 'precalculated closest indices' method that @ParfenovS describes above. Also, currently
getmatrix()unnecessarily recalculates the collision coefficients for each iteration - they could be done only in the first iteration of the solver loop instateq()and then stored because unlike the radiative part, they don't change from one iteration to another.imckstewart commentedon Jun 7, 2016
I guess Sergey can tell us.
tlunttil commentedon Jun 7, 2016
I earlier said that arrays knu and dust could also possibly be replaced with on-the-fly interpolation, but that's harder to do without a performance hit. I didn't realise then that these arrays are needed for each step of raytracing (through
sourceFunc_cont()). Luckily these arrays only scale with the number of radiative transitions, so they're usually a lot smaller than the collision coefficient arrays.tlunttil commentedon Jun 7, 2016
I ended up trying this myself. With the default model a single-thread test took 7 m 57 s with on-the-fly interpolation and 7 m 59 s with the old version, so no difference in running time. Old version used ~90 MB of memory vs ~40 MB with interpolation. Results looked identical.
imckstewart commentedon Jun 7, 2016
Actually one member each of knu and dust is used in the raytracing: the one closest to the image frequency.
Re your test: great! Let's see a PR. 😄
tlunttil commentedon Jun 7, 2016
By raytracing I meant that thing in
photon(), which does need knu and dust for the frequency of each line, i.e. the whole array. It should be possible to replace even them with a (small) lookup+interpolation kludge, but I'd rather not do that (yet). It is possible that also knu and dust are very large arrays if the molecule happens to have a large number of allowed radiative transitions (in the worst case ~ 0.5*nlev^2), but at least none of the lamda molecules seems to have nearly that many.I need to tidy up the code a bit and make it more robust (for instance it currently doesn't work with more than one collision partner), but maybe tomorrow ;)
allegroLeiden commentedon Oct 10, 2016
This thread has been pretty quiet for a while, but I want to ping at least @tlunttil to see if he ever came to a conclusion about the best form for the new grid struct? I have been thinking about it some more and have changed my ideas a bit. The place I am at now can be summarized as follows:
par->pIntensity; the remainder malloc'd topar->ncell.linksor whatever to store information (either the olda0,a1etc or @tlunttil 's newv1,v2etc) which pertains purely to the edges/links between grid points, plus a 3rd auxiliary array to connect grid to edge and thus replaceneigh.tlunttil commentedon Oct 10, 2016
Sorry, I haven't thought about this for a while.