In [1]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
os.chdir('..')
import proppy as pp
Source¶
In [2]:
nr_particles = 1*10**4
source_pos = np.array([0.0, 0.0, 0.0], dtype=np.float32)
energy = 10**12 # eV
sphere = 10**14 # [m]
def point_source(nr_particles=nr_particles):
return pp.PointSourceIsotropic(energy, source_pos, nr_particles)
def sphere_source(nr_particles=nr_particles):
radius = sphere
source = pp.SphereSourceIsotropic(energy, source_pos, nr_particles, radius)
return source
Propagator¶
In [3]:
step_size = 1.0*10**12
def propagator(nr_steps=1*10**4, step_size=step_size, diff=1.5*10**21):
speed_of_light = 3*10**8 # [m/s]
diffusion_coefficient = diff # [m^2/s]
mfp_iso = 3*diffusion_coefficient/speed_of_light
mfp = np.array([mfp_iso, mfp_iso, mfp_iso], dtype=np.float32) # [m]
prop = pp.IsotropicPropagator(mfp, nr_steps, step_size)
return prop
Observer¶
In [4]:
def time_evolution_observer(nr_steps):
substeps = [False, False, True] # observe only steps (no substeps)
min_step = 1
max_step = nr_steps
nr_obs_steps = 30
observer = pp.TimeEvolutionObserverLog(min_step, max_step, nr_obs_steps, substeps)
return observer
def sphere_observer(deactivate=True):
substeps = [False, False, True] # observe only steps (no substeps)
spheres = [sphere]
return pp.SphericalObserver(substeps, spheres, on_detection_deactivate=deactivate)
Simulation & Analyze setup¶
In [5]:
def run_simulation(source, observer, propagator, file_name):
sim = pp.Simulation()
sim.add_source(source)
sim.add_propagator(propagator)
sim.add_observer(observer)
%time sim.run_simulation()
sim.save_data(file_name)
In [6]:
def analyze(file_name, diff, nr_particles=nr_particles):
plt.figure(figsize=(5,3))
df = pd.read_pickle(file_name+'.pkl')
bins = 20
trajectory_lengths = df['d']
d = trajectory_lengths/10**14
hist, bins = np.histogram(d, bins=bins)
logbins = np.logspace(np.log10(min(d)),np.log10(max(d)),len(bins))
plt.hist(d, bins=logbins, alpha=0.5, label='$\kappa =$ {:.1e}m$^2$/s'.format(diff))
plt.axvline(x=1, color='k', ls='--', label='plasmoid radius')
plt.title('total # particles = {:.0e}'.format(nr_particles))
plt.xlabel('D/{:.0e}m'.format(sphere))
plt.ylabel('# particles')
plt.loglog()
plt.legend()
plt.show()
plt.figure(figsize=(5,3))
bins = 50
trajectory_lengths = df['d']
d = trajectory_lengths/10**14
hist, bins = np.histogram(d, bins=bins)
linbins = np.linspace((min(d)),(max(d)),len(bins))
plt.hist(d, bins=linbins, alpha=0.5, label='$\kappa =$ {:.1e}m$^2$/s'.format(diff))
plt.axvline(x=1, color='k', ls='--', label='plasmoid radius')
plt.title('total # particles = {:.0e}'.format(nr_particles))
plt.xlabel('D/{:.0e}m'.format(sphere))
plt.ylabel('# particles')
plt.legend()
plt.show()
Simulation¶
Point source with time evolution observer¶
All particles start in the center with isotropic directions (i,i,i), with i=-1 or i =1. Observe all particles n times (with log spacing time intervalls). At each observed time, the running diffusion coefficients (1. plot) and the particle distribution (2. and 3. plots) can be plottet. The latter show that particles start in the center and propagate in positive or negative direction (2. plot) and diffuse after some time (plateau in 1. plot) and Gaus-distribution in the 3. plot.
In [7]:
diff = 10**21 #[m^2/s]
In [8]:
file_0 = 'data/data_tut_5_point_time_kappa'
nr_steps = 3*10**3
run_simulation(point_source(nr_particles=10**3), time_evolution_observer(nr_steps=nr_steps), propagator(step_size=step_size, diff=diff, nr_steps=nr_steps), file_0)
df = pd.read_pickle(file_0+'.pkl')
sta = pp.Statistics(df)
df_kappas = sta.plot_diffusion_coefficients(n_points_plateau=5)
steps = sorted(list(set(df['i'])))
print('\n-----------------------------------')
print('particle distribution after '+str(int(steps[0]))+'. step')
sta.plot_distribution('x', steps[0], 50, None)
print('\n-----------------------------------')
print('particle distribution after last step')
sta.plot_distribution('x', steps[-1], 50, None)
df
Observer initialized Propagator initialized start simulation CPU times: user 6.18 s, sys: 14 ms, total: 6.2 s Wall time: 6.21 s saved data file data/data_tut_5_point_time_kappa.pkl init statistics plotting class
diffusion coefficients computed between 9.94e+14m and 3.00e+15m with 5 data points
kappa_{xx}: 9.95e+20 m²/s +- 1.91e+19 m²/s
kappa_{yy}: 9.06e+20 m²/s +- 1.53e+19 m²/s
kappa_{zz}: 9.43e+20 m²/s +- 1.01e+19 m²/s
-----------------------------------
particle distribution after 1. step
----------------------------------- particle distribution after last step
Out[8]:
| id | i | d | x | y | z | phi | pitch_angle | radius | sub_step | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.0 | 1.0 | 1.000000e+12 | -5.773503e+11 | 5.773503e+11 | -5.773503e+11 | 2.228828 | 1.545519 | -1.0 | 2.0 |
| 1 | 0.0 | 2.0 | 2.000000e+12 | -1.154701e+12 | 1.154701e+12 | -1.154701e+12 | 2.228828 | 1.545519 | -1.0 | 2.0 |
| 2 | 0.0 | 3.0 | 3.000000e+12 | -1.732051e+12 | 5.773503e+11 | -1.732051e+12 | 2.228828 | 1.545519 | -1.0 | 2.0 |
| 3 | 0.0 | 5.0 | 5.000000e+12 | -2.886751e+12 | -5.773503e+11 | -2.886751e+12 | 2.228828 | 1.545519 | -1.0 | 2.0 |
| 4 | 0.0 | 6.0 | 6.000000e+12 | -2.309401e+12 | -1.154701e+12 | -3.464101e+12 | 2.228828 | 1.545519 | -1.0 | 2.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 26995 | 999.0 | 994.0 | 9.940000e+14 | 1.558852e+14 | -2.193930e+13 | 9.353090e+13 | 1.567741 | 2.267743 | -1.0 | 2.0 |
| 26996 | 999.0 | 1310.0 | 1.310000e+15 | 6.004442e+13 | 9.237605e+12 | 1.316363e+14 | 1.567741 | 2.267743 | -1.0 | 2.0 |
| 26997 | 999.0 | 1727.0 | 1.727000e+15 | -9.814945e+12 | 3.637307e+13 | 1.195118e+14 | 1.567741 | 2.267743 | -1.0 | 2.0 |
| 26998 | 999.0 | 2276.0 | 2.276000e+15 | -1.223986e+14 | 5.080682e+13 | 1.408740e+14 | 1.567741 | 2.267743 | -1.0 | 2.0 |
| 26999 | 999.0 | 3000.0 | 3.000000e+15 | -1.639681e+14 | -1.847520e+13 | 1.812888e+14 | 1.567741 | 2.267743 | -1.0 | 2.0 |
27000 rows × 10 columns
Spherical source with time evolution observer¶
All particles start at a random position within a 3d sphere with isotropic directions (i,i,i), with i=-1 or i =1. Here, we show the initial particle distribution (after 1. step) and the final distribution. It shifts towards a Gaussian distribution over time. Effect is better visible with higher statistics.
In [9]:
file_1 = 'data/data_tut_5_source_time_distr'
nr_steps=3*10**3
run_simulation(sphere_source(nr_particles=10**3), time_evolution_observer(nr_steps=nr_steps), propagator(step_size=step_size, diff=diff, nr_steps=nr_steps), file_1)
df = pd.read_pickle(file_1+'.pkl')
sta = pp.Statistics(df)
steps = sorted(list(set(df['i'])))
print('\n-----------------------------------')
print('particle distribution after '+str(int(steps[0]))+'. step')
sta.plot_distribution('x', steps[0], 50, None)
print('\n-----------------------------------')
print('particle distribution after last step')
sta.plot_distribution('x', steps[-1], 50, None)
Observer initialized Propagator initialized start simulation CPU times: user 6.52 s, sys: 514 ms, total: 7.03 s Wall time: 7.28 s saved data file data/data_tut_5_source_time_distr.pkl init statistics plotting class ----------------------------------- particle distribution after 1. step
----------------------------------- particle distribution after last step
Point source with spherical observer¶
All particles start in the center with isotropic directions (i,i,i), with i=-1 or i =1. Observe particles only when they cross the spherical observer. Take observed particles out of the simulation, because they are treated as escaped. The plots show ne number of escaping particles per bin as a function of the travelled distance in plasmoid radii.
In [10]:
file_2 = 'data/data_tut_5_point_sphere_hist'
run_simulation(point_source(), sphere_observer(), propagator(step_size=step_size, diff=diff, nr_steps=1*10**4), file_2)
analyze(file_2, diff)
Observer initialized Propagator initialized start simulation CPU times: user 9.36 s, sys: 30.6 ms, total: 9.39 s Wall time: 9.41 s saved data file data/data_tut_5_point_sphere_hist.pkl
Spherical source with spherical observer¶
Same as above but with homogeneouse particle distribution within complete plasmoid.
In [11]:
file_3 = 'data/data_tut_5_sphere_sphere'
run_simulation(sphere_source(), sphere_observer(), propagator(step_size=step_size, diff=diff, nr_steps=1*10**4), file_3)
analyze(file_3, diff)
Observer initialized Propagator initialized start simulation CPU times: user 4.17 s, sys: 402 µs, total: 4.17 s Wall time: 4.17 s saved data file data/data_tut_5_sphere_sphere.pkl
Point source with spherical observer without escape¶
Same as example "Point source with spherical observer" with one change: particles are not dectivated after observation and can therefore be observed several times.
In [12]:
file_4 = 'data/data_tut_5_sphere_sphere_active'
run_simulation(point_source(nr_particles=5), sphere_observer(deactivate=False), propagator(step_size=step_size, diff=diff, nr_steps=1*10**4), file_4)
analyze(file_4, diff, nr_particles=5)
Observer initialized Propagator initialized start simulation CPU times: user 213 ms, sys: 16.1 ms, total: 229 ms Wall time: 228 ms saved data file data/data_tut_5_sphere_sphere_active.pkl
Source comparison¶
Comparison between a point source and a sphere source. Note that particles need trajectory legths that are larger than the plasmoid radius to escape with a cantral point source.
In [13]:
colors = ['r', 'dodgerblue']
labels = ['point source', 'spherical source']
bins = 100
logbins = np.logspace(np.log10(step_size/sphere), np.log10(step_size*nr_steps/sphere),bins)
for i, file in enumerate([file_2, file_3]):
df = pd.read_pickle(file+'.pkl')
trajectory_lengths = df['d']
d = trajectory_lengths/sphere
plt.hist(d, bins=logbins, histtype=u'step', edgecolor=colors[i], linewidth=1., facecolor="None", label=labels[i])
plt.axvline(x=step_size/sphere, color='grey', ls='--', label='step length')
plt.axvline(x=1, color='k', ls='--', label='plasmoid radius')
plt.title('total # particles injected = {:.0e}'.format(nr_particles), size=14)
plt.xlabel('$l_\mathrm{traj}$'+' / {:.0e}m'.format(sphere), size=14)
plt.ylabel('# particles', size=14)
plt.loglog()
plt.legend(loc = "upper left")
#plt.savefig('traj_lengths_comparison_rwpropa.pdf')
plt.show()
The issues with the initial bins for the spherical source result is caused by the fixed step size which can't resolve the small scales shown in the log-scaling of the x-axis. See the same result in a lin-lin-plot below.
In [14]:
colors = ['r', 'dodgerblue']
labels = ['point source', 'spherical source']
bins = 100
linbins = np.linspace(step_size/sphere, (step_size*nr_steps/sphere),bins)
for i, file in enumerate([file_2, file_3]):
df = pd.read_pickle(file+'.pkl')
trajectory_lengths = df['d']
d = trajectory_lengths/sphere
plt.hist(d, bins=linbins, histtype=u'step', edgecolor=colors[i], linewidth=1., facecolor="None", label=labels[i])
plt.axvline(x=step_size/sphere, color='grey', ls='--', label='step length')
plt.axvline(x=1, color='k', ls='--', label='plasmoid radius')
plt.title('total # particles injected = {:.0e}'.format(nr_particles), size=14)
plt.xlabel('$l_\mathrm{traj}$'+' / {:.0e}m'.format(sphere), size=14)
plt.ylabel('# particles', size=14)
plt.legend(loc = "upper right")
plt.show()
In [ ]: