Conversation
The new member 'numRays' will allow diagnostic output from raytracing.
This has been done to prepare for adding a new raytrace algorithm.
This algorithm addresses a shortcoming of the existing algorithm in cases where there is a large number of grid points included in a single image pixel. The previous algorithm sampled such regions poorly. The effect of the new algorithm can be seen using the standard model file: the pixels towards the centre exhibit less variation from pixel to pixel. Selection of the new algorithm is mediated by a new (optional) parameter raytraceAlgorithm.
Its main idea has now been implemented in the new raytrace algorithm.
This function was made a little more readable and perhaps efficient, and renamed to write3Dfits() in preparation for adding a 2D version.
At present this is just called to write diagnostic images showing the number of rays per pixel, invoked by setting the image unit to the new value of 5. It could also be called to wrote continuum images, instead of writing a cube with just 1 layer as at present.
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.
|
It would surely be quite a bit of work. Besides, I think internally (precision, performance, etc.) qhull is quite good, it's just that its C interface is a 'bit' awkward (C++ interface looks a lot better, but there are also other options for convex hull/Voronoi/Delaunay computations in C++). |
|
Yes, well, probably Lime should have been written in c++ to begin with, but I am quite glad personally that it wasn't, because I hate c++. ;) qhull works when it works but I hate things which silently fail. |
|
The problem with qhull is also that it is slow and it can take most of the time when constructing a grid. It would be great if there will be parallelized (at least partially) realization of triangulation. |
|
I think with the HCKB thing I will wait til Sebastien gives me the go-ahead to merge, because some of these conversions have already been done in master. |
|
One very weird thing about qhull is how slowly it works in LIME. The equivalent command line programs are a lot faster. I even tried dumping points from LIME to disc and running command line qhull in case LIME grids were somehow especially hard to triangulate, but even then there was a huge difference (a factor of 100 or more I think). I'm not sure if LIME calls qhull in some completely wrong way or if the command line program didn't actually do what I thought it would. I even considered writing code to use qhull as an external program from LIME (which actually is what the documentation recommends in most cases). So I'd dump the point coordinates to a file, run qhull with system(), read facet dump from the output file, and reconstruct from that the needed data structures. Unfortunately writing such a piece of code didn't (doesn't) fell very exciting. However, I may do that anyway, because all the benchmarking stuff in #147 really needs a lot of points. |
Adding options in the Makefile may break compilation with some compilers. So I would rather keep plain ANSI-C (89?) until we implement a proper build system (with either cmake or autoconf/automake) than can detect which options each compiler supports. |
From @tlunttil comments it seems that this PR is ready to be merged, so please go ahead. |
|
@tlunttil Could you please open a new issue about the qhull performances issues? |
|
Strictly speaking, those _Bool types should be removed if we are to stick to C89, because they only came in C99. I'll open a new issue for qhull stuff. |
Ok - I don't get warnings with my older gcc but I guess that is because that version of gcc recognized _Bool, it just wasn't in the standard. I had not realized that. I can revert them all to ints if you like. |
|
gcc often supports (by default) some non-standard extensions. I guess that older gcc has a default language mode that is C89/C90 with some C99 (and other) extensions. Personally I prefer the booleans and C99, but I don't have very a strong opinion about this. I remember I already saw booleans in some other (already merged) PRs, but I didn't think anything of it then. |
|
Oh, great. 17 modules requiring a manual merge. "I'm just going out for a while; I may be some time." |
|
Well, the merged stuff runs, but the results are not right, so I guess I won't be uploading today. I am about to knock off anyway, it is stinking hot right under the roof here. (Why are astronomy departments always right under the roof? Parkes, Leicester, Jodrell, Cape Town, now here. Only in Bonn was I on the ground floor.) |
|
I can't find any more places where I could replace stuff by HCKB. |
It looks like the changes had already been done in master. |
|
At lines 218 and 428 of raytrace.c there is a multiplication by md[0].norminv. '0' is surely wrong here but to use the formally correct molecule would require some fiddly changes and needs some thought, which is why I have left it. The effect is probably small, but this serves to flag that we have not yet, formally speaking, got rid of all the traps for multi-species use (issue #21). |
|
I'm not sure what exactly is the meaning of these norm/norminv things. Any comments? |
|
I thought it was a linewidth correction but it appears to be correction for BB spectrum. (shrug) |
|
I'll have closer look, I'm currently very confused as to what they are trying to accomplish. |
|
I don't think norm/norminv are doing anything useful at all. norm is just a scaling that is applied in one place and then removed in another with norminv. Maybe they were supposed to prevent over/underflows, but they're not really needed nor correctly used for that anyway. |
- The raytrace() routine has been modified. Since lime-rt#145 a linear intra-cell interpolation has been performed within raytrace() for those image pixels which fell within a 2D projected cell significantly larger than the pixel dimensions. A vital part of this interpolation was the identification of the 2D cell in which the pixel fell. This has been until now performed by calling the qhull function qh_findbestfacet(). This function however was not returning accurate results for some pixels, which resulted in artifacts as reported in lime-rt#162. The call to qh_findbestfacet() has now been replaced with code which makes use of the raythrucells code to identify the relevant 2D cell. This seems to work ok. - I also changed the code in raytrace() which added projected points on a circle in the image plane of radius equal to the model radius. I had been using the point density averaged across the whole image plane to calculate the number of circle points needed, but because of the high density in the centre (for the default model anyway), this gave a number of circle points much higher than necessary. The routine has been modified so that only points in the outer 1/3 of the projected model are counted.
5 participants
This is intended to be the final word (at least in the framework of enhancements discussed for 1.6) on raytrace - or my final word anyway. Here I have merged and thus superseded the code which is in PRs #79, #96 and #109, also I have added new code which is a modified version of the fast raytrace algorithm which was in LIME 1.4 but which was (for various reasons) not propagated to 1.5. This PR thus answers the issue raised by @christianbrinch in #95 .