#gradientbundleanalysis

GradientBundleAnalysis

Point-wise integration

We finished and implemented the point-wise integration. It works as described below:

Volume zones to be integrated must exist as closed, finite element surfaces

Simply loop over every cell (cubic cell composed of 8 nodes) in the full volume, determining if the cell is

contained entirely outside the FE volume, which is performed with an additional time-saving step of disqualifying cells that are outside of the bounding cube of the FE volume;

contained entirely within the FE volume; or

partially contained within the FE volume

If cell is entirely within, add the variable value (i.e. rho) for the center of the cell to the integral

if cell is partially contained in FE volume, subsample the cell according to a specified resolution, each time adding the variable value for the centers of the subcells if they are interior to the FE volume

This method works, but it takes a very high subsampling resolution to get decent accuracy, resulting is very, very long run-times.

FE volume integration

I have another idea for a much faster integration scheme that has the additional benefit of lacking any kind of user-specified resolution.

The method is pretty simple: Take the closed, FE surfaces and turn them into FE volumes composed of polyhedra. Once the volumes are created, Tecplot 360’s built-in integration can be used, although I would prefer to also implement my own integrator to prevent reliance on the Tecplot 360 Integration add-on.

My basic idea for the algorithm to make a FE volume out of the closed FE surface is as follows:

Get a complete list (IsInterior) of all Full Volume Zone nodes that are interior to the FE surface. The list will be in the form of a vector<bool> of length iMax • jMax • kMax so that lookup time is constant.

Make another list (ElemDirCreated) of vector<vector<bool> >, also of length iMax • jMax • kMax, but only the elements i for which IsInterior[i] is true will have their inner vector<bool> resized to length 8. These 8 bools will be used to keep track of whether FE volume elements have been created for the 8 directions from each interior node.

With each interior node, check it’s 8 values in ElemDirCreated.

If the element has been created, move on

If the element hasn’t been created, then check the interior-ness of the node’s 26 neighbors.

If a neighbor is interior, then simply use that node for a node of the new element.

If a neighbor is not interior, then use the intersection point of the line connecting the node and its neighbor as a node in the new element.

Need to think about this part of the algorithm carefully such that cases where only one corner of a cell is cutoff by a surface result in three nodes on the surface corresponding to the intersections between the surface and the lines between the one exterior node and its three neighbors in the cell.

Since node and element numbers serve no purpose other than to provide information about the connectivity of the FE volume, these can simply be generated as the elements are built.

Need to look into the possibility of using a set number of sides for the FE volume, as this will allow the structure of the connectivity list to be predicted.

Fully interior cells will be hexahedral, but…

I just realized this won’t work…

This method, and the current point-wise integration method, won’t work because of the assumption that, if all 8 nodes of a cell are not contained within a volume of interest, then the cell does not contribute to the volume’s density. This assumes that there will never be a skinny volume shooting through the center of the cell, which is a bad assumption.

I need to use a method of creating FE volumes that starts by using the existing information known about the volume. It will likely be impossibly non-trivial to utilize the discreet data in the system. Instead, I’ll probably devise a method that uses nothing but interpolated data, which isn’t that bad anyways.

Side note: Max and min number of sides for a hexahedron coincident with N planes is 6 + N and 6 - N respectively1

via https://dayone.me/29GPzC6

That’s the best thing for me to remember when coding. I found the problem with the compiler optimizations: I turned them on in addition to making changes to the code that killed the performance. I neglected to test the new code before messing with the optimization settings, and, like a moron, I thought the compiler was to blame…

I had setup some bad logic for a counter that was supposed to make gradient paths finish after stalling for 10 iterations, but instead the counter was never allowed to actually iterate. It would iterate and then get reset to 0 because the gradient path hadn’t failed yet. I fixed that and now the code is lightning fast again, and with all the optimizations enabled I get a smaller file; just over 400KB for #GradientBundleAnalysis. Pretty cool, minus my juvenile rage.

#GradientBundleAnalysis

Parallel version is finished and seems to work great. With compiler optimizations enabled, the speedup vs the serial version for mesh levels 3 and 4 are 22 and 26 respectively. That’s right, speedup of 26! The parallel version does much less work, uses my cheaper gradient paths, and get’s direct access to dataset values, in addition to running the bulk of the code in parallel.

There was a funny bug in some systems, where the gradient path will stall for a few iterations, so I made it so the path has to stall more than some number of iterations before it assumes that it’s fallen into a cage. When run on N2, the bug resulted in a sweet-looking effect that looked like a mushroom cloud. I fixed it though so there’s no more mushroom clouds (yay!).

