Hi all,

I'm trying to figure out how to get this as absolutely optimised as possible. This part of my code wasn't quite as critical before, but now I have optimised the rest of the code significantly, and this part now takes up 40% of the time when it used to take up like 10% of the time. So now I can justify spending more time on optimising it.

I've tried my best to get rid of unnecessary calculations, which helped a bit, but now I'm thinking I'd try to change the basic algorithm itself. However, I'm finding it hard to think up a more efficient way to do this, and it's a difficult topic to google, because I'm only coming up with graphical line-drawing algorithms, which aren't quite optimised for the same thing (especially in 2D!)

Basically, the set-up is I have a nested oct-tree of cells. That is, each cell is a square 2x2x2 block of subcells. Each of those subcells (or "nodes") can be empty (null), or itself contain another cell which can contain another 2x2x2 block of subcells and so on. Each of the cells may contain a list of particles. The idea is to propagate a ray through this tree structure, and build up a list of lists of particles for all of the cells and subcells that the ray intersects.

The basic algorithm is pretty much the standard line drawing algorithm, through each 2x2x2 block. I propagate through subcells one at a time, and if the subcell exists, I calculate the intersection with the subcell wall and pass it into the subcell for it to continue the line drawing on a deeper refinement level.

However, I feel like I may be doing more calculations than I need to here. A 2x2x2 block is a simple structure, and I'm wondering if there are simplifications that could be used to calculate which subcells are hit and what the points of intersection are, without having to tick forward one cell at a time. However, I can't think up anything that doesn't seem to involve just as much, or even more, calculations than the naïve algorithm.

Is there a more clever way to do this, or is this about as good as it gets? I'm attaching the relevant part of the code here, if that helps:

bool absorbtree::p_in_cell(double r[3]) {
    for ( int ii = 0 ; ii<3 ; ii   ) {
        if ( r[ii]<0. || r[ii]>=1. ) {
            return false;
        }
    }
    return true;
}


bool absorbtree::p_in_subcell(double r[3], int ipt[3]) {
    for ( int ii = 0 ; ii<3 ; ii   ) {
        if ( ipt[ii]==0 ) {
            if ( r[ii]<0. || r[ii]>=.5 ) return false;
        } else {
            if ( r[ii]<.5 || r[ii]>=1. ) return false;
        }
    }
    return true;
}

bool absorbtree::beam_subcells(double rline[3], double rtarget[3], double ddr[3], int idir[3], std::vector<int>** lol, unsigned int* lolsize, int ThisTask) {
    // add local particles
    if ( this->myP.size()>0 ) {
        lol[*lolsize] = &this->myP;
        (*lolsize)  ;
    }

    // beam through the 2x2x2 block on this level, check what subcells (if they exist) get hit
    int iline_local[3];
    bool ingrid,hitp;
    int dir_choice;
    double dt[3];
    int inode;

    for ( int ii=0 ; ii<3 ; ii   ) {
        if ( rline[ii]==.5 ) {
            if ( idir[ii]>0 ) {
                iline_local[ii] = 1;
            } else {
                iline_local[ii] = 0;
            }
        } else if ( rline[ii]>.5 ) {
            iline_local[ii] = 1;
        } else {
            iline_local[ii] = 0;
        }
    }


    hitp = false; // have I hit the destination particle?
    ingrid = true; // is the beam still in this cell?

    // Loop through line
    while(ingrid && !hitp) {
        // add this subcell, propagate to lower levels, check if r_incell at lower level

        inode=iline_local[0] (iline_local[1]<<1) (iline_local[2]<<2);

        // propagate line to subcell if it exists, otherwise just check if we hit that particle
        if ( this->nodes[inode] ) {
            double rline_out[3],rtarget_out[3];
            for ( int ii=0 ; ii<3 ; ii   ) {
                rline_out[ii] = rline[ii]*2-iline_local[ii];
                rtarget_out[ii] = rtarget[ii]*2-iline_local[ii];
            }
            hitp = this->nodes[inode]->beam_subcells(rline_out,rtarget_out,ddr,idir,lol,lolsize,ThisTask);
        } else {
            hitp = this->p_in_subcell(rtarget,iline_local);
        }

        // calculate "dt" to each edge
        // go in the direction with the lowest "dt"
        for ( int ii=0 ; ii<3 ; ii   ) {
            if ( idir[ii]==1 )         dt[ii] = ((iline_local[ii] 1)*.5-rline[ii])/ddr[ii]; // "right" side of this cell is "left" cell of next cell
            else if ( idir[ii]==-1 )   dt[ii] = (iline_local[ii]*.5-rline[ii])/ddr[ii]; // "left" side of this cell
            else dt[ii] = std::numeric_limits<double>::max();//(rtreei(iline[ii] idir[ii],ii,ilevel)-liner[ii])/ddr[ii];
        }
        dir_choice = 0;
        for ( int ii=1 ; ii<3 ; ii   ) {
            if ( dt[ii]<dt[dir_choice] ) {
                dir_choice = ii;
            }
        }

        // update line location - both the integer cell and the physical location
        iline_local[dir_choice] =idir[dir_choice];

        for ( int ii=0 ; ii<3 ; ii   ) {
            rline[ii] =dt[dir_choice]*ddr[ii];
        }

        // check if we've left the grid
        for ( int ii=0 ; ii<3 ; ii   ) {
            if ( iline_local[ii]<0 || iline_local[ii]>1 ) {
                ingrid = false;
            }
        }
    }

    if ( hitp ) {
        return true;
    } else {
        return this->p_in_cell(rtarget); // This shouldn't happen, but if because of cumulative errors we find that the particle is in this cell, but not in any of the subcells we checked, then we do want to stop
    }

}