From 792e91f44177d49cfbe60a4aca62075350f65743 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Sun, 9 Jun 2024 12:48:41 -0400 Subject: [PATCH] Add distance weighted average interpolator --- doc/modules/changes/20240609_gassmoeller | 8 + .../interpolator/distance_weighted_average.h | 72 +++++++++ .../interpolator/distance_weighted_average.cc | 139 ++++++++++++++++++ ...interpolator_distance_weighted_average.prm | 122 +++++++++++++++ .../screen-output | 31 ++++ .../statistics | 26 ++++ 6 files changed, 398 insertions(+) create mode 100644 doc/modules/changes/20240609_gassmoeller create mode 100644 include/aspect/particle/interpolator/distance_weighted_average.h create mode 100644 source/particle/interpolator/distance_weighted_average.cc create mode 100644 tests/particle_interpolator_distance_weighted_average.prm create mode 100644 tests/particle_interpolator_distance_weighted_average/screen-output create mode 100644 tests/particle_interpolator_distance_weighted_average/statistics diff --git a/doc/modules/changes/20240609_gassmoeller b/doc/modules/changes/20240609_gassmoeller new file mode 100644 index 00000000000..f2a3a519d02 --- /dev/null +++ b/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. +
+(Rene Gassmoeller, 2024/06/09) 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..21e75eaa63a --- /dev/null +++ b/include/aspect/particle/interpolator/distance_weighted_average.h @@ -0,0 +1,72 @@ +/* + 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 + . + */ + +#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: + /** + * Update function. Called once at the beginning of each time step. + */ + void update() override; + + /** + * Return the cell-wise evaluated particle properties at the given @p 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; + + private: + /** + * Cached information. + */ + std::vector::active_cell_iterator>> vertex_to_cell_map; + }; + } + } +} + +#endif diff --git a/source/particle/interpolator/distance_weighted_average.cc b/source/particle/interpolator/distance_weighted_average.cc new file mode 100644 index 00000000000..206cffdb23b --- /dev/null +++ b/source/particle/interpolator/distance_weighted_average.cc @@ -0,0 +1,139 @@ +/* + 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 + . + */ + +#include + +#include +#include + +namespace aspect +{ + namespace Particle + { + namespace Interpolator + { + namespace internal + { + double weight (const double distance, const double cell_width) + { + // linear weight function (hat function) + return std::max(1.0 - (distance/cell_width),0.0); + } + } + + + + 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 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> cell_properties(n_interpolate_positions, + std::vector(n_particle_properties, + numbers::signaling_nan())); + + // 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::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 = 0.5 * 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; + } + } + } +} + + +// 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.") + } + } +} diff --git a/tests/particle_interpolator_distance_weighted_average.prm b/tests/particle_interpolator_distance_weighted_average.prm new file mode 100644 index 00000000000..2fccd19853d --- /dev/null +++ b/tests/particle_interpolator_distance_weighted_average.prm @@ -0,0 +1,122 @@ +# 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 Number of particles = 100000 + set Time between data output = 70 + set Data output format = none + 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 + end + end + end +end diff --git a/tests/particle_interpolator_distance_weighted_average/screen-output b/tests/particle_interpolator_distance_weighted_average/screen-output new file mode 100644 index 00000000000..368ef7c0901 --- /dev/null +++ b/tests/particle_interpolator_distance_weighted_average/screen-output @@ -0,0 +1,31 @@ + +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. + Rebuilding Stokes preconditioner... + Solving Stokes system... 34+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... 34+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 + + + diff --git a/tests/particle_interpolator_distance_weighted_average/statistics b/tests/particle_interpolator_distance_weighted_average/statistics new file mode 100644 index 00000000000..4dec4a581b8 --- /dev/null +++ b/tests/particle_interpolator_distance_weighted_average/statistics @@ -0,0 +1,26 @@ +# 1: Time step number +# 2: Time (seconds) +# 3: Time step size (seconds) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of degrees of freedom for all compositions +# 8: Iterations for temperature solver +# 9: Iterations for composition solver 1 +# 10: Iterations for Stokes solver +# 11: Velocity iterations in Stokes preconditioner +# 12: Schur complement iterations in Stokes preconditioner +# 13: RMS velocity (m/s) +# 14: Max. velocity (m/s) +# 15: Minimal value for composition advection_field +# 16: Maximal value for composition advection_field +# 17: Global mass for composition advection_field +# 18: Minimal value for composition advection_particle +# 19: Maximal value for composition advection_particle +# 20: Global mass for composition advection_particle +# 21: Minimal value for composition advection_particle2 +# 22: Maximal value for composition advection_particle2 +# 23: Global mass for composition advection_particle2 +# 24: Number of advected particles +0 0.000000000000e+00 0.000000000000e+00 256 2467 1089 6912 0 0 33 35 102 1.81487746e-04 4.04466712e-04 0.00000000e+00 1.00000000e+00 1.82454287e-01 0.00000000e+00 9.99999998e-01 1.82788208e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00 100000 +1 7.000000000000e+01 7.000000000000e+01 256 2467 1089 6912 0 5 33 35 105 3.27581259e-04 7.32023329e-04 -7.59476082e-03 1.02334847e+00 1.82464209e-01 0.00000000e+00 9.99999998e-01 1.82936297e-01 -2.46643937e-04 3.16607913e-04 -1.01316101e-08 100000