Particle deletion using point density function

benchmarks/particle_distribution/calculate_removal_algorithms_statistics.py

@@ -0,0 +1,140 @@
+import pandas as pd
+import sys
+import os
+
+sys.path.append(os.path.join(os.path.dirname(__file__), '..', '..'))
+from contrib.python.scripts.aspect_data import read_statistics
+
+# Random selection / no kernel function, oscillating velocity
+#random_df = pd.read_table('output-random/statistics',names=column_names,comment='#',engine='python',sep='\s+')
+random_df = read_statistics('output-random/statistics')
+random_oscillating_data = [
+    'Random, oscillating V',
+    random_df['Cell Score Standard Deviation: '].mean(),
+    random_df['Average particle distribution score: '].mean(),
+    random_df['Maximal particle distribution score: '].mean()
+]
+
+# Gaussian kernel function, oscillating velocity
+gaussian_df = read_statistics('output-gaussian/statistics')
+gaussian_oscillating_data = [
+    'Gaussian, oscillating V',
+    gaussian_df['Cell Score Standard Deviation: '].mean(),
+    gaussian_df['Average particle distribution score: '].mean(),
+    gaussian_df['Maximal particle distribution score: '].mean(),
+]
+
+# Cutoff-w1 kernel function, oscillating velocity
+cutoffw1_df = read_statistics('output-cutoff-w1/statistics')
+cutoffw1_oscillating_data = [
+    'Cutoff-w1, oscillating V',
+    cutoffw1_df['Cell Score Standard Deviation: '].mean(),
+    cutoffw1_df['Average particle distribution score: '].mean(),
+    cutoffw1_df['Maximal particle distribution score: '].mean()
+]
+
+# Cutoff-c1 kernel function, oscillating velocity
+cutoffc1_df = read_statistics('output-cutoff-c1/statistics')
+cutoffc1_oscillating_data = [
+    'Cutoff-c1, oscillating V',
+    cutoffc1_df['Cell Score Standard Deviation: '].mean(),
+    cutoffc1_df['Average particle distribution score: '].mean(),
+    cutoffc1_df['Maximal particle distribution score: '].mean()
+]
+
+# Uniform kernel function, oscillating velocity
+uniform_df = read_statistics('output-uniform/statistics')
+uniform_oscillating_data = [
+    'Uniform, oscillating V',
+    uniform_df['Cell Score Standard Deviation: '].mean(),
+    uniform_df['Average particle distribution score: '].mean(),
+    uniform_df['Maximal particle distribution score: '].mean()
+]
+
+# Triangular kernel function, oscillating velocity
+triangular_df = read_statistics('output-triangular/statistics')
+triangular_oscillating_data = [
+    'Triangular, oscillating V',
+    triangular_df['Cell Score Standard Deviation: '].mean(),
+    triangular_df['Average particle distribution score: '].mean(),
+    triangular_df['Maximal particle distribution score: '].mean()
+]
+
+# Random kernel function, constant velocity
+random_constant_df = read_statistics('output-random-constant-velocity/statistics')
+random_constant_data = [
+    'Random, constant V',
+    random_constant_df['Cell Score Standard Deviation: '].mean(),
+    random_constant_df['Average particle distribution score: '].mean(),
+    random_constant_df['Maximal particle distribution score: '].mean()
+]
+
+# Gaussian kernel function, constant velocity
+gaussian_constant_df = read_statistics('output-gaussian-constant-velocity/statistics')
+gaussian_constant_data = [
+    'Gaussian, constant V',
+    gaussian_constant_df['Cell Score Standard Deviation: '].mean(),
+    gaussian_constant_df['Average particle distribution score: '].mean(),
+    gaussian_constant_df['Maximal particle distribution score: '].mean()
+]
+
+# Cutoff-w1 kernel function, constant velocity
+cutoffw1_constant_df = read_statistics('output-cutoff-w1-constant-velocity/statistics')
+cutoffw1_constant_data = [
+    'Cutoff-w1, constant V',
+    cutoffw1_constant_df['Cell Score Standard Deviation: '].mean(),
+    cutoffw1_constant_df['Average particle distribution score: '].mean(),
+    cutoffw1_constant_df['Maximal particle distribution score: '].mean(),
+]
+
+# Cutoff-c1 kernel function, constant velocity
+cutoffc1_constant_df = read_statistics('output-cutoff-c1-constant-velocity/statistics')
+cutoffc1_constant_data = [
+    'Cutoff-c1, constant V',
+    cutoffc1_constant_df['Cell Score Standard Deviation: '].mean(),
+    cutoffc1_constant_df['Average particle distribution score: '].mean(),
+    cutoffc1_constant_df['Maximal particle distribution score: '].mean()
+]
+
+# Uniform kernel function, constant velocity
+uniform_constant_df = read_statistics('output-uniform-constant-velocity/statistics')
+uniform_constant_data = [
+    'Uniform, constant V',
+    uniform_constant_df['Cell Score Standard Deviation: '].mean(),
+    uniform_constant_df['Average particle distribution score: '].mean(),
+    uniform_constant_df['Maximal particle distribution score: '].mean(),
+]
+
+# Triangular kernel function, constant velocity
+triangular_constant_df = read_statistics('output-triangular-constant-velocity/statistics')
+triangular_constant_data = [
+    'Triangular, constant V',
+    triangular_constant_df['Cell Score Standard Deviation: '].mean(),
+    triangular_constant_df['Average particle distribution score: '].mean(),
+    triangular_constant_df['Maximal particle distribution score: '].mean()
+]
+
+output_data_array = [
+    random_oscillating_data,
+    gaussian_oscillating_data,
+    cutoffw1_oscillating_data,
+    cutoffc1_oscillating_data,
+    uniform_oscillating_data,
+    triangular_oscillating_data,
+    random_constant_data,
+    gaussian_constant_data,
+    cutoffw1_constant_data,
+    cutoffc1_constant_data,
+    uniform_constant_data,
+    triangular_constant_data
+]
+
+column_names_output = [
+    'Particle Removal Algorithm',
+    'Time Averaged Cell Score Standard Deviation',
+    'Time Averaged Mean Score',
+    'Time Averaged Maximum Score',
+]
+
+output_dataframe = pd.DataFrame(output_data_array,columns=column_names_output)
+output_dataframe.to_csv('removal_algorithm_comparison_data.csv',index=False)

