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
|
#!/usr/bin/env python3
#
# Copyright 2021 Remi Lehe
#
# This file is part of WarpX.
#
# License: BSD-3-Clause-LBNL
"""
This script tests the interaction between 2 electron macroparticles.
The macroparticles are moving towards each other and slow down due to the
repelling force between charges of the same sign. The particles should
be moving slow enough for relativistic effects to be negligible. In
this case the conservation of energy gives:
beta(t)^2 = beta(0)^2 + 2 * w * re * log(d(0)/d(t))
where:
re is the classical electron radius
w is the weight of the individual macroparticles
d is the distance between them
beta is the velocity normalized by the speed of light
"""
import numpy as np
from scipy.constants import m_e, c, physical_constants
import sys, re
import yt
import glob
import matplotlib.pyplot as plt
yt.funcs.mylog.setLevel(0)
# Check plotfile name specified in command line
last_filename = sys.argv[1]
filename_radical = re.findall('(.*?)\d+/*$', last_filename)[0]
# Loop through files, and extract the position and velocity of both particles
x1 = []
x2 = []
beta1 = []
beta2 = []
for filename in sorted(glob.glob(filename_radical + '*')):
print(filename)
ds = yt.load(filename)
ad = ds.all_data()
x1.append( float(ad[('electron1','particle_position_x')]) )
x2.append( float(ad[('electron2','particle_position_x')]) )
beta1.append( float(ad[('electron1','particle_momentum_x')])/(m_e*c) )
beta2.append( float(ad[('electron2','particle_momentum_x')])/(m_e*c) )
# Convert to numpy array
x1 = np.array(x1)
x2 = np.array(x2)
beta1 = np.array(beta1)
beta2 = np.array(beta2)
# Plot velocities, compare with theory
w = 5.e12
re = physical_constants['classical electron radius'][0]
beta_th = np.sqrt( beta1[0]**2 - 2*w*re*np.log( (x2[0]-x1[0])/(x2-x1) ) )
plt.plot( beta1, '+', label='Particle 1' )
plt.plot( -beta2, 'x', label='Particle 2' )
plt.plot( beta_th, '*', label='Theory' )
plt.legend(loc=0)
plt.xlabel('Time (a.u.)')
plt.ylabel('Normalized velocity')
plt.savefig('Comparison.png')
# Check that the results are close to the theory
assert np.allclose( beta1[1:], beta_th[1:], atol=0.01 )
assert np.allclose( -beta2[1:], beta_th[1:], atol=0.01 )
# Run checksum regression test
sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
import checksumAPI
test_name = last_filename[:-9]
checksumAPI.evaluate_checksum(test_name, last_filename)
|
