New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Nearest neighbor particle interpolator #1527
Nearest neighbor particle interpolator #1527
Conversation
1a9e096
to
1b8e4b7
Compare
/run-tests |
This seems like something for @gassmoeller :-) |
I was on it yesterday, but got sidetracked again and again. Will try to finish the review now. @jperryhouts you will need to rebase to current master at some point, because the tester had a hickup yesterday. |
1b8e4b7
to
d1eb0b7
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Very nice! I have three suggestions for improving the speed (of varying usefulness). Please address at least the top and bottom comment, the one about saving the call to std::sqrt
is up to you. The rest is really well done, and will be very useful.
@@ -0,0 +1,163 @@ | |||
/* | |||
Copyright (C) 2017 by the authors of the ASPECT code. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is this our usual indentation or is there typo (double space) here?
for (typename std::multimap<types::LevelInd, Particle<dim> >::const_iterator particle = particle_range.first; | ||
particle != particle_range.second; ++particle) | ||
{ | ||
const double dist = positions[pos_idx].distance(particle->second.get_location()); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you are super concerned about performance you can replace the call to distance
by something like (positions[pos_idx] - particle->second.get_location()).norm_square()
to save one call to std::sqrt
. Not important, I just need to document I have reade the code 😄.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is way less important than my comment above.
const unsigned int n_particles = std::distance(particle_range.first,particle_range.second); | ||
const unsigned int n_particle_properties = particles.begin()->second.get_properties().size(); | ||
|
||
std::vector<std::vector<double> > point_properties; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You already know how large this vector of vectors is going to be, and you can do all of the memory allocation here, instead of inside of the inner loops below the following way:
std::vector<double> temp(n_particle_properties,0.0);
point_properties (n_positions.size(),temp);
Then below you can always fill point_properties[pos_idx][i]
instead of calling push_back()
. This should speed things up significantly (I do not know how fast or slow this algorithm currently is, likely faster than cell average, but no reason to waste efficiency).
nearest_neighbor = particle->second; | ||
} | ||
} | ||
std::vector<double> properties (n_particle_properties,0.0); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You will not need this vector anymore after the changes suggested above
nearest_neighbor_cell = i; | ||
} | ||
} | ||
std::vector<double> properties (n_particle_properties,0.0); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
same here. no need after the changes
} | ||
} | ||
std::vector<double> properties (n_particle_properties,0.0); | ||
const std::vector<double> neighbor_props = properties_at_points(particles, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think you can just assign the return value to point_properties
. It is the job of the recursively called function to select the appropriate components already.
d1eb0b7
to
c1955fc
Compare
Thanks @gassmoeller for the quick (and specific) review! I addressed all the comments, so I think this should be ready to go as long as the tester is happy with it. |
I just remembered, could you write an entry for modules/changes? This is a contribution that definitely deserves an entry over there. |
While you're all at it reviewing patches, here's one more!
This adds a particle interpolator which gets properties from the nearest particle in the same cell. If the cell is empty it gets the nearest particle from the nearest cell. This helps maintain sharp boundaries between compositional fields.