diff options
Diffstat (limited to 'Examples/Modules')
| -rwxr-xr-x | Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py | 53 | ||||
| -rwxr-xr-x | Examples/Modules/RigidInjection/analysis_rigid_injection_LabFrame.py | 70 | ||||
| -rw-r--r-- | Examples/Modules/RigidInjection/inputs | 88 | ||||
| -rw-r--r-- | Examples/Modules/RigidInjection/inputs.BoostedFrame | 49 | ||||
| -rw-r--r-- | Examples/Modules/RigidInjection/inputs.LabFrame | 41 | ||||
| -rw-r--r-- | Examples/Modules/ionization/inputs.bf.rt | 61 | ||||
| -rw-r--r-- | Examples/Modules/ionization/inputs.rt | 54 | ||||
| -rwxr-xr-x | Examples/Modules/ionization/ionization_analysis.py | 77 |
8 files changed, 405 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 diff --git a/Examples/Modules/ionization/inputs.bf.rt b/Examples/Modules/ionization/inputs.bf.rt new file mode 100644 index 000000000..2fcb2fedc --- /dev/null +++ b/Examples/Modules/ionization/inputs.bf.rt @@ -0,0 +1,61 @@ +max_step = 420 +amr.n_cell = 16 800 +amr.max_grid_size = 64 +amr.blocking_factor = 16 +amr.plot_int = 10000 +geometry.coord_sys = 0 +geometry.is_periodic = 1 0 +geometry.prob_lo = -5.e-6 -40.e-6 +geometry.prob_hi = 5.e-6 0.e-6 +amr.max_level = 0 + +algo.maxwell_fdtd_solver = ckc +warpx.do_pml = 1 +warpx.cfl = .999 +warpx.do_moving_window = 1 +warpx.moving_window_dir = z +warpx.moving_window_v = 1.0 +warpx.gamma_boost = 2. +warpx.boost_direction = z + +particles.nspecies = 2 +particles.species_names = electrons ions + +ions.mass = 2.3428415e-26 +ions.charge = q_e +ions.injection_style = nuniformpercell +ions.num_particles_per_cell_each_dim = 2 2 +ions.zmin = 0. +ions.zmax = 50.e-6 +ions.profile = constant +ions.density = 1. +ions.momentum_distribution_type = constant +ions.do_field_ionization = 1 +ions.ionization_initial_level = 2 +ions.ionization_product_species = electrons +ions.physical_element = N +ions.do_continuous_injection=1 + +electrons.mass = m_e +electrons.charge = -q_e +electrons.injection_style = nuniformpercell +electrons.num_particles_per_cell_each_dim = 2 2 +electrons.zmin = 0. +electrons.zmax = 50.e-6 +electrons.profile = constant +electrons.density = 2. +electrons.momentum_distribution_type = constant +electrons.do_continuous_injection=1 + +lasers.nlasers = 1 +lasers.names = laser1 +laser1.profile = Gaussian +laser1.position = 0. 0. -1.e-6 +laser1.direction = 0. 0. 1. +laser1.polarization = 1. 0. 0. +laser1.e_max = 7.224e12 # a0=1.8 +laser1.profile_waist = 1.e10 +laser1.profile_duration = 26.685e-15 +laser1.profile_t_peak = 60.e-15 +laser1.profile_focal_distance = 0 +laser1.wavelength = 0.8e-6 diff --git a/Examples/Modules/ionization/inputs.rt b/Examples/Modules/ionization/inputs.rt new file mode 100644 index 000000000..e3e4622a6 --- /dev/null +++ b/Examples/Modules/ionization/inputs.rt @@ -0,0 +1,54 @@ +max_step = 1600 +amr.n_cell = 16 800 +amr.max_grid_size = 64 +amr.blocking_factor = 16 +amr.plot_int = 10000 +geometry.coord_sys = 0 +geometry.is_periodic = 1 0 +geometry.prob_lo = -5.e-6 0.e-6 +geometry.prob_hi = 5.e-6 20.e-6 +amr.max_level = 0 + +algo.maxwell_fdtd_solver = ckc +warpx.do_pml = 1 +warpx.cfl = .999 + +particles.nspecies = 2 +particles.species_names = electrons ions + +ions.mass = 2.3428415e-26 +ions.charge = q_e +ions.injection_style = nuniformpercell +ions.num_particles_per_cell_each_dim = 2 1 +ions.zmin = 5.e-6 +ions.zmax = 15.e-6 +ions.profile = constant +ions.density = 1. +ions.momentum_distribution_type = constant +ions.do_field_ionization = 1 +ions.ionization_initial_level = 2 +ions.ionization_product_species = electrons +ions.physical_element = N + +electrons.mass = m_e +electrons.charge = -q_e +electrons.injection_style = nuniformpercell +electrons.num_particles_per_cell_each_dim = 1 1 +electrons.zmin = 1. +electrons.zmax = 15.e-6 +electrons.profile = constant +electrons.density = 2. +electrons.momentum_distribution_type = constant + +lasers.nlasers = 1 +lasers.names = laser1 +laser1.profile = Gaussian +laser1.position = 0. 0. 3.e-6 +laser1.direction = 0. 0. 1. +laser1.polarization = 1. 0. 0. +laser1.e_max = 7.224e12 # a0=1.8 +laser1.profile_waist = 1.e10 +laser1.profile_duration = 26.685e-15 +laser1.profile_t_peak = 60.e-15 +laser1.profile_focal_distance = 0 +laser1.wavelength = 0.8e-6 diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/ionization_analysis.py new file mode 100755 index 000000000..b94541f90 --- /dev/null +++ b/Examples/Modules/ionization/ionization_analysis.py @@ -0,0 +1,77 @@ +#! /usr/bin/env python + +""" +This script tests the result of the ionization module in WarpX. + +Input files inputs.rt and inputs.bf.rt are used to reproduce the test from +Chen, JCP, 2013, figure 2 (in the lab frame and in a boosted frame, +respectively): a plane-wave laser pulse propagates through a +uniform N2+ neutral plasma and further ionizes the Nitrogen atoms. This test +checks that, after the laser went through the plasma, ~32% of Nitrogen +ions are N5+, in agreement with theory from Chen's article. +""" + +import sys +import yt +import numpy as np +import scipy.constants as scc +yt.funcs.mylog.setLevel(0) + +# Open plotfile specified in command line, and get ion's ionization level. +filename = sys.argv[1] +ds = yt.load( filename ) +ad = ds.all_data() +ilev = ad['ions', 'particle_ionization_level'].v + +# Fraction of Nitrogen ions that are N5+. +N5_fraction = ilev[ilev == 5].size/float(ilev.size) + +print("Number of ions: " + str(ilev.size)) +print("Number of N5+ : " + str(ilev[ilev == 5].size)) +print("N5_fraction : " + str(N5_fraction)) + +do_plot = False +if do_plot: + import matplotlib.pyplot as plt + all_data_level_0 = ds.covering_grid(level=0,left_edge=ds.domain_left_edge, + dims=ds.domain_dimensions) + F = all_data_level_0['boxlib', 'Ex'].v.squeeze() + extent = [ ds.domain_left_edge[1], ds.domain_right_edge[1], + ds.domain_left_edge[0], ds.domain_right_edge[0] ] + ad = ds.all_data() + + # Plot ions with ionization levels + species = 'ions'; + xi = ad[species, 'particle_position_x'].v + zi = ad[species, 'particle_position_y'].v + ii = ad[species, 'particle_ionization_level'].v + plt.figure(figsize=(10,10)) + plt.subplot(211) + plt.imshow(np.abs(F), extent=extent, aspect='auto', + cmap='magma', origin='default') + plt.colorbar() + for lev in range(int(np.max(ii)+1)): + select = (ii == lev) + plt.scatter(zi[select],xi[select],s=.2, + label='ionization level: ' + str(lev)) + plt.legend() + plt.title("abs(Ex) (V/m) and ions") + plt.xlabel("z (m)") + plt.ylabel("x (m)") + plt.subplot(212) + plt.imshow(np.abs(F), extent=extent, aspect='auto', + cmap='magma', origin='default') + plt.colorbar() + + # Plot electrons + species = 'electrons'; + if species in [x[0] for x in ds.field_list]: + xe = ad[species, 'particle_position_x'].v + ze = ad[species, 'particle_position_y'].v + plt.scatter(ze,xe,s=.1,c='r',label='electrons') + plt.title("abs(Ex) (V/m) and electrons") + plt.xlabel("z (m)") + plt.ylabel("x (m)") + plt.savefig("image_ionization.pdf", bbox_inches='tight') + +assert ((N5_fraction > 0.30) and (N5_fraction < 0.34)) |