Bilinear Least Squares Particle Interpolation Scheme #1554

hlokavarapu · 2017-05-10T03:15:39Z

@gassmoeller, New bilinear particle interpolation scheme contributed by @yinghe616 , @egpuckett . We plan to incrementally update this scheme in the near future introducing limiters to prevent overshoot and undershoot as well as implement the case for 3D models.

gassmoeller · 2017-05-10T12:40:14Z

/run-tests

gassmoeller · 2017-05-10T12:41:50Z

Just to make sure I understood you correctly: You want to get a review and merge this PR, and later update the scheme? Or is this a work in progress PR?

gassmoeller

Very nice! I have a few comments about functionality and style, but I am looking forward to merging this! I suppose you have tested the accuracy and convergence in your paper? Maybe at some point we can add a benchmark to ASPECT as a follow-up PR, but for now the test is sufficient.

gassmoeller · 2017-05-10T21:18:43Z

source/particle/interpolator/bilinear_least_squares.cc

+            A(index,0) = 1;
+            A(index,1) = (position[0] - approximated_cell_midpoint[0])/cell_diameter;
+            A(index,2) = (position[1] - approximated_cell_midpoint[1])/cell_diameter;
+            A(index,3) = (position[0] - approximated_cell_midpoint[0]) * (position[1] - approximated_cell_midpoint[1])/std::pow(cell_diameter,2);


this will not work, because approximated_cell_midpoint is only set if the cell iterator is empty.

gassmoeller · 2017-05-10T21:19:42Z

source/particle/interpolator/bilinear_least_squares.cc

+                   ExcMessage("The particle property interpolator was not given any "
+                              "positions to evaluate the particle cell_properties at."));
+
+            approximated_cell_midpoint = std::accumulate (positions.begin(), positions.end(), Point<dim>())


See my comment below. This will need to be moved outside of the if-loop. You might as well make it const in that case:

const Point<dim> approximated_cell_midpoint = std:: ...

gassmoeller · 2017-05-10T21:20:55Z

source/particle/interpolator/bilinear_least_squares.cc

+
+        unsigned int property_index = 0;
+        for (unsigned int i=0; i < n_particle_properties; i++)
+          if (selected_properties[i])


So this interpolator only works if only one property is selected in the component_mask? Please check that this is the case at the beginning of the function, currently you silently only compute the last selected component.

gassmoeller · 2017-05-10T21:21:20Z

source/particle/interpolator/bilinear_least_squares.cc

+        r = 0;
+
+        unsigned int index = 0;
+        double cell_diameter = found_cell->diameter();


gassmoeller · 2017-05-10T21:23:31Z

source/particle/interpolator/bilinear_least_squares.cc

+        A.Tmmult(B, A, false);
+        A.Tvmult(c_ATr,r);
+        B_inverse.compute_svd();
+        B_inverse.compute_inverse_svd(threshold);


Do you need both of the calls above? I am not expert in Least squares, and even less in SVD, so please check both calls are necessary and add a comment that explains what is happening in these four lines.

gassmoeller · 2017-05-10T21:23:56Z

source/particle/interpolator/bilinear_least_squares.cc