benchmarks/particle_distribution/doc/README.md

@@ -0,0 +1,46 @@
+# Particle distribution benchmark
+*This section was contributed by Jarett Baker-Dunn and Rene Gassmöller*
+
+This benchmark is designed to compare different particle removal algorithms. In particular, it is designed to test the use of a point density function (PDF) of particle positions to remove particles from the part of the cell which has the highest density of particles. The benchmark measures how well using a PDF to remove particles functions compared to randomly selecting particles to remove. This comparison matters because randomly deleting particles was the only method available in ASPECT before version 3.1.0, and randomly selecting particles to delete could cause uneven particle distributions to develop in some models where particles were advected across cell refinement boundaries.
+
+The PDF is generated using kernel density estimation. Kernel density estimation is a method for generating PDFs from sets of data using a kernel function and a bandwidth value. The PDFs generated through kernel density estimation are not defined continuously, instead, they are defined at discrete points from which the kernel function has been summed.
+
+At each point at which the PDF is defined, the distance between the point in question and every point in the dataset serves as the input to the kernel function. The output of the kernel function between the point in question and each datapoint is summed and scaled by the bandwidth. The bandwidth affects how closely the PDF reflects the measured data. It is important that the bandwidth is small enough to result in a meaningful PDF, but not so small that the function is overfit to the data.
+
+In ASPECT's particle manager, the kernel function is defined at the location of each particle in the cell, and the particles in the cell serve as the data points on which the kernel function operates. This enables the manager to select the particles with the highest point density values when deleting particles, which works to prevent excessive clustering of particles which can occur with simpler removal algorithms.
+
+The benchmark tests all of the available kernel functions, under constant and oscillating velocity. It also compares the PDF-based algorithm to the first particle removal algorithm implemented, which randomly selected particles to delete. It is important to compare any new particle removal algorithms to random selection because random selection was the only particle removal algorithm available in ASPECT until the PDF based algorithm was added. By comparing the PDF based algorithm to  random selection, this benchmark intends to show that the PDF based algorithm avoids the uneven particle distributions which sometimes resulted from random particle removal.
+
+Models using both constant and oscillating velocity are used in this benchmark. The models prescribing constant velocity make for the simplest possible setup which can reliably generate uneven distributions of particles, which is useful because their simplicity makes it very easy to isolate the effects of different particle removal algorithms. On the other hand models using oscillating velocities more closely resemble the models which are actually used to simulate geodynamic settings. Additionally, oscillating velocity has the advantage that it moves particles across a cell refinement boundary in two directions which can make the test more robust. For these reasons every kernel function and random selection are tested under both velocity schemes.
+
+
+
+The kernel functions tested are:
+
+deal.II's cutoff-w1 function,
+deal.II's cutoff-c1 function,
+a gaussian function,
+a uniform function,
+a triangular function.
+
+The gaussian kernel function emulates a gaussian distribution, so that the value it returns given a certain distance is the value of a gaussian distribution at that distance from the center of the curve.
+
+The uniform kernel function returns a constant value as long as the distance it is given is less than the bandwidth.
+
+The value of the triangular kernel function scales linearly with distance, so that the value returned decreases at a constant rate as the distance input increases. The ratio between increasing distance and decreasing return value is 1:1. If the kernel function were to be graphed, the slope of the triangle's edges would be 1.
+
+The cutoff-w1 kernel function is implemented by the deal.II library. It is similar to the triangular kernel function in that the return value decreases with increasing distance, but it is not linear.
+
+After testing each kernel function and evaluating their performance using the particle distribution score and particles distributions statistics postprocessors, the cutoff-w1 function was chosen as the default function since it performed the best numerically.
+
+
+```{figure-md} fig:plot
+<img src="../ave-score-oscillate.png" alt="Plot of average cell score over time."  width="100%"/>
+
+Figure shows a plot of the average cell score over time. The score of each cell is a measure of how clustered particles are within that cell. The score is calculated using the Particle Distribution Score postprocessor, which uses a histogram based method. A score of 0 represents a cell with a perfectly even distribution of internal particles, while a score of 1 represents a cell with the highest clustering of particles. Lower scores signal a more uniform distribution of particles, which is desirable. In the figure, the random deletion algorithm consistently shows the highest (worst) average score while the PDF algorithm using deal.II kernel functions generally show the lowest (best) score.
+```
+
+
+
+
+The benchmarks work by advecting particles across a cell refinement boundary, causing ASPECT to remove excess particles as particles move from smaller, finer cells into larger, coarser cells. Under constant velocity, each particle crosses the refinement boundary only once and in one direction. Under oscillating velocity, each particles crosses the refinement boundary multiple times and in two directions (as long as the particle is not removed by the particle manager).

