Conversation
This is the first of many commits I plan for an eventual pull request which will provide LIME with a smoother raytracing capability. This added functionality requires quite a lot of new code and will also involve some cleaning up of existing code. The smoother raytracing algorithm works by linear interpolation of level populations inside the Delaunay cells of the model grid. This first commit lays the groundwork for reading and using the cell information by defining several new structs to allow convenient storage and manipulation of this data. Structs added were: cell, pop2, gridInterp, gAuxType, intersectType, faceType and triangle2D. None of these are yet active so as yet they cause no change in the task function. Also added in this commit is a new user parameter traceRayAlgorithm. This will be used to invoke the new raytracing algorithm; at present it has no effect. Finally, I have commented out the function template for FastExp() in lime.h and added the variable type of the input argument of greetings_parallel() in messages.c just to avoid the warning messages that gcc issues now on compilation. This is just for present convenience, I plan to make these minor fixes the subject of a separate PR.
This attribute belongs in struct populations rather than struct grid (directly). Moving it involved making changes to 18 lines of code over 7 modules. I checked afterward that the code runs still ok with the default model. I also slightly rearranged attributes x and vel in the struct template for struct grid, and added a new attribute (as yet unused) B. These changes have no further effects on the code.
This is the first of several commits which make changes to the function qhull(). The only one of these which introduces new functionality will be an option to read and return the information about the Delaunay cells. However I plan to use this necessity to reform the whole function, starting with the name. 'qhull' is a bad name because it is not immediately explanatory and because it can be confused with the library of that name.
I replaced the input argument par with two, namely numDims and numPoints. This helps to make the function more generic. An additional variable ppi of type matching numPoints was added to loop over the latter quantity. Iteration variables of types other than the upper limit variable can cause strange errors!
The block starting at line 168 in grid.c is rearranged, the if test being inverted such that the code exits if qh_new_qhull() returns a non-zero status. This involves unindenting all the code (~32 lines) which was in the block before. No other functional changes have been made to these lines in the present commit.
The array 'simplex' is renamed to 'pointIdsThisFacet' for added clarity. Additional variables idI and idJ have been added. All three of these have been made type unsigned long, and variable 'id' is also now of the same type.
Three more arguments, to mediate the return of the cell information. These are not yet used, but the changed interface required some additions to every function which calls delaunay().
This is not yet called anywhere however.
This is a purely cosmetic change, but a useful one, because searching for g returns too many false positives. I plan to do this elsewhere too.
I also rearranged the variable declarations a little. The only substantive change is the explicit cast of the return value of qh_pointid() to unsigned long. This is the last of the delaunay-focussed commits.
To assist the interpolative raytracing algorithm I plan to introduce a function calcLineAmpInterp() to return the spectral line amplitude relative to line centre, a quantity which in LIME is usually signified by the variable named 'vfac'. This function is therefore analogous to velocityspline2() in raytrace.c and the similarly-named functions in photon.c. This planned new addition is a good excuse to change the names and interfaces of all three existing 'velocityspline' functions to be more similar. In the present commit I have changed the names as follows: velocityspline to calcLineAmpSpline velocityspline_lin to calcLineAmpLinear velocityspline2 to calcLineAmpSample I also made some of the arguments 'const'. This has no effect on the functioning, it is just more conservative practice. The final change is the addition of a calcLineAmpLinear() template to lime.h
This contains the code for finding the cells traversed by a ray on the one hand and interpolating within Delaunay cells on the other. None of this code is yet used, but it compiles at this stage.
This function is called by sourceFunc_pol() when the user parameter par->polarization is set. In the old version of stokesangles(), the user-supplied function magfield() was called for an arbitrary position in the model. In the planned interpolation raytracing scheme however it is possible to interpolate the B field between the vertices of any Delaunay tetrahedron. It is also desirable to remove the calls to user-supplied functions from any sections of code traversed in the case that either par.doPregrid or par.restart are set (although this is not yet achieved within the alterative raytrace code, as will be seen when this is added). For these reasons I re-wrote stokesangles() to accept a precalculated value of the B field as an input argument. These values are now also calculated via calls to magfield() at the grid cell loci within grid.c. At some point the storage of grid data needs to be improved, i.e. to include both velocity and B-field values. Until that time the par.doPregrid option will not work, not even after the present commit, even though it moves things in the right direction.
The changes are as follows: . Changes to the arguments, to make it easier to call this routine from the planned interpolated raytrace; the change from 'struct grid *g' to 'struct pop2 gm' also allows much more compact and clear expressions. . Variable name 'angles' changed to more descriptive 'trigFuncs'. . All elements of snu[3] are now explicitly defaulted to zero. The interface changes have necessitated changes in the call to sourceFunc_pol() within traceray() and the addition of a new argument gAux (whose values are set within raytrace()) to the traceray() call. I removed raytrace_1.4() at this point. This obsolete chunk of code was supposed to have already been removed in a previous PR.
I just changed sourceFunc_line() and sourceFunc_cont() to bring them more into line with sourceFunc_pol(). This necessitated changes in photon(), getjbar() and traceray(). In the last of these, the call to sourceFunc_cont() is slightly rearranged, with new variables contMolI and contLineI being set before the start of the loop.
These duplicate the functionality of sourceFunc_cont() and sourceFunc_line() but are made necessary, or at least desirable, by the present poor way in which grid data is stored in memory structures. When this is fixed (which it will hopefully be eventually), the *_raytrace versions can be merged into the previous ones.
There is a test near the start of the function to see whether the ray will intersect the model sphere at all; if so, the following large block of ~110 statements (containing the entire functionality of the function) is entered. It is obviously much less spendthrift of indents to simply exit the function if the ray does _not_ pass through the model sphere. The code has been rearranged to implement this change. All the statements which were formerly within the block have been unindented by 2 spaces; there are no further changes in the present commit.
For better searchability.
The substantive changes are: . Where par->polarization and FASTEXP was defined, one of the two necessary exponentials was still calculated by the native exp() function. This has been changed now such that Fastexp() is called for both. . sourceFunc_cont_raytrace() is no longer (unnecessarily) called once per channel, but only once; the results being stored and added to jnu and alpha for each new channel.
I also changed the interface of traceray(). This was inadvertent and should have been in a separate commit.
In the previous raytracing algorithm the values of the sink point atributes are not accessed. However they are accessed in the interpolation routine. Because of this it is necessary to be a bit more careful about setting them. Thus I have set their: . nmol to zero in LTE(); . pops to L.T.E with the CMB in molinit(); . vel to the user-specified values (instead of zero) in getVelosplines().
With this commit, the new algorithm is essentially complete. I'll make 1 more with some minor cosmetic changes.
Here I have: . Improved some variable names in raytrace(). . Ordered all the function templates alphabetically (or at least, more so) in lime.h. . Removed some tab indents in predefgrid.c which were bugging me.
|
Do you know how much slower approximately this new raytracing is compared to version 1.5? |
|
No because the one time I tried 1.5, I gave up on the raytracing and killed it before it finished. The same model took 20 seconds to run in 1.4. Sent from my iPhone On 26. maj 2016, at 20.22, tlunttil <notifications@github.commailto:notifications@github.com> wrote: Do you know how much slower approximately this new raytracing is compared to version 1.5? — |
Do you mean release version 1.5 or one with the changes in this PR included? I haven't tried this PR yet, but raytracing in release 1.5 shouldn't be that terribly slow (though definitely slower than in 1.4 if the image size is large). |
|
The release version I suppose. Sent from my iPhone On 26. maj 2016, at 22.53, tlunttil <notifications@github.commailto:notifications@github.com> wrote: No because the one time I tried 1.5, I gave up on the raytracing and killed it before it finished. The same model took 20 seconds to run in 1.4. Do you mean release version 1.5 or one with the changes in this PR included? I haven't tried this PR yet, but raytracing in release 1.5 shouldn't be that terribly slow (though definitely slower than in 1.4 if the image size is large). — |
|
Hmm, that sounds strange. I wonder if it's a bug rather than a 'feature' in version 1.5 raytracing. I don't don't think a model that only takes 20 seconds in 1.4 should take very long even in 1.5. Can you send the model files so I could have a look? |
|
There seems to be a bit of a non-sequitur between @tlunttil 's first comment (which refers to the PR) and @christianbrinch 's reply to that (which refers to 1.5). The probable reason why 1.5 is slower has been remarked on already by @ParfenovS in thread #95. On my desktop, for the master branch, which is still essentially 1.5, the time to run the standard model is about 13 min. Raytracing takes 25 sec of this. For branch 'raytrace_smooth' (the source of the present PR) with |
|
Yeah, what I originally wanted to know is what you just told. Ok, so raytracing in this PR is a factor of 5-ish slower than release in 1.5. That doesn't sound too bad and it's still only a small fraction of the total run time. I have some further comments and questions about 1.4 vs 1.5 raytracing, but it's better to continue that discussion in #95. |
I had left a comment in the part of traceray() which dealt with polarized continuum, but this comment only indicated that I did not understand the code. I do now, so I have removed the comment.
|
Superseded by #145. |
5 participants
The current raytracing algorithm leads to step-function artefacts in the channel maps it produces. These occur because the algorithm, although continuously sampling velocity at relatively finely-spaced intervals along the path of a ray, uses a single set of values for level populations throughout an entire Voronoi cell. In outer regions of an image where the Voronoi cells are likely to be significantly bigger than the image pixels, this leads to the step artefacts which can be seen in the image which is generated by the default example model file. The algorithm subject of the present PR removes these step artefacts by interpolating (to 1st order) the level populations within each Delaunay tetrahedron. This new algorithm has necessitated a lot of extra or altered code, from a new qhull() function (renamed to 'delaunay') which can return the necessary information about the tetrahedra, through routines to work out the sequence of cells passed through by any ray, to the new traceray_smooth() function to actually implement the algorithm.
To aid in checking these many changes I have divided the push into many smaller commits. I checked the intermediate code at several points with valgrind. The final code passes valgrind when the new algorithm is invoked, is backward compatible (because it defaults to the old), and with the old algorithm selected, still produces an identical cube to before.
The following points are worth noting: