Use singular value decomposition rather than Gauss Jordan to compute …

…the inverse of a singular matrix

doc/modules/to-1.5.0.h

@@ -7,13 +7,6 @@
  *
  * <ol>
  *
- * <li> New: Support for two particle interpolation schemes (bilinear, and
- * biquadratic) has been added. Using the least squares method, a bilinear
- * or biquadratic function is approximated and evaluated onto compositional
- * fields with a simple overshoot and undershoot correction.
- * <br>
- * (Harsha Lokavarapu, Ying He, Gerry Puckett, 2017/02/27)
- *
  * <li> New: Blankenbach benchmarks were added.
  * <br>
  * (Timo Heister, 2017/02/16)
@@ -87,20 +80,8 @@
  * Added tests and benchmarks for the new formulation, and updated
  * the according particle property in the same way. Also updated
  * and extended the cookbook description in the manual.
-<<<<<<< HEAD:doc/modules/to-1.5.0.h
  * <br>
  * (Rene Gassmoeller, 2017/01/19)
-=======
- * <br> (Rene Gassmoeller, 2017/01/19)
-<<<<<<< HEAD:doc/modules/to-1.5.0.h
-=======
- * fields with a simple overshoot and undershoot correction. Support for
- * 2D models only, has been implemented.
- * <br> (Harsha Lokavarapu, Ying He, Gerry Puckett, 2016/12/28)
->>>>>>> Updated code in response to PR #1333 comments
->>>>>>> 7808038... Updated code in response to PR #1333 comments:doc/modules/changes.h
-=======
->>>>>>> 50e1602... Update particle_at_properties function call for bilinear and biquadratic particle interpolation schemes:doc/modules/changes.h
  *
  * <li> New: Added a "heat flux densities" postprocessor.
  * <br>

include/aspect/particle/interpolator/bilinear_least_squares.h

