2 files changed, 203 insertions, 0 deletions
diff --git a/Examples/Modules/resampling/analysis_leveling_thinning.py b/Examples/Modules/resampling/analysis_leveling_thinning.py
new file mode 100755
index 000000000..dc96b0138
--- /dev/null
+++ b/Examples/Modules/resampling/analysis_leveling_thinning.py
@@ -0,0 +1,138 @@
+#! /usr/bin/env python
+
+# Copyright 2020 Neil Zaim
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+## In this test, we check that leveling thinning works as expected on two simple cases. Each case
+## corresponds to a different particle species.
+
+import yt
+import numpy as np
+import sys
+from scipy.special import erf
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+fn_final = sys.argv[1]
+fn0 = fn_final[:-4] + '0000'
+
+ds0 = yt.load(fn0)
+ds = yt.load(fn_final)
+
+ad0 = ds0.all_data()
+ad = ds.all_data()
+
+numcells = 16*16
+t_r = 1.3 # target ratio
+relative_tol = 1.e-13 # tolerance for machine precision errors
+
+
+#### Tests for first species ####
+# Particles are present in all simulation cells and all have the same weight
+
+ppc = 400.
+numparts_init = numcells*ppc
+
+w0 = ad0['resampled_part1','particle_weight'].to_ndarray() # weights before resampling
+w = ad['resampled_part1','particle_weight'].to_ndarray() # weights after resampling
+
+# Renormalize weights for easier calculations
+w0 = w0*ppc
+w = w*ppc
+
+# Check that initial number of particles is indeed as expected
+assert(w0.shape[0] == numparts_init)
+# Check that all initial particles have the same weight
+assert(np.unique(w0).shape[0] == 1)
+# Check that this weight is 1 (to machine precision)
+assert(abs(w0[0] - 1) < relative_tol)
+
+# Check that the number of particles after resampling is as expected
+numparts_final = w.shape[0]
+expected_numparts_final = numparts_init/t_r**2
+error = np.abs(numparts_final - expected_numparts_final)
+std_numparts_final = np.sqrt(numparts_init/t_r**2*(1.-1./t_r**2))
+# 5 sigma test that has an intrisic probability to fail of 1 over ~2 millions
+print("difference between expected and actual final number of particles (1st species): " + str(error))
+print("tolerance: " + str(5*std_numparts_final))
+assert(error<5*std_numparts_final)
+
+# Check that the final weight is the same for all particles (to machine precision) and is as
+# expected
+final_weight = t_r**2
+assert(np.amax(np.abs(w-final_weight)) < relative_tol*final_weight )
+
+
+#### Tests for second species ####
+# Particles are only present in a single cell and have a gaussian weight distribution
+# Using a single cell makes the analysis easier because leveling thinning is done separately in
+# each cell
+
+ppc = 100000.
+numparts_init = ppc
+
+w0 = ad0['resampled_part2','particle_weight'].to_ndarray() # weights before resampling
+w = ad['resampled_part2','particle_weight'].to_ndarray() # weights after resampling
+
+# Renormalize and sort weights for easier calculations
+w0 = np.sort(w0)*ppc
+w = np.sort(w)*ppc
+
+## First we verify that the initial distribution is as expected
+
+# Check that the mean initial weight is as expected
+mean_initial_weight = np.average(w0)
+expected_mean_initial_weight = 2.*np.sqrt(2.)
+error = np.abs(mean_initial_weight - expected_mean_initial_weight)
+expected_std_initial_weight = 1./np.sqrt(2.)
+std_mean_initial_weight = expected_std_initial_weight/np.sqrt(numparts_init)
+# 5 sigma test that has an intrisic probability to fail of 1 over ~2 millions
+print("difference between expected and actual mean initial weight (2nd species): " + str(error))
+print("tolerance: " + str(5*std_mean_initial_weight))
+assert(error<5*std_mean_initial_weight)
+
+# Check that the initial weight variance is as expected
+variance_initial_weight = np.var(w0)
+expected_variance_initial_weight = 0.5
+error = np.abs(variance_initial_weight - expected_variance_initial_weight)
+std_variance_initial_weight = expected_variance_initial_weight*np.sqrt(2./numparts_init)
+# 5 sigma test that has an intrisic probability to fail of 1 over ~2 millions
+print("difference between expected and actual variance of initial weight (2nd species): " + str(error))
+print("tolerance: " + str(5*std_variance_initial_weight))
+
+## Next we verify that the resampling worked as expected
+
+# Check that the level weight value is as expected from the initial distribution
+level_weight = w[0]
+assert(np.abs(level_weight - mean_initial_weight*t_r) < level_weight*relative_tol)
+
+# Check that the number of particles at the level weight is the same as predicted from analytic
+# calculations
+numparts_leveled = np.argmax(w > level_weight) # This returns the first index for which
+# w > level_weight, which thus corresponds to the number of particles at the level weight
+expected_numparts_leveled = numparts_init/(2.*t_r)*(1+erf(expected_mean_initial_weight*(t_r-1.)) \
+                            -1./(np.sqrt(np.pi)*expected_mean_initial_weight)* \
+                            np.exp(-(expected_mean_initial_weight*(t_r-1.))**2))
+error = np.abs(numparts_leveled - expected_numparts_leveled)
+std_numparts_leveled = np.sqrt(expected_numparts_leveled - numparts_init/np.sqrt(np.pi)/(t_r* \
+                       expected_mean_initial_weight)**2*(np.sqrt(np.pi)/4.* \
+                       (2.*expected_mean_initial_weight**2+1.)*(1.-erf(expected_mean_initial_weight* \
+                       (t_r-1.)))-0.5*np.exp(-(expected_mean_initial_weight*(t_r-1.))**2* \
+                       (expected_mean_initial_weight*(t_r+1.)))))
+# 5 sigma test that has an intrisic probability to fail of 1 over ~2 millions
+print("difference between expected and actual number of leveled particles (2nd species): " + str(error))
+print("tolerance: " + str(5*std_numparts_leveled))
+
+numparts_unaffected = w.shape[0] - numparts_leveled
+numparts_unaffected_anticipated = w0.shape[0] - np.argmax(w0 > level_weight)
+# Check that number of particles with weight higher than level weight is the same before and after
+# resampling
+assert(numparts_unaffected == numparts_unaffected_anticipated)
+# Check that particles with weight higher than level weight are unaffected by resampling.
+assert(np.all(w[-numparts_unaffected:] == w0[-numparts_unaffected:]))
+
+test_name = fn_final[:-9] # Could also be os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, fn_final)
diff --git a/Examples/Modules/resampling/inputs_leveling_thinning b/Examples/Modules/resampling/inputs_leveling_thinning
new file mode 100644
index 000000000..c761a75dd
--- /dev/null
+++ b/Examples/Modules/resampling/inputs_leveling_thinning
@@ -0,0 +1,65 @@
+max_step = 8
+amr.n_cell = 16 16
+amr.blocking_factor = 8
+amr.max_grid_size = 8
+geometry.is_periodic = 1   1
+geometry.prob_lo     = 0.  0.
+geometry.prob_hi     = 16. 16.
+amr.max_level = 0
+
+warpx.do_pml = 0
+
+particles.species_names = resampled_part1 resampled_part2
+
+my_constants.pi = 3.141592653589793
+
+# First particle species. The particles are distributed throughout the simulation box and all have
+# the same weight.
+resampled_part1.species_type = electron
+resampled_part1.injection_style = NRandomPerCell
+resampled_part1.num_particles_per_cell = 400
+resampled_part1.profile = constant
+resampled_part1.density = 1.
+resampled_part1.momentum_distribution_type = constant
+resampled_part1.do_not_deposit = 1
+resampled_part1.do_not_gather = 1
+resampled_part1.do_not_push = 1
+resampled_part1.do_resampling = 1
+resampled_part1.resampling_algorithm = leveling_thinning
+resampled_part1.resampling_algorithm_target_ratio = 1.3
+# This should trigger resampling at timestep 5 only in this case
+resampled_part1.resampling_trigger_intervals = 5
+# This should trigger resampling at first timestep
+resampled_part1.resampling_trigger_max_avg_ppc = 395
+
+# Second particle species. The particles are only present in one cell and have a Gaussian weight
+# distribution.
+resampled_part2.species_type = electron
+resampled_part2.injection_style = NRandomPerCell
+resampled_part2.num_particles_per_cell = 100000
+resampled_part2.zmin = 0.
+resampled_part2.zmax = 1.
+resampled_part2.xmin = 0.
+resampled_part2.xmax = 1.
+resampled_part2.profile = parse_density_function
+# Trick to get a Gaussian weight distribution is to do a Box-Muller transform using the position
+# within the cell as the two random variables. Here, we have a distribution with standard deviation
+# sigma = 1/sqrt(2) and mean 4*sigma.
+resampled_part2.density_function(x,y,z) = 2.*sqrt(2.)+sqrt(-log(x))*cos(2*pi*z)
+resampled_part2.momentum_distribution_type = constant
+resampled_part2.do_not_deposit = 1
+resampled_part2.do_not_gather = 1
+resampled_part2.do_not_push = 1
+resampled_part2.do_resampling = 1
+resampled_part2.resampling_algorithm = leveling_thinning
+resampled_part2.resampling_algorithm_target_ratio = 1.3
+# This should trigger resampling at timestep 7 only in this case. The rest is here to test the
+# intervals parser syntax.
+resampled_part2.resampling_trigger_intervals = 100 :120:3, 7:7:7
+# This should not trigger resampling: this species only has ~391 ppc on average.
+resampled_part2.resampling_trigger_max_avg_ppc = 395
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.period = 8
+diag1.diag_type = Full