+
+        for (typename std::vector<Point<dim> >::const_iterator itr = positions.begin(); itr != positions.end(); ++itr, ++index_positions)
+          {
+            Point<dim> support_point = *itr;


gassmoeller · 2017-05-10T21:24:04Z

source/particle/interpolator/bilinear_least_squares.cc

+        for (typename std::vector<Point<dim> >::const_iterator itr = positions.begin(); itr != positions.end(); ++itr, ++index_positions)
+          {
+            Point<dim> support_point = *itr;
+            double interpolated_value = c[0] +


gassmoeller · 2017-05-10T21:25:11Z

source/particle/interpolator/bilinear_least_squares.cc

+    namespace Interpolator
+    {
+      ASPECT_REGISTER_PARTICLE_INTERPOLATOR(BilinearLeastSquares,
+                                            "bilinear",


would you mind changing the input file name to "bilinear least squares"? Somebody might have the idea to add a "bilinear nearest neighbor" or something similar later on.

gassmoeller · 2017-05-10T21:25:40Z

source/particle/interpolator/bilinear_least_squares.cc

+                                            "Interpolates particle properties onto a vector of points using a "
+                                            "bilinear least squares method. Currently only 2D models are "
+                                            "supported. Note that because no limiter is utilized this may "
+                                            "result in overshoot and/or undershoot of interpolated property. ")


remove the last space in the last line

bangerth · 2017-05-12T12:38:21Z

It doesn't currently compile :-(

The error message is this:

/home/dealii/source/source/particle/interpolator/bilinear_least_squares.cc: In instantiation of â€˜std::vector<std::vector<double> > aspect::Particle::Interpolator::BilinearLeastSquares<dim>::properties_at_points(const std::multimap<std::pair<int, int>, aspect::Particle::Particle<dim> >&, const std::vector<dealii::Point<dim> >&, const dealii::ComponentMask&, const typename dealii::parallel::distributed::Triangulation<dim>::active_cell_iterator&) const [with int dim = 2; typename dealii::parallel::distributed::Triangulation<dim>::active_cell_iterator = dealii::TriaActiveIterator<dealii::CellAccessor<2, 2> >]â€™:
/home/dealii/source/source/particle/interpolator/bilinear_least_squares.cc:151:7:   required from here
/home/dealii/source/source/particle/interpolator/bilinear_least_squares.cc:45:21: error: no matching function for call to â€˜dealii::StandardExceptions::ExcMessage::ExcMessage()â€™
                     ExcMessage());
                     ^~~~~~~~~~~~

tjhei · 2017-05-13T12:59:36Z

source/particle/interpolator/bilinear_least_squares.cc

+         */
+        A.Tmmult(B, A, false);
+        A.Tvmult(c_ATr,r);
+        B_inverse.compute_svd();


can you remove this line? compute_inverse_svd will call it automatically.

Sure, removed the line.

tjhei · 2017-05-13T13:02:51Z

source/particle/interpolator/bilinear_least_squares.cc

+                                            "Interpolates particle properties onto a vector of points using a "
+                                            "bilinear least squares method. Currently only 2D models are "
+                                            "supported. Note that because no limiter is utilized this may "
+                                            "result in overshoot and/or undershoot of interpolated property.")


can you add a note that this feature requires deal.II configured with blas/lapack?

tjhei · 2017-05-13T13:03:19Z

tests/particle_interpolator_bilinear_solkz.prm

@@ -0,0 +1,117 @@
+# A description of the SolKZ benchmark using active particles


please update the description here.

tjhei · 2017-05-13T13:05:38Z

source/particle/interpolator/bilinear_least_squares.cc

+         * which usually is not a square matrix. Therefore, we solve Ax=r by
+         * solving A^TAx= A^Tr.
+         */
+        A = 0;


is A initialized to 0 automatically? If yes, remove this.

A is initialized to 0. I checked that this was true by printing out the matrix to std::cout

bangerth · 2017-05-13T18:53:00Z

source/particle/interpolator/bilinear_least_squares.cc

+        dealii::LAPACKFullMatrix<double> A(n_particles, matrix_dimension);
+        Vector<double> r(n_particles);
+
+        /**


we usually use // comments inside functions; but even in cases where one may want to use /* ... */, we never use /** -- these are reserved for doxygen comments that document functions and classes (typically in .h files), not regular commentary.

hlokavarapu · 2017-05-13T20:59:35Z

PR is ready to be reviewed. I have resolved all code review comments thus far.

gassmoeller

2 comments, after that this is ready

gassmoeller · 2017-05-13T21:07:17Z

source/particle/interpolator/bilinear_least_squares.cc

+          {
+            const Point<dim> support_point = *itr;
+            // TODO: Removed const
+            double interpolated_value = c[0] +


please fix the TODO in the line above

gassmoeller · 2017-05-13T21:31:25Z

source/particle/interpolator/bilinear_least_squares.cc

+      {
+        const unsigned int n_particle_properties = particles.begin()->second.get_properties().size();
+
+        unsigned int property_index = selected_properties.first_selected_component(selected_properties.size());


This can be const.
Also does this mean the algorithm only works for one selected component? If so please make sure that is the case by inserting the following line (and check that it compiles):

AssertThrow(selected_properties.n_selected_components(n_particle_properties),ExcNotImplemented());

bangerth · 2017-05-13T23:59:57Z

Can you squash it all into one commit?

hlokavarapu changed the title ~~Particle interpolation scheme bilinear~~ Bilinear Particle Interpolation Scheme May 10, 2017

hlokavarapu changed the title ~~Bilinear Particle Interpolation Scheme~~ Bilinear Least Squares Particle Interpolation Scheme May 10, 2017

hlokavarapu force-pushed the particle_interpolation_scheme_bilinear branch 2 times, most recently from b4a7054 to c9943cb Compare May 10, 2017 12:59

gassmoeller reviewed May 10, 2017

View reviewed changes

hlokavarapu force-pushed the particle_interpolation_scheme_bilinear branch from 118bca2 to 576bea0 Compare May 11, 2017 01:04

hlokavarapu force-pushed the particle_interpolation_scheme_bilinear branch 4 times, most recently from 3280953 to 6bea2d6 Compare May 13, 2017 00:41

tjhei requested changes May 13, 2017

View reviewed changes

bangerth reviewed May 13, 2017

View reviewed changes

gassmoeller reviewed May 13, 2017

View reviewed changes

Added particle interpolation scheme bilinear least squares

36586c0

hlokavarapu force-pushed the particle_interpolation_scheme_bilinear branch from ce8e74c to 36586c0 Compare May 14, 2017 01:14

gassmoeller added the ready to merge label May 14, 2017

gassmoeller merged commit 22db212 into geodynamics:master May 14, 2017

hlokavarapu mentioned this pull request May 18, 2017

Fix bug; Update test case; Update limiting scheme #1749

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Bilinear Least Squares Particle Interpolation Scheme #1554

Bilinear Least Squares Particle Interpolation Scheme #1554

hlokavarapu commented May 10, 2017

gassmoeller commented May 10, 2017

gassmoeller commented May 10, 2017

gassmoeller left a comment

gassmoeller May 10, 2017

gassmoeller May 10, 2017

gassmoeller May 10, 2017

gassmoeller May 10, 2017

gassmoeller May 10, 2017

gassmoeller May 10, 2017

gassmoeller May 10, 2017

gassmoeller May 10, 2017

gassmoeller May 10, 2017

bangerth commented May 12, 2017

tjhei May 13, 2017

hlokavarapu May 13, 2017

tjhei May 13, 2017

tjhei May 13, 2017

tjhei May 13, 2017

hlokavarapu May 13, 2017

bangerth May 13, 2017

hlokavarapu commented May 13, 2017

gassmoeller left a comment

gassmoeller May 13, 2017

gassmoeller May 13, 2017

bangerth commented May 13, 2017

		@@ -0,0 +1,117 @@
		# A description of the SolKZ benchmark using active particles

Bilinear Least Squares Particle Interpolation Scheme #1554

Bilinear Least Squares Particle Interpolation Scheme #1554

Conversation

hlokavarapu commented May 10, 2017

gassmoeller commented May 10, 2017

gassmoeller commented May 10, 2017

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bangerth commented May 12, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hlokavarapu commented May 13, 2017

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bangerth commented May 13, 2017