@@ -1,5 +1,5 @@
 /*
- Copyright (C) 2016 by the authors of the ASPECT code.
+ Copyright (C) 2017 by the authors of the ASPECT code.
 
  This file is part of ASPECT.
 
@@ -31,7 +31,7 @@ namespace aspect
     namespace Interpolator
     {
       /**
-       * Return the interpolated properties of all tracers of the given cell using bilinear least squares method.
+       * Return the interpolated properties of all particles of the given cell using bilinear least squares method.
        * Currently, only the two dimensional model is supported.
        *
        * @ingroup ParticleInterpolators

source/particle/interpolator/bilinear_least_squares.cc

@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2016 by the authors of the ASPECT code.
+  Copyright (C) 2017 by the authors of the ASPECT code.
 
  This file is part of ASPECT.
 
@@ -24,8 +24,6 @@
 #include <deal.II/base/signaling_nan.h>
 #include <deal.II/lac/full_matrix.templates.h>
 
-#include <fstream>
-
 namespace aspect
 {
   namespace Particle
@@ -73,8 +71,8 @@ namespace aspect
         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> > cell_properties(positions.size(), std::vector<double>(n_particle_properties,
-                                                          numbers::signaling_nan<double>()));
+        std::vector<std::vector<double> > cell_properties(positions.size(),
+                                                          std::vector<double>(n_particle_properties, numbers::signaling_nan<double>()));
 
         unsigned int property_index = 0;
         for (unsigned int i=0; i < n_particle_properties; i++)
@@ -85,88 +83,50 @@ namespace aspect
                     ExcMessage("At least one cell contained no particles. The 'bilinear'"
                                "interpolation scheme does not support this case. "));
 
+
         const unsigned int matrix_dimension = 4;
-        dealii::FullMatrix<double> A(n_particles, matrix_dimension);
+        dealii::LAPACKFullMatrix<double> A(n_particles, matrix_dimension);
+        Vector<double> r(n_particles);
         A = 0;
+        r = 0;
 
         unsigned int index = 0;
+        double cell_diameter = found_cell->diameter();
         for (typename std::multimap<types::LevelInd, Particle<dim> >::const_iterator particle = particle_range.first;
              particle != particle_range.second; ++particle, ++index)
-          {
-            const Point<dim> position = particle->second.get_location();
-            A(index,0) = 1;
-            A(index,1) = (position[0] - approximated_cell_midpoint[0])/found_cell->diameter();
-            A(index,2) = (position[1] - approximated_cell_midpoint[1])/found_cell->diameter();
-            A(index,3) = (position[0] - approximated_cell_midpoint[0]) * (position[1] - approximated_cell_midpoint[1])/std::pow(found_cell->diameter(),2);
-          }
-
-        dealii::FullMatrix<double> B(matrix_dimension, matrix_dimension);
-        A.Tmmult(B, A, false);
-        dealii::FullMatrix<double> B_inverse(B);
-
-          std::ofstream debug;
-          debug.open("singular_matrix.out", std::ofstream::app);
-
-          debug << "===========Matrix A============" << std::endl;
-          B.print_formatted(debug, 16, true, 0, "0", 1, 0);
-
-
-
-        dealii::FullMatrix<double> r(matrix_dimension,1);
-        r = 0;
-
-        double max_value_for_particle_property = (particle_range.first)->second.get_properties()[property_index];
-        double min_value_for_particle_property = (particle_range.first)->second.get_properties()[property_index];
-
-        for (typename std::multimap<types::LevelInd, Particle<dim> >::const_iterator particle = particle_range.first;
-             particle != particle_range.second; ++particle)
           {
             const double particle_property_value = particle->second.get_properties()[property_index];
-            const Point<dim> position = particle->second.get_location();
-
-            r(0,0) += particle_property_value;
-            r(1,0) += particle_property_value * (position[0] - approximated_cell_midpoint[0])/found_cell->diameter();
-            r(2,0) += particle_property_value * (position[1] - approximated_cell_midpoint[1])/found_cell->diameter();
-            r(3,0) += particle_property_value * (position[0] - approximated_cell_midpoint[0]) * (position[1] - approximated_cell_midpoint[1])/std::pow(found_cell->diameter(),2);
+            r[index] = particle_property_value;
 
-            if (max_value_for_particle_property < particle_property_value)
-              max_value_for_particle_property = particle_property_value;
-            if (min_value_for_particle_property > particle_property_value)
-              min_value_for_particle_property = particle_property_value;
+            const Point<dim> position = particle->second.get_location();
+            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);
           }
 
+        dealii::LAPACKFullMatrix<double> B(matrix_dimension, matrix_dimension);
+        dealii::LAPACKFullMatrix<double> B_inverse(B);
 
-          debug << "===========Matrix f============" << std::endl;
-          r.print_formatted(debug, 16, true, 0, "0", 1, 0);
+        Vector<double> c_ATr(matrix_dimension);
+        Vector<double> c(matrix_dimension);
 
-          debug.close();
+        const double threshold = 1e-15;
+        unsigned int index_positions = 0;
 
-        dealii::FullMatrix<double> c(matrix_dimension,1);
-        c = 0;
-          try
-          {
-              B_inverse.gauss_jordan();
-          }
-          catch (...)
-          {
-              std::cout << "Level:" << found_cell->level() << "index:" << found_cell->index() << std::endl;
-          }
-        B_inverse.mmult(c, r);
+        A.Tmmult(B, A, false);
+        A.Tvmult(c_ATr,r);
+        B_inverse.compute_svd();
+        B_inverse.compute_inverse_svd(threshold);
+        B_inverse.vmult(c, c_ATr);
 
-        unsigned int index_positions = 0;
         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,0) +
-                                        c(1,0)*(support_point[0] - approximated_cell_midpoint[0])/found_cell->diameter() +
-                                        c(2,0)*(support_point[1] - approximated_cell_midpoint[1])/found_cell->diameter() +
-                                        c(3,0)*(support_point[0] - approximated_cell_midpoint[0])*(support_point[1] - approximated_cell_midpoint[1])/std::pow(found_cell->diameter(),2);
-
-            if (interpolated_value > max_value_for_particle_property)
-              interpolated_value = max_value_for_particle_property;
-            else if (interpolated_value < min_value_for_particle_property)
-              interpolated_value = min_value_for_particle_property;
-
+            double interpolated_value = c[0] +
+                                        c[1]*(support_point[0] - approximated_cell_midpoint[0])/cell_diameter +
+                                        c[2]*(support_point[1] - approximated_cell_midpoint[1])/cell_diameter +
+                                        c[3]*(support_point[0] - approximated_cell_midpoint[0])*(support_point[1] - approximated_cell_midpoint[1])/std::pow(cell_diameter,2);
             cell_properties[index_positions][property_index] = interpolated_value;
           }
         return cell_properties;
@@ -187,7 +147,8 @@ namespace aspect
                                             "bilinear",
                                             "Interpolates particle properties onto a vector of points using a "
                                             "bilinear least squares method. Currently only 2D models are "
-                                            "supported. ")
+                                            "supported. Note that because no limiter is utilized this may "
+                                            "result in overshoot and/or undershoot of interpolated property. ")
     }
   }
 }