diff --git a/benchmarks/rigid_shear/steady-state/run_all_models.sh b/benchmarks/rigid_shear/steady-state/run_all_models.sh
index 0182fd751e7..c739200f209 100755
--- a/benchmarks/rigid_shear/steady-state/run_all_models.sh
+++ b/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
diff --git a/include/aspect/particle/interpolator/distance_weighted_average.h b/include/aspect/particle/interpolator/distance_weighted_average.h
new file mode 100644
index 00000000000..9cee84693df
--- /dev/null
+++ b/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
+ .
+ */
+
+#ifndef _aspect_particle_interpolator_distance_weighted_average_h
+#define _aspect_particle_interpolator_distance_weighted_average_h
+
+#include
+#include
+
+#include
+
+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
+ class DistanceWeightedAverage : public Interface, public aspect::SimulatorAccess
+ {
+ public:
+ void update() override;
+ /**
+ * Return the cell-wise evaluated properties of the bilinear least squares function at the positions.
+ */
+ std::vector>
+ properties_at_points(const ParticleHandler &particle_handler,
+ const std::vector> &positions,
+ const ComponentMask &selected_properties,
+ const typename parallel::distributed::Triangulation::active_cell_iterator &cell) const override;
+
+ // avoid -Woverloaded-virtual:
+ using Interface::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 global_maximum_particle_properties;
+ std::vector global_minimum_particle_properties;
+
+ /**
+ * Cached information.
+ */
+ std::vector::active_cell_iterator>> vertex_to_cell_map;
+ };
+ }
+ }
+}
+
+#endif
diff --git a/include/aspect/particle/property/interface.h b/include/aspect/particle/property/interface.h
index ea92f560f7e..b98ab5a151c 100644
--- a/include/aspect/particle/property/interface.h
+++ b/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
diff --git a/include/aspect/particle/world.h b/include/aspect/particle/world.h
index 42b4e1f0a07..af7db29f67d 100644
--- a/include/aspect/particle/world.h
+++ b/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.
*
diff --git a/source/particle/interpolator/distance_weighted_average.cc b/source/particle/interpolator/distance_weighted_average.cc
new file mode 100644
index 00000000000..ae1bf1ef6a5
--- /dev/null
+++ b/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
+ .
+ */
+
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+
+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
+ void
+ DistanceWeightedAverage::update()
+ {
+ vertex_to_cell_map = GridTools::vertex_to_cell_map(this->get_triangulation());
+ }
+
+
+ template
+ std::vector>
+ DistanceWeightedAverage::properties_at_points(const ParticleHandler &particle_handler,
+ const std::vector> &positions,
+ const ComponentMask &selected_properties,
+ const typename parallel::distributed::Triangulation::active_cell_iterator &cell) const
+ {
+ const Point approximated_cell_midpoint = std::accumulate (positions.begin(), positions.end(), Point())
+ / static_cast (positions.size());
+
+ typename parallel::distributed::Triangulation::active_cell_iterator found_cell;
+
+ if (cell == typename parallel::distributed::Triangulation::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> cell_properties(n_interpolate_positions,
+ std::vector(n_particle_properties,
+ numbers::signaling_nan()));
+
+ 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::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 integrated_weight(n_interpolate_positions,0.0);
+
+ for (const auto ¤t_cell: cell_and_neighbors)
+ {
+ const typename ParticleHandler::particle_iterator_range particle_range =
+ particle_handler.particles_in_cell(current_cell);
+
+ for (const auto &particle: particle_range)
+ {
+ const ArrayView 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
+ void
+ DistanceWeightedAverage::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::numeric_limits::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::numeric_limits::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
+ void
+ DistanceWeightedAverage::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 &particle_postprocessor =
+ this->get_postprocess_manager().template get_matching_postprocessor>();
+ 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.")
+ }
+ }
+}
diff --git a/source/particle/property/interface.cc b/source/particle/property/interface.cc
index f0cce585756..f434d086ebc 100644
--- a/source/particle/property/interface.cc
+++ b/source/particle/property/interface.cc
@@ -347,6 +347,19 @@ namespace aspect
+ template
+ void
+ Manager::update ()
+ {
+ // Update all property plugins
+ for (const auto &p : property_list)
+ {
+ p->update();
+ }
+ }
+
+
+
template
void
Manager::initialize_one_particle (typename ParticleHandler::particle_iterator &particle) const
diff --git a/source/particle/world.cc b/source/particle/world.cc
index 6acba0d9b8e..3dc20c20e05 100644
--- a/source/particle/world.cc
+++ b/source/particle/world.cc
@@ -81,6 +81,20 @@ namespace aspect
connect_to_signals(this->get_signals());
}
+
+
+ template
+ void
+ World::update()
+ {
+ generator->update();
+ integrator->update();
+ interpolator->update();
+ property_manager->update();
+ }
+
+
+
template
const Property::Manager &
World::get_property_manager() const
@@ -1336,6 +1350,7 @@ namespace aspect
if (SimulatorAccess *sim = dynamic_cast*>(interpolator.get()))
sim->initialize_simulator (this->get_simulator());
interpolator->parse_parameters(prm);
+ interpolator->initialize();
}
}
}
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index e2ce9bfd73d..a68a166cf99 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -637,6 +637,9 @@ namespace aspect
if (prescribed_stokes_solution.get())
prescribed_stokes_solution->update();
+ if (particle_world.get() != nullptr)
+ particle_world->update();
+
// do the same for the traction boundary conditions and other things
// that end up in the bilinear form. we update those that end up in
// the constraints object when calling compute_current_constraints()