Add new interpolator

benchmarks/rigid_shear/steady-state/run_all_models.sh

@@ -6,13 +6,13 @@
 processes=8
 ASPECT_EXEC="../plugin/aspect"
 
-for stokes_degree in 2 3; do
+for stokes_degree in 2; do
   for discontinuous_pressure in false; do
-    for refinement in 3 4 5 6; do
-      for particles_per_direction in 4 5 6 7 10 15 20 32 45 64 80; do
+    for refinement in 3 4 5; do
+      for particles_per_direction in 7; do
         for generator in "reference cell"; do
-        for interpolator in "cell average" "bilinear least squares"; do
-        for integrator in "rk2" "rk4"; do
+        for interpolator in "distance weighted average" "cell average" "bilinear least squares"; do
+        for integrator in "rk2"; do
 
         echo "subsection Discretization" > current.prm
         echo "  set Stokes velocity polynomial degree = $stokes_degree" >> current.prm

include/aspect/particle/interpolator/distance_weighted_average.h

@@ -0,0 +1,95 @@
+/*
+ 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 LICENSE.  If not see
+ <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef _aspect_particle_interpolator_distance_weighted_average_h
+#define _aspect_particle_interpolator_distance_weighted_average_h
+
+#include <aspect/particle/interpolator/interface.h>
+#include <aspect/simulator_access.h>
+
+#include <deal.II/grid/grid_tools_cache.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 DistanceWeightedAverage : public Interface<dim>, public aspect::SimulatorAccess<dim>
+      {
+        public:
+          void update() override;
+          /**
+           * Return the cell-wise evaluated properties of the bilinear least squares function at the positions.
+           */
+          std::vector<std::vector<double>>
+          properties_at_points(const ParticleHandler<dim> &particle_handler,
+                               const std::vector<Point<dim>> &positions,
+                               const ComponentMask &selected_properties,
+                               const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const override;
+
+          // avoid -Woverloaded-virtual:
+          using Interface<dim>::properties_at_points;
+
+          /**
+           * Declare the parameters this class takes through input files.
+           */
+          static
+          void
+          declare_parameters (ParameterHandler &prm);
+
+          /**
+           * Read the parameters this class declares from the parameter file.
+           */
+          void
+          parse_parameters (ParameterHandler &prm) override;
+
+        private:
+          /**
+           * Variables related to a limiting scheme that prevents overshoot and
+           * undershoot of interpolated particle properties based on global max
+           * and global min for each propery.
+           */
+          bool use_global_valued_limiter;
+
+          /**
+           * For each interpolated particle property, a global max and global
+           * min are stored as elements of vectors.
+           */
+          std::vector<double> global_maximum_particle_properties;
+          std::vector<double> global_minimum_particle_properties;
+
+          /**
+           * Cached information.
+          */
+          std::vector<std::set<typename Triangulation<dim>::active_cell_iterator>> vertex_to_cell_map;
+      };
+    }
+  }
+}
+
+#endif

include/aspect/particle/property/interface.h

@@ -533,6 +533,13 @@ namespace aspect
           void
           initialize ();
 
+          /**
+           * Update function. This function is called once at the
+           * beginning of each timestep.
+           */
+          void
+          update ();
+
           /**
            * Initialization function for particle properties. This function is
            * called once for each of the particles of a particle

include/aspect/particle/world.h

@@ -166,6 +166,11 @@ namespace aspect
          */
         void initialize();
 
