Bilinear Least Squares Particle Interpolation Scheme

doc/modules/changes/20170513_egp-cig-reu

@@ -0,0 +1,7 @@
+New particle interpolation scheme bilinear least squares
+using singular value decomposition algorithm. Currently
+only 2D models are supported. We chose a simple overshoot
+and undershoot correction scheme, mainly to truncate based
+on local cell maximum and minimum property value.
+<br>
+(Elbridge Gerry Puckett, Ying He, Harsha Lokavarapu 2017/05/13)

include/aspect/particle/interpolator/bilinear_least_squares.h

@@ -0,0 +1,62 @@
+/*
+ Copyright (C) 2017 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING.  If not see
+ <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef _aspect_particle_interpolator_bilinear_least_squares_h
+#define _aspect_particle_interpolator_bilinear_least_squares_h
+
+#include <aspect/particle/interpolator/interface.h>
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+  namespace Particle
+  {
+    namespace Interpolator
+    {
+      /**
+       * 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
+       */
+      template <int dim>
+      class BilinearLeastSquares : public Interface<dim>, public aspect::SimulatorAccess<dim>
+      {
+        public:
+          /**
+           * Return the cell-wise evaluated properties of the bilinear least squares function at the positions.
+           *
+           * @copydoc aspect::Particle::Interpolator::Interface::properties_at_points()
+           */
+          virtual
+          std::vector<std::vector<double> >
+          properties_at_points(const std::multimap<types::LevelInd, Particle<dim> > &particles,
+                               const std::vector<Point<dim> > &positions,
+                               const ComponentMask &selected_properties,
+                               const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const;
+
+          // avoid -Woverloaded-virtual:
+          using Interface<dim>::properties_at_points;
+      };
+    }
+  }
+}
+
+#endif

source/particle/interpolator/bilinear_least_squares.cc