#GradientBundleAnalysis bugs

Here’s a quick shot showing the problems with 3- and 4-sided FV volumes (that’s right, they’re both screwed up…)

#GradientBundleAnalysis

Making great progress. Almost finished with the parallel version, which is way faster than the serial version for several reasons. Right now I’m just ironing out some bugs (mmm, burnt bugs…) related to FE volume creation. Granted, they’re annoying bugs whose cause I have no clue, but at least I know where this is all happening.

#Mathematica

I took some code that Amanda had made for comparing stream traces and turned it into an easy-to-use module so she can work faster (see pic). I’ll make another version of the tool that works for arbitrary stream traces too so that Erica can perform the same type of analysis.

For the arbitrary case, you’ll lose information about the stream traces in order to plot 3-d stuff in a 1-d plot. I’m thinking of replacing the x-axis with distance from a plane normal to the inter-CP axis, and the y-axis with distance from the inter-CP axis.

Grad path

I’ve now finished all the necessary methods for grad paths in GBA. Now I need to make a real default constructor (that takes no arguments) and a SetAll method that uses the current default constructor to set everything.

After that the only code remaining is for making FE volume zones from sets of grad paths. Here I’m not sure if I should make a FE zone class, but I think I won’t right now. A FE zone class in GBA would likely only be useful in GBA, and would become a two-face class once I tried to use it in Bondalyzer. Maybe not though. It’d be nice to be able to just pass a vector of grad paths to a FE zone class and have it put them together into a nice surface. Then pass a bunch of FE zone objects to another FE zone and have it make a nice volume out of it. I’ll mull it over, but for now I’ll just do that FE zone creation part in the main function.

Grad paths

Got the GSL ODE solver-based grad path working ok. Added more functionality than will be necessary with GBA, but I'm planning on replacing Bondalyzer's grad paths with this eventually, so the functionality will be used.

Next step is to add methods for resampling and resample/concatenation with other grad paths. For concatenation, need to add another constructor that simply takes pointers to a Bondalyzer grad path zone and makes a CSM_GradPath_c out of it, then concatenation can just be between CSM_GradPath_c’s.

Parallel

Something I hadn’t explicitly though about: When you fetch a raw pointer, Tec360 must read in ALL the data to memory. In some cases this can flood the memory and then you’re stuck waiting on your OS to write it all out to virtual memory. Need to keep this in mind when Bondalyzer goes parallel.

For GBA this will be an issue with larger systems, since 7 pointers (3 for XYZ, 3 for Gradient, and 1 for Rho) are necessary in order to do a grad path. Additionally, there’ll be NumGradPaths X PointsPerGradPath floats necessary. Still think it’ll be way faster than the serial version, and have a smaller memory footprint than my first idea at parallel with stream traces.

Making progress. Today I wrote the function that the GSL ODE solver will use to make the gradient path. It needs to interpolate between grid points so that the solver can query arbitrary locations. So the function takes the location, figures out the indices of the corners of the cell that contains it, and interpolates from those corners to get a value for gradient at the location.

Tomorrow I’ll write the actual grad path function that uses the ODE solver. Right now It’ll be a simple version just to see if everything’s working right. Later I’ll add more robustness like the different use cases (terminating at a CP, an arbitrary point, a value of Rho, or at the system boundary).

I’m hoping I can use the ODE driver object that does all the ODE solver stuff for me. I should just be able to take small enough time steps so that I can see if I’m approaching the end of the system (or other termination criteria) and handle everything myself.

I’m a little concerned that, in the even of a collision with the system boundaries, that I won’t be able to have the ODE solver tell be the new position, should that position be out of bounds. I need the new position so that the last step still reflects the behavior of the system. On the other hand, the behavior at the boundary is usually boring anyways, and it’ll probably be just as well to go in the same direction as the previous step and interpolate that direction until the boundary collision.

Another thing I’d like to keep in mind is how to use the GSL ODE solver with periodic conditions such that the gradient path can pac-man from one side of the system to the other. Again, shouldn’t be such a big deal, but it’ll freak out the solver since the coordinates will be changing drastically. Another option is to let the solver continue moving in the same direction, with locations that are farther and farther “outside” the boundaries. Then the function I provide the solver can just deal with those crazy locations, and I’ll use the periodicity of the system to take the resulting locations and put them in the correct place in the actual system. So the solver gets to continue in the same direction at all times, and I do the work of putting the results in the actual system outside of the solver.