Merge pull request #5815 from gassmoeller/distance_weighted_particle_…

…average

Add distance weighted particle interpolator

doc/modules/changes/20240609_gassmoeller

@@ -0,0 +1,8 @@
+Added: A new particle interpolator plugin 'distance weighted average' that
+computes the interpolated properties as an average of particle properties
+weighted by a distance function. This plugin is similarly accurate
+as the 'cell average' interpolator, but provides smoother results and is
+less susceptible to jumps in interpolated properties if particles move
+small distances.
+<br>
+(Rene Gassmoeller, 2024/06/09)

include/aspect/particle/interpolator/distance_weighted_average.h

@@ -0,0 +1,70 @@
+/*
+ Copyright (C) 2024 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
+       * a distance weighted averaging. The weight function is a hat function.
+       *
+       * @ingroup ParticleInterpolators
+       */
+      template <int dim>
+      class DistanceWeightedAverage : public Interface<dim>, public aspect::SimulatorAccess<dim>
+      {
+        public:
+          /**
+           * Initialize function. Called once at the beginning of the model.
+           */
+          void initialize() override;
+
+          /**
+           * Return the cell-wise evaluated particle properties at the given @p 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;
+
+        private:
+          /**
+           * Cached information that stores information about the grid so that we
+           * do not need to recompute it every time properties_at_points() is called.
+          */
+          std::unique_ptr<GridTools::Cache<dim>> grid_cache;
+      };
+    }
+  }
+}
+
+#endif

source/particle/interpolator/distance_weighted_average.cc

@@ -0,0 +1,151 @@
+/*
+  Copyright (C) 2024 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 <deal.II/base/signaling_nan.h>
+#include <deal.II/grid/grid_tools.h>
+
+namespace aspect
+{
+  namespace Particle
+  {
+    namespace Interpolator
+    {
+      namespace internal
+      {
+        double weight (const double distance, const double interpolation_range)
+        {
+          // linear weight function (hat function)
+          return std::max(1.0 - (distance/interpolation_range),0.0);
+        }
+      }
+
+
+
+      template <int dim>
+      void
+      DistanceWeightedAverage<dim>::initialize()
+      {
+        grid_cache = std::make_unique<GridTools::Cache<dim>>(this->get_triangulation(), this->get_mapping());
+      }
+
+
+
+      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 unsigned int n_interpolate_positions = positions.size();
+        const unsigned int n_particle_properties = particle_handler.n_properties_per_particle();
+
+        // Create with signaling NaNs
+        std::vector<std::vector<double>> cell_properties(n_interpolate_positions,
+                                                          std::vector<double>(n_particle_properties,
+                                                                              numbers::signaling_nan<double>()));
+
+        // Set requested properties to 0.0
+        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;
+
+        const auto &vertex_to_cell_map = grid_cache->get_vertex_to_cell_map();
+
+        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());
+          }
+
+        // Average over all particles that are within half a cell diameter.
+        // This distance strikes a balance between required number of particles per
+        // cell and accuracy of the interpolation. This also assumes we can find
+        // most particles within this distance in neighbor cells. If we do not find
+        // some particles in range (e.g. because the neighbor is refined) this
+        // will slightly affect the accuracy of the interpolation, but we accept that.
+        //
+        // TODO: This could be made dependent on the number of particles per cell, the more
+        // particles, the smaller the interpolation range to increase accuracy.
+        const double interpolation_range = 0.5 * 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_range);
+
+                    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;
+      }
+    }
+  }
+}
+
+
+// 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.")
+    }
+  }
+}

tests/particle_interpolator_distance_weighted_average.prm

@@ -0,0 +1,125 @@
+# A test that performs interpolation from particle properties to
+# compositional fields that are flagged as 'particle' advected
+# fields.
+
+# MPI: 2
+
+set Dimension                              = 2
+set End time                               = 70
+set Use years in output instead of seconds = false
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent  = 0.9142
+    set Y extent  = 1.0000
+  end
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = left, right
+  set Zero velocity boundary indicators       = bottom, top
+end
+
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Reference density             = 1010
+    set Viscosity                     = 1e2
+    set Thermal expansion coefficient = 0
+    set Density differential for compositional field 1 = -10
+  end
+end
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 10
+  end
+end
+
+############### Parameters describing the temperature field
+# Note: The temperature plays no role in this model
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0
+  end
+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 = 3
+  set Names of fields = advection_field, advection_particle, advection_particle2
+  set Compositional field methods = field, particles, particles
+  set Mapped particle properties = advection_particle2:velocity [1], advection_particle:function
+end
+
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Variable names      = x,z
+    set Function constants  = pi=3.1415926
+    set Function expression = 0.5*(1+tanh((0.2+0.02*cos(pi*x/0.9142)-z)/0.02));0.0;0.0
+  end
+end
+
+subsection Discretization
+  set Use discontinuous composition discretization = true
+end
+
+############### Parameters describing the discretization
+
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Strategy                           = composition
+  set Initial global refinement          = 4
+  set Time steps between mesh refinement = 0
+  set Coarsening fraction                = 0.05
+  set Refinement fraction                = 0.3
+end
+
+############### Parameters describing what to do with the solution
+
+subsection Postprocess
+  set List of postprocessors = velocity statistics, composition statistics,particles
+
+  subsection Visualization
+    set Interpolate output = false
+    set Time between graphical output = 70
+  end
+
+  subsection Particles
+    set Time between data output = 70
+    set Data output format = none
+  end
+end
+
+subsection Particles
+  set List of particle properties = velocity, function, initial composition
+  set Interpolation scheme = distance weighted average
+  set Update ghost particles = true
+  set Particle generator name = random uniform
+
+  subsection Function
+    set Variable names      = x,z
+    set Function expression = 0.5*(1+tanh((0.2+0.02*cos(pi*x/0.9142)-z)/0.02))
+  end
+
+  subsection Generator
+    subsection Probability density function
+      set Variable names      = x,z
+      set Function expression = x*x*z
+      set Number of particles = 100000
+    end
+  end
+end

tests/particle_interpolator_distance_weighted_average/screen-output

@@ -0,0 +1,30 @@
+
+Number of active cells: 256 (on 5 levels)
+Number of degrees of freedom: 10,468 (2,178+289+1,089+2,304+2,304+2,304)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Skipping temperature solve because RHS is zero.
+   Advecting particles...  done.
+   Solving advection_field system ... 0 iterations.
+   Solving Stokes system... 12+0 iterations.
+
+   Postprocessing:
+     RMS, max velocity:            0.000181 m/s, 0.000404 m/s
+     Compositions min/max/mass:    0/1/0.1825 // 0/1/0.1828 // 0/0/0
+     Number of advected particles: 100000
+
+*** Timestep 1:  t=70 seconds, dt=70 seconds
+   Skipping temperature solve because RHS is zero.
+   Advecting particles...  done.
+   Solving advection_field system ... 5 iterations.
+   Solving Stokes system... 12+0 iterations.
+
+   Postprocessing:
+     RMS, max velocity:            0.000328 m/s, 0.000732 m/s
+     Compositions min/max/mass:    -0.007595/1.023/0.1825 // 0/1/0.1829 // -0.0002466/0.0003166/-1.013e-08
+     Number of advected particles: 100000
+
+Termination requested by criterion: end time
+
+
+