@@ -0,0 +1,180 @@
+/*
+  Copyright (C) 2017 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING.  If not see
+ <http://www.gnu.org/licenses/>.
+ */
+
+#include <aspect/particle/interpolator/bilinear_least_squares.h>
+
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/base/signaling_nan.h>
+#include <deal.II/lac/full_matrix.templates.h>
+
+namespace aspect
+{
+  namespace Particle
+  {
+    namespace Interpolator
+    {
+      template <int dim>
+      std::vector<std::vector<double> >
+      BilinearLeastSquares<dim>::properties_at_points(const std::multimap<types::LevelInd, Particle<dim> > &particles,
+                                                      const std::vector<Point<dim> > &positions,
+                                                      const ComponentMask &selected_properties,
+                                                      const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const
+      {
+        const unsigned int n_particle_properties = particles.begin()->second.get_properties().size();
+
+        const unsigned int property_index = selected_properties.first_selected_component(selected_properties.size());
+
+        AssertThrow(property_index != numbers::invalid_unsigned_int,
+                    ExcMessage("Internal error: the particle property interpolator was "
+                               "called without a specified component to interpolate."));
+
+        AssertThrow(dim == 2,
+                    ExcMessage("Currently, the particle interpolator 'bilinear' is only supported for 2D models."));
+
+        AssertThrow(selected_properties.n_selected_components(n_particle_properties) == 1,
+                    ExcNotImplemented("Interpolation of multiple components is not supported."));
+
+        const Point<dim> approximated_cell_midpoint = std::accumulate (positions.begin(), positions.end(), Point<dim>())
+                                                      / static_cast<double> (positions.size());
+
+        typename parallel::distributed::Triangulation<dim>::active_cell_iterator found_cell;
+
+        if (cell == typename parallel::distributed::Triangulation<dim>::active_cell_iterator())
+          {
+            // We can not simply use one of the points as input for find_active_cell_around_point
+            // because for vertices of mesh cells we might end up getting ghost_cells as return value
+            // instead of the local active cell. So make sure we are well in the inside of a cell.
+            Assert(positions.size() > 0,
+                   ExcMessage("The particle property interpolator was not given any "
+                              "positions to evaluate the particle cell_properties at."));
+
+
+            found_cell =
+              (GridTools::find_active_cell_around_point<> (this->get_mapping(),
+                                                           this->get_triangulation(),
+                                                           approximated_cell_midpoint)).first;
+          }
+        else
+          found_cell = cell;
+
+        const types::LevelInd cell_index = std::make_pair<unsigned int, unsigned int> (found_cell->level(),found_cell->index());
+        const std::pair<typename std::multimap<types::LevelInd, Particle<dim> >::const_iterator,
+              typename std::multimap<types::LevelInd, Particle<dim> >::const_iterator> particle_range = particles.equal_range(cell_index);
+
+
+        std::vector<std::vector<double> > cell_properties(positions.size(),
+                                                          std::vector<double>(n_particle_properties,
+                                                                              numbers::signaling_nan<double>()));
+
+        const unsigned int n_particles = std::distance(particle_range.first,particle_range.second);
+
+        AssertThrow(n_particles != 0,
+                    ExcMessage("At least one cell contained no particles. The 'bilinear'"
+                               "interpolation scheme does not support this case. "));
+
+
+        // Noticed that the size of matrix A is n_particles x matrix_dimension
+        // which usually is not a square matrix. Therefore, we solve Ax=r by
+        // solving A^TAx= A^Tr.
+        const unsigned int matrix_dimension = 4;
+        dealii::LAPACKFullMatrix<double> A(n_particles, matrix_dimension);
+        Vector<double> r(n_particles);
+        r = 0;
+
+        // Keep track of the local min and max in order to correct
+        // overshoot and undershoot of the interpolated particle property.
+        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];
+
+        unsigned int index = 0;
+        const 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 double particle_property_value = particle->second.get_properties()[property_index];
+            r[index] = 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);
+
+            if (min_value_for_particle_property > particle_property_value)
+              min_value_for_particle_property = particle_property_value;
+            if (max_value_for_particle_property < particle_property_value)
+              max_value_for_particle_property = particle_property_value;
+          }
+
+        dealii::LAPACKFullMatrix<double> B(matrix_dimension, matrix_dimension);
+        dealii::LAPACKFullMatrix<double> B_inverse(B);
+
+        Vector<double> c_ATr(matrix_dimension);
+        Vector<double> c(matrix_dimension);
+
+        const double threshold = 1e-15;
+        unsigned int index_positions = 0;
+
+        // Matrix A can be rank deficient if it does not have full rank, therefore singular.
+        // To circumvent this issue, we solve A^TAx=A^Tr by using singular value
+        // decomposition (SVD).
+        A.Tmmult(B, A, false);
+        A.Tvmult(c_ATr,r);
+        B_inverse.compute_inverse_svd(threshold);
+        B_inverse.vmult(c, c_ATr);
+
+        for (typename std::vector<Point<dim> >::const_iterator itr = positions.begin(); itr != positions.end(); ++itr, ++index_positions)
+          {
+            const Point<dim> support_point = *itr;
+            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;
+
+            // Overshoot and undershoot correction of interpolated particle property.
+            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;
+          }
+        return cell_properties;
+      }
+    }
+  }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+  namespace Particle
+  {
+    namespace Interpolator
+    {
+      ASPECT_REGISTER_PARTICLE_INTERPOLATOR(BilinearLeastSquares,
+                                            "bilinear least squares",
+                                            "Interpolates particle properties onto a vector of points using a "
+                                            "bilinear least squares method. Currently only 2D models are "
+                                            "supported. Note that deal.II must be configured with BLAS/LAPACK.")
+    }
+  }
+}

tests/particle_interpolator_bilinear_solkz.cc

@@ -0,0 +1,2 @@
+#include "../benchmarks/solkz/compositional_fields/solkz_compositional_fields.cc"
+

tests/particle_interpolator_bilinear_solkz.prm

