-
Notifications
You must be signed in to change notification settings - Fork 5
/
magnetized_cooler.py
369 lines (297 loc) · 14.6 KB
/
magnetized_cooler.py
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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
from __future__ import division
import warp as wp
import numpy as np
import h5py as h5
import os
from warp.data_dumping.openpmd_diag import ParticleDiagnostic
from warp.particles.singleparticle import TraceParticle
from warp.particles.extpart import ZCrossingParticles
# Use more-up-to-date local rswarp
import sys
sys.path.insert(1, '/home/vagrant/jupyter/rswarp')
from rswarp.diagnostics import FieldDiagnostic
from rswarp.cathode.injectors import UserInjectors
from rswarp.utilities.file_utils import cleanupPrevious
import rsoopic.h2crosssections as h2crosssections
sys.path.insert(1, '/home/vagrant/jupyter/rswarp/rswarp/ionization')
from ionization import Ionization
import crosssections as Xsect
# Seed set for reproduction
np.random.seed(967644)
diagDir = 'diags/hdf5'
diagFDir = ['diags/fields/magnetic', 'diags/fields/electric']
# Cleanup command if directories already exist
# Warp's HDF5 diagnostics will not overwrite existing files. This command cleans the diagnostic directory
# to allow for rerunning the simulation
#if wp.comm_world.rank == 0:
# cleanupPrevious(diagDir, diagFDir)
# try:
# os.remove('./diags/crossing_record.h5')
# except OSError:
# pass
####################
# General Parameters
####################
# Set the solver geometry to cylindrical
wp.w3d.solvergeom = wp.w3d.RZgeom
# Switches
particle_diagnostic_switch = True # Record particle data periodically
field_diagnostic_switch = True # Record field/potential data periodically
user_injection = True # Switches injection type
space_charge = True # Controls field solve on/off
simulateIonization = False # Include ionization in simulation
if user_injection:
# User injection thermionic_rz_injector method uses a r**2 scaling to distribute particles uniformly
variable_weight = False
else:
# Warp default injection uses radially weighted particles to make injection uniform in r
variable_weight = True
# Dimensions for the electron cooler
#pipe_radius = 0.1524 / 2. # m (Based on ECE specs)
pipe_radius = 0.03
cooler_length = 30. # m
cathode_temperature = 0.25 # eV
# Beam
beam_beta = 0.990813945176
beam_ke = wp.emass / wp.jperev * wp.clight**2 * (1. / np.sqrt(1-beam_beta**2) - 1.) # eV
#print '*** beam_ke, beam_gamma =', beam_ke, 1. / np.sqrt(1-beam_beta**2)
beam_current = 3. # A
beam_radius = 0.01 # m
##########################
# 3D Simulation Parameters
##########################
# Set cells (nx == ny in cylindrical coordinates)
wp.w3d.nx = 40
wp.w3d.ny = 40
wp.w3d.nz = 40000
# Set field boundary conditions (Warp only allows setting x and y together)
wp.w3d.bound0 = wp.neumann # Use neumann to remove field from electrode edge at boundary
wp.w3d.boundnz = wp.neumann
wp.w3d.boundxy = wp.neumann
# Set boundary dimensions
wp.w3d.xmmin = -pipe_radius
wp.w3d.xmmax = pipe_radius
wp.w3d.ymmin = -pipe_radius
wp.w3d.ymmax = pipe_radius
wp.w3d.zmmin = 0.0
wp.w3d.zmmax = cooler_length
# Set particle boundary conditions
wp.top.pbound0 = wp.absorb
wp.top.pboundnz = wp.absorb
# Type of particle push. 2 includes tan correction to Boris push.
wp.top.ibpush = 2
dx = (wp.w3d.xmmax - wp.w3d.xmmin) / wp.w3d.nx
dz = (wp.w3d.zmmax - wp.w3d.zmmin) / wp.w3d.nz # Because warp doesn't set w3d.dz until after solver instantiated
wp.top.dt = 0.75 * dz / (beam_beta * wp.clight) # Timestep set to cell crossing time with some padding
####################################
# Create Beam and Set its Parameters
####################################
# parameters for electron beam:
# total energy: 3.2676813658 MeV
# beta_e: 0.990813945176
wp.top.lrelativ = True
wp.top.relativity = True
ptcl_per_step = 400 # number of particles to inject on each step
wp.top.npinject = ptcl_per_step
beam = wp.Species(type=wp.Electron, name='Electron', lvariableweights=variable_weight)
if space_charge:
beam.ibeam = beam_current
else:
beam.ibeam = 0.0
wp.w3d.l_inj_exact = True # if true, position and angle of injected particle are
# computed analytically rather than interpolated
"""
A custom injector routine is used here to allow for a very relativistic beam to be injected directly into the simulation
because Warp's built in routine is not based on relativistic kinematics. Setting user_injection = False should work
well for lower energy beams.
"""
if user_injection:
wp.top.inject = 6 # Option 6 is a user specified input distribution
wp.top.ainject = beam_radius # Width of injection area in x
wp.top.binject = beam_radius # Wdith of injection area in y
injector = UserInjectors(beam, wp.w3d, wp.gchange, cathode_temperature=cathode_temperature,
cathode_radius=beam_radius,
ptcl_per_step=ptcl_per_step, accelerating_voltage=beam_ke, zmin_scale=0.545)
wp.installuserinjection(injector.thermionic_rz_injector)
else:
wp.top.inject = 1 # Option 1 is a constant injection (default: particles injected from xy plane at z=0)
wp.top.lhalfmaxwellinject = True # Use half maxwell in axis perpendicular to emission surface (full in other axes)
beam.a0 = beam_radius # Width of injected beam in x
beam.b0 = beam_radius # Width of injected beam in y
beam.ap0 = 0.0 # Width of injected beam in vx/vz
beam.bp0 = 0.0 # Width of injected beam in vy/vz
beam.vbeam = beam_beta * wp.clight
assert beam.vbeam < 0.1 * wp.clight, "Injection velocity > 0.1c. " \
"Constant injection does not use relativistic kinematics"
beam.vthz = np.sqrt(cathode_temperature * wp.jperev / beam.mass) # z thermal velocity
beam.vthperp = np.sqrt(cathode_temperature * wp.jperev / beam.mass) # x and y thermal velocity
wp.w3d.l_inj_rz = (wp.w3d.solvergeom == wp.w3d.RZgeom)
##############################
# Ionization of background gas
##############################
# Ionization e- + H2 -> H2+ + e- + e-
if simulateIonization is True:
# Particle species for emission products of an ionization event
h2plus = wp.Species(type=wp.Dihydrogen, charge_state=+1, name='H2+', weight=1000)
emittedelec = wp.Species(type=wp.Electron, name='emitted e-', weight=1000)
target_pressure = 0.4 # in Pa
target_temp = 273.0 # in K
target_density = target_pressure / wp.boltzmann / target_temp # in 1/m^3
# Instantiate Ionization
# Dimensions control where background gas reservoir will be present in the domain
ioniz = Ionization(
stride=100, # Number of particles allowed to be involved in ionization events per timestep
xmin=wp.w3d.xmmin,
xmax=wp.w3d.xmmax,
ymin=wp.w3d.ymmin,
ymax=wp.w3d.ymmax,
zmin=wp.w3d.zmmin,
zmax=wp.w3d.zmmax,
nx=wp.w3d.nx,
ny=wp.w3d.ny,
nz=wp.w3d.nz,
l_verbose=True
)
# add method used to add possible ionization events
h2xs = Xsect.H2IonizationEvent()
def xswrapper(vi):
return h2xs.getCrossSection(vi)
ioniz.add(
incident_species=beam,
emitted_species=[h2plus, emittedelec], # iterable of species created from ionization
#cross_section=h2crosssections.h2_ioniz_crosssection, # Cross section, can be float or function
cross_section=xswrapper,
#emitted_energy0=[0, h2crosssections.ejectedEnergy], # Energy of each emitted species, can be float or function
# or set to 'thermal' to create ions with a thermal energy spread set by temperature
emitted_energy0=['thermal', h2crosssections.ejectedEnergy], # Energy of each emitted species, can be float or function
emitted_energy_sigma=[0, 0], # Energy spread of emitted species (gives width of gaussian distribution)
temperature=[target_temp, None],
sampleEmittedAngle=h2crosssections.generateAngle,
writeAngleDataDir=False, # Write file recording statistics of angles
writeAnglePeriod=1000, # Period to write angle data, if used
l_remove_incident=False, # Remove incident macroparticles involved in ionization
l_remove_target=False, # Remove target macroparticles (only used if target_species set)
ndens=target_density # Background gas density (if target is reservoir of gas and not macroparticles)
)
####################
# Install Conductors
####################
# Create conductors that will represent electrodes placed inside vacuum vessel
pipe_voltage = 0.0 # Set main pipe section to ground
electrode_voltage = +2e2 # electrodes held at several V relative to main pipe
assert electrode_voltage < beam_ke, "Electrodes potential greater than beam KE."
pipe_radius = pipe_radius
electrode_length = 0.25
electrode_gap = 0.1
pipe_length = cooler_length - 2 * electrode_length - 2 * electrode_gap
z_positions = [0.0]
z_positions.append(z_positions[-1] + electrode_length)
z_positions.append(z_positions[-1] + electrode_gap)
z_positions.append(z_positions[-1] + pipe_length)
z_positions.append(z_positions[-1] + electrode_gap)
z_positions.append(z_positions[-1] + electrode_length)
conductors = []
#bottom_cap = wp.Box(xsize=pipe_radius, ysize=pipe_radius, zsize=dz, voltage=electrode_voltage, zcent=.5*dz)
bottom_cap = wp.Box(xsize=pipe_radius, ysize=pipe_radius, zsize=dz, voltage=pipe_voltage, zcent=.5*dz)
#conductors.append(bottom_cap)
left_electrode = wp.Box(xsize=2.*dx, ysize=2.*dx, zsize=dz, voltage=electrode_voltage, xcent=beam_radius+dx, zcent=2.5*dz)
#conductors.append(left_electrode)
entrance_electrode = wp.ZCylinderOut(voltage=electrode_voltage, radius=pipe_radius,
zlower=z_positions[0], zupper=z_positions[1])
#conductors.append(entrance_electrode)
beam_pipe = wp.ZCylinderOut(voltage=pipe_voltage, radius=pipe_radius, zlower=0.0, zupper=cooler_length)
# zlower=z_positions[2], zupper=z_positions[3])
conductors.append(beam_pipe)
#exit_electrode = wp.ZPlane(z0=cooler_length, zsign=-1., voltage=electrode_voltage)
exit_electrode = wp.ZCylinderOut(voltage=electrode_voltage, radius=pipe_radius,
zlower=z_positions[4], zupper=z_positions[5])
#conductors.append(exit_electrode)
right_electrode = wp.Box(xsize=2.*dx, ysize=2.*dx, zsize=dz, voltage=electrode_voltage, xcent=beam_radius+dx, zcent=cooler_length-2.5*dz)
#conductors.append(right_electrode)
#top_cap = wp.Box(xsize=pipe_radius, ysize=pipe_radius, zsize=dz, voltage=electrode_voltage, zcent=cooler_length-.5*dz)
top_cap = wp.Box(xsize=pipe_radius, ysize=pipe_radius, zsize=dz, voltage=pipe_voltage, zcent=cooler_length-.5*dz)
#conductors.append(top_cap)
##############################
# Install Ideal Solenoid Field
##############################
# electron cooler interaction takes place inside a solenoid
# for time being this is represented as an ideal magnetic field B = (0, 0, bz)
# External B-Field can be added by setting vector components at each cell
bz = np.zeros([wp.w3d.nx, wp.w3d.ny, wp.w3d.nz])
bz[:, :, :] = 1.0 # T
z_start = wp.w3d.zmmin
z_stop = wp.w3d.zmmax
# Add B-Field to simulation
wp.addnewbgrd(z_start, z_stop,
xs=wp.w3d.xmmin, dx=(wp.w3d.xmmax - wp.w3d.xmmin),
ys=wp.w3d.ymmin, dy=(wp.w3d.ymmax - wp.w3d.ymmin),
nx=wp.w3d.nx, ny=wp.w3d.ny, nz=wp.w3d.nz, bz=bz)
########################
# Register Field Solvers
########################
# prevent GIST plotting from starting upon setup
wp.top.lprntpara = False
wp.top.lpsplots = False
if space_charge:
# magnetostatic solver disabled for now, was causing mpi issues and is very slow
# solverB = MagnetostaticMG()
# solverB.mgtol = [0.01] * 3
# registersolver(solverB)
# Add 2D Field Solve, will be cylindrical solver based on setting of w3d.solvergeom
solverE = wp.MultiGrid2D()
wp.registersolver(solverE)
# Conductors must be registered after Field Solver is instantiated if we want them to impact field solve
for cond in conductors:
wp.installconductor(cond)
#print 'BCs: ', solverE.bounds
#sys.exit(0)
# Conductors set as scrapers will remove impacting macroparticles from the simulation
#scraper = wp.ParticleScraper(conductors, lsavecondid=1)
scraper = wp.ParticleScraper(conductors)
######################
# Particle Diagnostics
######################
#wp.top.lsavelostparticles = True
# HDF5 Particle/Field diagnostic options
if particle_diagnostic_switch:
particleperiod = 1000 # Particle diagnostic write frequency
particle_diagnostic_0 = ParticleDiagnostic(period=particleperiod, top=wp.top, w3d=wp.w3d, # Should always be set
# Include data from all existing species in write
species={species.name: species for species in wp.listofallspecies},
# Option for parallel write (if available on system)
comm_world=wp.comm_world, lparallel_output=False,
# `ParticleDiagnostic` automatically append 'hdf5' to path name
write_dir=diagDir[:-5])
wp.installafterstep(particle_diagnostic_0.write) # Write method is installed as an after-step action
if field_diagnostic_switch:
fieldperiod = 1000 # Field diagnostic write frequency
efield_diagnostic_0 = FieldDiagnostic.ElectrostaticFields(solver=solverE, top=wp.top, w3d=wp.w3d,
comm_world=wp.comm_world,
period=fieldperiod)
wp.installafterstep(efield_diagnostic_0.write)
# No B-Field diagnostic since magnetostatic solve turned off
# bfield_diagnostic_0 = FieldDiagnostic.MagnetostaticFields(solver=solverB, top=top, w3d=w3d,
# comm_world=comm_world,
# period=fieldperiod)
# installafterstep(bfield_diagnostic_0.write)
# Crossing Diagnostics
##zcross_l = ZCrossingParticles(zz=z_positions[1], laccumulate=1)
##zcross_r = ZCrossingParticles(zz=z_positions[4], laccumulate=1)
#zcross_l = ZCrossingParticles(zz=5.*dz, laccumulate=1)
#zcross_r = ZCrossingParticles(zz=cooler_length-5.*dz, laccumulate=1)
###########################
# Generate and Run PIC Code
###########################
#electrons_tracked_t0 = wp.Species(type=wp.Electron)
#tracer_count = 50
wp.derivqty() # Set derived beam properties if any are required
wp.package("w3d") # Use w3d solver/geometry package
wp.generate() # Allocate arrays, generate mesh, perform initial field solve
#wp.restart('magnetized_cooler3200000')
Nsteps = 60000
#particle_diagnostic_switch = False
#field_diagnostic_switch = False
while wp.top.it < Nsteps:
wp.step(1000)
# if wp.top.it % 100000 == 0 and wp.comm_world.rank == 0:
# wp.dump()