+        /**
+         * Update the particle world.
+        */
+       void update();
+
         /**
          * Get the particle property manager for this particle world.
          *

source/particle/interpolator/distance_weighted_average.cc

@@ -0,0 +1,248 @@
+/*
+  Copyright (C) 2017 - 2018 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 LICENSE.  If not see
+ <http://www.gnu.org/licenses/>.
+ */
+
+#include <aspect/particle/interpolator/distance_weighted_average.h>
+#include <aspect/postprocess/particles.h>
+#include <aspect/simulator.h>
+
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/base/signaling_nan.h>
+#include <deal.II/lac/full_matrix.templates.h>
+
+#include <boost/lexical_cast.hpp>
+
+namespace aspect
+{
+  namespace Particle
+  {
+    namespace Interpolator
+    {
+      namespace internal
+      {
+        double weight (const double distance, const double cell_width)
+        {
+          if (false)
+            //return std::max(std::pow(1.0 - std::pow((distance/cell_width),2.0),2.0),0.0);
+            return std::max(1.0 - std::pow((distance/cell_width),2.0),0.0); // quadratic weight function
+          else
+            return std::max(1.0 -(distance/cell_width),0.0); // linear weight function (hat function)
+        }
+      }
+
+
+
+      template <int dim>
+      void
+      DistanceWeightedAverage<dim>::update()
+      {
+        vertex_to_cell_map = GridTools::vertex_to_cell_map(this->get_triangulation());
+      }
+
+
+      template <int dim>
+      std::vector<std::vector<double>>
+      DistanceWeightedAverage<dim>::properties_at_points(const ParticleHandler<dim> &particle_handler,
+                                                         const std::vector<Point<dim>> &positions,
+                                                         const ComponentMask &selected_properties,
+                                                         const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const
+      {
+        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 unsigned int n_interpolate_positions = positions.size();
+        const unsigned int n_particle_properties = particle_handler.n_properties_per_particle();
+
+        std::vector<std::vector<double>> cell_properties(n_interpolate_positions,
+                                                          std::vector<double>(n_particle_properties,
+                                                                              numbers::signaling_nan<double>()));
+
+          for (unsigned int index_positions = 0; index_positions < n_interpolate_positions; ++index_positions)
+            for (unsigned int index_properties = 0; index_properties < n_particle_properties; ++index_properties)
+              if (selected_properties[index_properties])
+                cell_properties[index_positions][index_properties] = 0.0;
+
+        std::set<typename Triangulation<dim>::active_cell_iterator> cell_and_neighbors;
+
+        for (const auto v : cell->vertex_indices())
+          {
+            const unsigned int vertex_index = cell->vertex_index(v);
+            cell_and_neighbors.insert(vertex_to_cell_map[vertex_index].begin(),vertex_to_cell_map[vertex_index].end());
+          }
+
+        const double interpolation_support = found_cell->diameter();
+        std::vector<double> integrated_weight(n_interpolate_positions,0.0);
+
+        for (const auto &current_cell: cell_and_neighbors)
+          {
+            const typename ParticleHandler<dim>::particle_iterator_range particle_range =
+              particle_handler.particles_in_cell(current_cell);
+
+            for (const auto &particle: particle_range)
+              {
+                const ArrayView<const double> particle_properties = particle.get_properties();
+                unsigned int index_positions = 0;
+
+                for (const auto &interpolation_point: positions)
+                  {
+                    const double distance = particle.get_location().distance(interpolation_point);
+                    const double weight = internal::weight(distance,interpolation_support);
+
+                    for (unsigned int index_properties = 0; index_properties < particle_properties.size(); ++index_properties)
+                      if (selected_properties[index_properties])
+                        cell_properties[index_positions][index_properties] += weight * particle_properties[index_properties];
+
+                    integrated_weight[index_positions] += weight;
+                    ++index_positions;
+                  }
+              }
+          }
+
+          for (unsigned int index_positions = 0; index_positions < n_interpolate_positions; ++index_positions)
+          {
+            AssertThrow(integrated_weight[index_positions] > 0.0, ExcInternalError());
+
+            for (unsigned int index_properties = 0; index_properties < n_particle_properties; ++index_properties)
+              if (selected_properties[index_properties])
+                cell_properties[index_positions][index_properties] /= integrated_weight[index_positions];
+          }
+
+        return cell_properties;
+      }
+
+      template <int dim>
+      void
+      DistanceWeightedAverage<dim>::declare_parameters (ParameterHandler &prm)
+      {
+        prm.enter_subsection("Postprocess");
+        {
+          prm.enter_subsection("Particles");
+          {
+            prm.enter_subsection("Interpolator");
+            {
+              prm.enter_subsection("Bilinear least squares");
+              {
+                prm.declare_entry ("Global particle property maximum",
+                                   boost::lexical_cast<std::string>(std::numeric_limits<double>::max()),
+                                   Patterns::List(Patterns::Double ()),
+                                   "The maximum global particle property values that will be used as a "
+                                   "limiter for the bilinear least squares interpolation. The number of the input "
+                                   "'Global particle property maximum' values separated by ',' has to be "
+                                   "the same as the number of particle properties.");
+                prm.declare_entry ("Global particle property minimum",
+                                   boost::lexical_cast<std::string>(-std::numeric_limits<double>::max()),
+                                   Patterns::List(Patterns::Double ()),
+                                   "The minimum global particle property that will be used as a "
+                                   "limiter for the bilinear least squares interpolation. The number of the input "
+                                   "'Global particle property minimum' values separated by ',' has to be "
+                                   "the same as the number of particle properties.");
+                prm.declare_entry("Use limiter", "false",
+                                  Patterns::Bool (),
+                                  "Whether to apply a global particle property limiting scheme to the interpolated "
+                                  "particle properties.");
+
+              }
+              prm.leave_subsection();
+            }
+            prm.leave_subsection();
+          }
+          prm.leave_subsection();
+        }
+        prm.leave_subsection();
+      }
+
+      template <int dim>
+      void
+      DistanceWeightedAverage<dim>::parse_parameters (ParameterHandler &prm)
+      {
+        prm.enter_subsection("Postprocess");
+        {
+          prm.enter_subsection("Particles");
+          {
+            prm.enter_subsection("Interpolator");
+            {
+              prm.enter_subsection("Bilinear least squares");
+              {
+                use_global_valued_limiter = prm.get_bool("Use limiter");
+                if (use_global_valued_limiter)
+                  {
+                    global_maximum_particle_properties = Utilities::string_to_double(Utilities::split_string_list(prm.get("Global particle property maximum")));
+                    global_minimum_particle_properties = Utilities::string_to_double(Utilities::split_string_list(prm.get("Global particle property minimum")));
+
+                    const Postprocess::Particles<dim> &particle_postprocessor =
+                      this->get_postprocess_manager().template get_matching_postprocessor<const Postprocess::Particles<dim>>();
+                    const unsigned int n_property_components = particle_postprocessor.get_particle_world().get_property_manager().get_n_property_components();
+
+                    AssertThrow(global_minimum_particle_properties.size() == n_property_components,
+                                ExcMessage("Make sure that the size of list 'Global minimum particle property' "
+                                           "is equivalent to the number of particle properties."));
+
+                    AssertThrow(global_maximum_particle_properties.size() == n_property_components,
+                                ExcMessage("Make sure that the size of list 'Global maximum particle property' "
+                                           "is equivalent to the number of particle properties."));
+                  }
+              }
+              prm.leave_subsection();
+            }
+            prm.leave_subsection();
+          }
+          prm.leave_subsection();
+        }
+        prm.leave_subsection();
+      }
+    }
+  }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+  namespace Particle
+  {
+    namespace Interpolator
+    {
+      ASPECT_REGISTER_PARTICLE_INTERPOLATOR(DistanceWeightedAverage,
+                                            "distance weighted average",
+                                            "Interpolates particle properties onto a vector of points using a "
+                                            "distance weighed averaging method. The results show similar convergence "
+                                            "properties and dependence on number of PPC as the bilinear least squares method.")
+    }
+  }
+}