diff options
Diffstat (limited to 'Examples/Modules/RigidInjection')
5 files changed, 213 insertions, 88 deletions
diff --git a/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py b/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py new file mode 100755 index 000000000..497a30097 --- /dev/null +++ b/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py @@ -0,0 +1,53 @@ +#! /usr/bin/env python + +''' +Analysis script of a WarpX simulation of rigid injection in a boosted frame. + +A Gaussian electron beam starts from -5 microns, propagates rigidly up to +20 microns after which it expands due to emittance only (the focal position is +20 microns). The beam width is measured after ~50 microns, and compared with +the theory (with a 5% error allowed). + +The simulation runs in a boosted frame, and the analysis is done in the lab +frame, i.e., on the back-transformed diagnostics. +''' + +import sys, os, yt, glob +import numpy as np +import scipy.constants as scc +import read_raw_data +yt.funcs.mylog.setLevel(0) + +# filename = sys.argv[1] + +def get_particle_field(snapshot, species, field): + fn = snapshot + '/' + species + files = glob.glob(os.path.join(fn, field + '_*')) + files.sort() + all_data = np.array([]) + for f in files: + data = np.fromfile(f) + all_data = np.concatenate((all_data, data)) + return all_data + +# Read data from back-transformed diagnostics +snapshot = './lab_frame_data/snapshot00001' +header = './lab_frame_data/Header' +allrd, info = read_raw_data.read_lab_snapshot(snapshot, header) +z = np.mean( get_particle_field(snapshot, 'beam', 'z') ) +w = np.std ( get_particle_field(snapshot, 'beam', 'x') ) + +# initial parameters +z0 = 20.e-6 +w0 = 1.e-6 +theta0 = np.arcsin(0.1) + +# Theoretical beam width after propagation if rigid ON +wth = np.sqrt( w0**2 + (z-z0)**2*theta0**2 ) +error = np.abs((w-wth)/wth) + +# Print error and assert small error +print("Beam position: " + str(z)) +print("Beam width : " + str(w)) +print("error: " + str(error)) +assert( error < 0.03 ) diff --git a/Examples/Modules/RigidInjection/analysis_rigid_injection_LabFrame.py b/Examples/Modules/RigidInjection/analysis_rigid_injection_LabFrame.py new file mode 100755 index 000000000..86214ad72 --- /dev/null +++ b/Examples/Modules/RigidInjection/analysis_rigid_injection_LabFrame.py @@ -0,0 +1,70 @@ +#! /usr/bin/env python + +''' +Analysis script of a WarpX simulation of rigid injection. + +A Gaussian electron beam starts from -5 microns, propagates rigidly up to +20 microns after which it expands due to emittance only (the focal position is +20 microns). The beam width is measured after ~50 microns, and compared with +the theory (with a 5% error allowed). + +As a help to the user, the script also compares beam width to the theory in +case rigid injection is OFF (i.e., the beam starts expanding from -5 microns), +in which case a warning is raised. +''' + +import sys +import yt +import numpy as np +import scipy.constants as scc +yt.funcs.mylog.setLevel(0) + +filename = sys.argv[1] + +# WarpX headers include more data when rigid injection is used, +# which gives an error with the last yt release. +# To avoid this issue, the three last lines of WarpXHeader are removed if +# needed. +def remove_rigid_lines(plotfile, nlines_if_rigid): + header_name = plotfile + '/WarpXHeader' + f = open(header_name, 'r') + file_lines = f.readlines() + nlines = len(file_lines) + f.close() + if nlines == nlines_if_rigid: + f = open(header_name, 'w') + f.writelines(file_lines[:-3]) + f.close() + +# Remove rigid injection header lines +remove_rigid_lines(filename, 18) +# Read beam parameters +ds = yt.load( filename ) +ad = ds.all_data() +# Beam longitudinal position +z = np.mean(ad['beam', 'particle_position_y'].v) +# Beam width +w = np.std(ad['beam', 'particle_position_x'].v) + +# initial parameters +z0 = 20.e-6 +z0_no_rigid = -5.e-6 +w0 = 1.e-6 +theta0 = np.arcsin(0.1) + +# Theoretical beam width after propagation if rigid OFF +# Inform the user if rigid injection simply off (just to be kind) +wth_no_rigid = np.sqrt( w0**2 + (z-z0_no_rigid)**2*theta0**2 ) +error_no_rigid = np.abs((w-wth_no_rigid)/wth_no_rigid) +if ( error_no_rigid < 0.05): + print("error no rigid: " + str(error_no_rigid)) + print("Looks like the beam defocuses as if rigid injection were OFF") + +# Theoretical beam width after propagation if rigid ON +wth = np.sqrt( w0**2 + (z-z0)**2*theta0**2 ) +error = np.abs((w-wth)/wth) +# Print error and assert small error +print("Beam position: " + str(z)) +print("Beam width : " + str(w)) +print("error: " + str(error)) +assert( error < 0.05 ) diff --git a/Examples/Modules/RigidInjection/inputs b/Examples/Modules/RigidInjection/inputs deleted file mode 100644 index 3c0d52ed0..000000000 --- a/Examples/Modules/RigidInjection/inputs +++ /dev/null @@ -1,88 +0,0 @@ -# Maximum number of time steps -max_step = 60 - -# number of grid points -amr.n_cell = 32 32 32 - -# Maximum allowable size of each subdomain in the problem domain; -# this is used to decompose the domain for parallel calculations. -amr.max_grid_size = 16 - -# Maximum level in hierarchy (for now must be 0, i.e., one level in total) -amr.max_level = 0 - -amr.plot_int = 1 # How often to write plotfiles. "<= 0" means no plotfiles. - -# Geometry -geometry.coord_sys = 0 # 0: Cartesian -geometry.is_periodic = 1 1 0 # Is periodic? -geometry.prob_lo = -2. -2. -4. # physical domain -geometry.prob_hi = 2. 2. 4. - -# Verbosity -warpx.verbose = 1 - -# Algorithms - -# interpolation -interpolation.nox = 3 -interpolation.noy = 3 -interpolation.noz = 3 - -# CFL -warpx.cfl = 1.0 - -# Information about the particle species -particles.nspecies = 1 -particles.species_names = electrons -particles.rigid_injected_species = electrons - - - -# -# The electron species information -# - -electrons.charge = -q_e -electrons.mass = m_e -electrons.injection_style = "gaussian_beam" -electrons.x_rms = 0.1 -electrons.y_rms = 0.1 -electrons.z_rms = 0.1 -electrons.x_m = 0. -electrons.y_m = 0. -electrons.z_m = -1.5 -electrons.npart = 1000 -electrons.q_tot = -8.010883097437485e-07 - -electrons.profile = "constant" -electrons.density = 1 -electrons.momentum_distribution_type = "gaussian" -electrons.ux_m = 0. -electrons.uy_m = 0. -electrons.uz_m = 3. -electrons.ux_th = 0.01 -electrons.uy_th = 0.01 -electrons.uz_th = 0.01 - -electrons.xmin = -2 -electrons.xmax = 2 -electrons.ymin = -2 -electrons.ymax = 2 -electrons.zmin = -2 -electrons.zmax = 2 - -electrons.zinject_plane = 0. -electrons.projected = true -electrons.focused = false - -warpx.do_pml = 0 - -# Moving window -warpx.do_moving_window = 0 -warpx.moving_window_dir = z -warpx.moving_window_v = 1.0 # in units of the speed of light - -# Boosted frame -warpx.gamma_boost = 1.5 -warpx.boost_direction = z diff --git a/Examples/Modules/RigidInjection/inputs.BoostedFrame b/Examples/Modules/RigidInjection/inputs.BoostedFrame new file mode 100644 index 000000000..c7a60f14f --- /dev/null +++ b/Examples/Modules/RigidInjection/inputs.BoostedFrame @@ -0,0 +1,49 @@ +# stop_time = 1.5e-13 + +warpx.zmax_plasma_to_compute_max_step = 50.e-6 +warpx.gamma_boost = 5. +warpx.boost_direction = z +warpx.do_boosted_frame_diagnostic = 1 +warpx.num_snapshots_lab = 2 +warpx.dt_snapshots_lab = 1.8679589331096515e-13 + +amr.n_cell = 256 512 # 32 64 +amr.max_grid_size = 512 +amr.blocking_factor = 16 +amr.max_level = 0 +amr.plot_int = 100000 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 1 0 # Is periodic? +geometry.prob_lo = -50.e-6 -10.e-6 +geometry.prob_hi = 50.e-6 0.e-6 +warpx.cfl = .999 +warpx.do_pml = 1 +warpx.do_moving_window = 1 +warpx.moving_window_dir = z +warpx.moving_window_v = 1.0 # in units of the speed of light +warpx.serialize_ics = 1 +particles.nspecies = 1 +particles.species_names = beam +particles.rigid_injected_species = beam +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "gaussian_beam" +beam.x_rms = 1.e-6 +beam.y_rms = 1.e-6 +beam.z_rms = .5e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = -5.e-6 +beam.npart = 10000 +beam.q_tot = -1.e-20 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 1000. +beam.ux_th = 100. +beam.uy_th = 100. +beam.uz_th = 0. +beam.zinject_plane = 20.e-6 +beam.rigid_advance = true +beam.projected = true +beam.focused = false diff --git a/Examples/Modules/RigidInjection/inputs.LabFrame b/Examples/Modules/RigidInjection/inputs.LabFrame new file mode 100644 index 000000000..730b46e01 --- /dev/null +++ b/Examples/Modules/RigidInjection/inputs.LabFrame @@ -0,0 +1,41 @@ +stop_time = 1.5e-13 +amr.n_cell = 32 64 +amr.max_grid_size = 256 +amr.blocking_factor = 16 +amr.max_level = 0 +amr.plot_int = 10000 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 1 0 # Is periodic? +geometry.prob_lo = -50.e-6 -10.e-6 +geometry.prob_hi = 50.e-6 0.e-6 +warpx.cfl = .999 +warpx.do_pml = 1 +warpx.do_moving_window = 1 +warpx.moving_window_dir = z +warpx.moving_window_v = 1.0 # in units of the speed of light +warpx.serialize_ics = 1 +particles.nspecies = 1 +particles.species_names = beam +particles.rigid_injected_species = beam +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "gaussian_beam" +beam.x_rms = 1.e-6 +beam.y_rms = 1.e-6 +beam.z_rms = .5e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = -5.e-6 +beam.npart = 2000 +beam.q_tot = -1.e-20 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 1000. +beam.ux_th = 100. +beam.uy_th = 100. +beam.uz_th = 0. +beam.zinject_plane = 20.e-6 +beam.rigid_advance = true +beam.projected = true +beam.focused = false |