blob: 50fec6282008a6a1cdda41d04415dfac127c1d99 (
plain) (
blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
|
#!/usr/bin/env python3
# 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 os
import sys
import numpy as np
from scipy.special import erf
import yt
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 intrinsic 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 intrinsic 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 intrinsic 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 intrinsic 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 = os.path.split(os.getcwd())[1]
checksumAPI.evaluate_checksum(test_name, fn_final)
|