@@ -0,0 +1,116 @@
+# A test based on the SolKz benchmark that performs bilinear
+# least squares interpolation from particle properties to
+# compositional fields that are flagged as 'particle' advected fields.
+
+set Additional shared libraries            = ./libparticle_interpolator_bilinear_solkz.so
+
+############### Global parameters
+
+set Dimension                              = 2
+
+set Start time                             = 0
+set End time                               = 0
+
+set Output directory                       = output
+
+set Pressure normalization                 = volume
+
+############### Parameters describing the model
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 1
+    set Y extent = 1
+  end
+end
+
+subsection Model settings
+  set Prescribed velocity boundary indicators =
+  set Tangential velocity boundary indicators = left, right, bottom, top
+  set Zero velocity boundary indicators       =
+end
+
+subsection Material model
+  set Model name = SolKzCompositionalMaterial
+end
+
+subsection Gravity model
+  set Model name = vertical
+end
+
+############### Parameters describing the temperature field
+
+subsection Boundary temperature model
+  set Model name = box
+end
+
+subsection Initial temperature model
+  set Model name = perturbed box
+end
+
+############### Parameters describing the discretization
+
+subsection Discretization
+  set Stokes velocity polynomial degree       = 2
+  set Use locally conservative discretization = false
+  set Use discontinuous composition discretization = true
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement              = 0
+  set Initial global refinement                = 5
+end
+
+############### Parameters describing the compositional field
+# Note: The compositional field is what drives the flow
+# in this example
+
+subsection Compositional fields
+  set Number of fields = 2
+  set Names of fields = density_comp, viscosity_comp
+  set Compositional field methods = particles, particles
+  set Mapped particle properties = density_comp:function[0], viscosity_comp:function[1]
+end
+
+subsection Initial composition model
+  set Model name = function
+  subsection Function
+    set Variable names      = x,z
+    set Function constants  = pi=3.1415926, eta_b=1e6
+    set Function expression = -1 * sin(2*z) * cos(3*pi*x); exp(log(eta_b)*z)
+  end
+end
+
+############### Parameters describing what to do with the solution
+
+subsection Postprocess
+  set List of postprocessors = particles, SolKzPostprocessor
+
+  subsection Particles
+    set Number of particles = 16384
+    set Time between data output = 0
+    set Data output format = vtu
+    set List of particle properties = function
+    set Integration scheme = rk2
+    set Interpolation scheme = cell average
+    set Maximum particles per cell = 16384
+    set Update ghost particles = true
+
+    subsection Function
+      set Number of components = 2
+      set Variable names      = x, z
+      set Function constants  = pi=3.1415926, eta_b=1e6
+      set Function expression =  -1 * sin(2*z) * cos(3*pi*x); exp(log(eta_b)*z)
+    end
+
+    set Particle generator name = reference cell
+
+    subsection Generator
+      subsection Reference cell
+        set Number of particles per cell per direction = 4
+      end
+    end
+  end
+end

tests/particle_interpolator_bilinear_solkz/screen-output

@@ -0,0 +1,24 @@
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+
+Loading shared library <./libparticle_interpolator_bilinear_solkz.so>
+
+Number of active cells: 1,024 (on 6 levels)
+Number of degrees of freedom: 32,196 (8,450+1,089+4,225+9,216+9,216)
+
+*** Timestep 0:  t=0 years
+   Skipping temperature solve because RHS is zero.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 46+0 iterations.
+
+   Postprocessing:
+     Writing particle output:       output-particle_interpolator_bilinear_solkz/particles/particles-00000
+     Errors u_L1, p_L1, u_L2, p_L2: 2.396693e-07, 3.813761e-04, 3.798509e-07, 1.430329e-03
+
+Termination requested by criterion: end time
+
+
++---------------------------------------------+------------+------------+
++---------------------------------+-----------+------------+------------+
++---------------------------------+-----------+------------+------------+
+