benchmarks/particle_distribution/plot_all_removal_algorithms.gnuplot

@@ -0,0 +1,93 @@
+plot "output-cutoff-w1/statistics" using 2:23 title "deal.ii-cutoff-w1", \
+     "output-cutoff-c1/statistics" using 2:23 title "deal.ii-cutoff-c1", \
+     "output-gaussian/statistics" using 2:23 title "gaussian", \
+     "output-triangular/statistics" using 2:23 title "triangular", \
+     "output-uniform/statistics" using 2:23 title "uniform", \
+     "output-random/statistics" using 2:23 title "random"
+set style data linespoints
+set xl "Time"
+set ylabel "Standard Deviation of Cell Scores"
+
+set terminal png
+set output "cell-stdev-mean-oscillate.png"
+replot
+
+
+
+plot "output-cutoff-w1-constant-velocity/statistics" using 2:23 title "deal.ii-cutoff-w1 constant V", \
+     "output-cutoff-c1-constant-velocity/statistics" using 2:23 title "deal.ii-cutoff-c1 constant V", \
+     "output-gaussian-constant-velocity/statistics" using 2:23 title "gaussian constant V", \
+     "output-triangular-constant-velocity/statistics" using 2:23 title "triangular constant V", \
+     "output-uniform-constant-velocity/statistics" using 2:23 title "uniform constant V", \
+     "output-random-constant-velocity/statistics" using 2:23 title "random constant V"
+set style data linespoints
+set xl "Time"
+set ylabel "Standard Deviation of Cell Scores"
+
+set terminal png
+set output "cell-stdev-mean-constant.png"
+replot
+
+
+
+plot "output-cutoff-w1/statistics" using 2:21 title "deal.ii-cutoff-w1", \
+     "output-cutoff-c1/statistics" using 2:21 title "deal.ii-cutoff-c1", \
+     "output-gaussian/statistics" using 2:21 title "gaussian", \
+     "output-triangular/statistics" using 2:21 title "triangular", \
+     "output-uniform/statistics" using 2:21 title "uniform", \
+     "output-random/statistics" using 2:21 title "random"
+set style data linespoints
+set xl "Time"
+set ylabel "Average Cell Score"
+
+set terminal png
+set output "ave-score-oscillate.png"
+replot
+
+
+
+plot "output-cutoff-w1-constant-velocity/statistics" using 2:21 title "deal.ii-cutoff-w1 constant V", \
+     "output-cutoff-c1-constant-velocity/statistics" using 2:21 title "deal.ii-cutoff-c1 constant V", \
+     "output-gaussian-constant-velocity/statistics" using 2:21 title "gaussian constant V", \
+     "output-triangular-constant-velocity/statistics" using 2:21 title "triangular constant V", \
+     "output-uniform-constant-velocity/statistics" using 2:21 title "uniform constant V", \
+     "output-random-constant-velocity/statistics" using 2:21 title "random constant V"
+set style data linespoints
+set xl "Time"
+set ylabel "Average Cell Score"
+
+set terminal png
+set output "ave-score-constant.png"
+replot
+
+
+
+plot "output-cutoff-w1/statistics" using 2:22 title "deal.ii-cutoff-w1", \
+     "output-cutoff-c1/statistics" using 2:22 title "deal.ii-cutoff-c1", \
+     "output-gaussian/statistics" using 2:22 title "gaussian", \
+     "output-triangular/statistics" using 2:22 title "triangular", \
+     "output-uniform/statistics" using 2:22 title "uniform", \
+     "output-random/statistics" using 2:22 title "random"
+set style data linespoints
+set xl "Time"
+set ylabel "Maximum Cell Score"
+
+set terminal png
+set output "max-score-oscillate.png"
+replot
+
+
+
+plot "output-cutoff-w1-constant-velocity/statistics" using 2:22 title "deal.ii-cutoff-w1 constant V", \
+     "output-cutoff-c1-constant-velocity/statistics" using 2:22 title "deal.ii-cutoff-c1 constant V", \
+     "output-gaussian-constant-velocity/statistics" using 2:22 title "gaussian constant V", \
+     "output-triangular-constant-velocity/statistics" using 2:22 title "triangular constant V", \
+     "output-uniform-constant-velocity/statistics" using 2:22 title "uniform constant V", \
+     "output-random-constant-velocity/statistics" using 2:22 title "random constant V"
+set style data linespoints
+set xl "Time"
+set ylabel "Maximum Cell Score"
+
+set terminal png
+set output "max-score-constant.png"
+replot

benchmarks/particle_distribution/removal_algorithm_benchmarks.prm

@@ -0,0 +1,115 @@
+# This parameter file is used as a starting point by the "run_all_deletion_algorithms.sh"
+# script. The script builds on this .prm file to run all particle removal algorithms with all
+# available kernel function.
+
+set Dimension                              = 2
+set Use years instead of seconds           = false
+set End time                               = 6
+set Nonlinear solver scheme                = single Advection, no Stokes
+set Output directory                       = output-default
+set Resume computation                     = false
+set Maximum time step                      = 0.03
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set Y periodic = true
+    set X extent = 0.5
+    set Y extent = 5
+    set Y repetitions = 10
+  end
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Variable names      = x,y
+    set Function constants  =
+    set Function expression = 1
+  end
+end
+
+subsection Boundary temperature model
+  set List of model names = box
+  set Fixed temperature boundary indicators   = left, right
+
+  subsection Box
+    set Left temperature   = 0
+    set Right temperature  = 0
+  end
+end
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 0.0
+  end
+end
+
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Reference density             = 1
+    set Reference specific heat       = 1
+    set Reference temperature         = 0
+    set Thermal conductivity          = 1e-5
+    set Thermal expansion coefficient = 0
+    set Viscosity                     = 0
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement              = 2
+  set Initial global refinement                = 1
+  set Time steps between mesh refinement       = 0
+
+  set Strategy = minimum refinement function
+
+  subsection Minimum refinement function
+    set Coordinate system   = cartesian
+    set Variable names      = x,y
+    set Function expression = if ( y<=4, 2, 4)
+  end
+
+end
+
+subsection Particles
+  set Particle generator name =  reference cell
+  set Maximum particles per cell = 19
+  set Load balancing strategy = remove and add particles
+  set Particle removal algorithm = point density function
+  subsection Generator
+    subsection Reference cell
+      set Number of particles per cell per direction = 4
+    end
+  end
+end
+
+subsection Postprocess
+  set List of postprocessors = velocity statistics, temperature statistics, heat flux statistics, visualization, particles, particle distribution score,particle distribution statistics
+
+  subsection Visualization
+    set Time between graphical output = .001
+    set Interpolate output = false
+  end
+
+  subsection Particles
+    set Time between data output = 0.001
+    set Data output format       = vtu
+  end
+
+end
+
+
+subsection Prescribed Stokes solution
+  set Model name = function
+  subsection Velocity function
+    set Variable names      = x,y,t
+    set Function constants  = velConstant=-0
+    set Function expression = 0; (-0.5*sin(pi*t)) +velConstant
+  end
+end